Triplets of supermassive black holes: Astrophysics, Gravitational Waves and Detection
Pau Amaro-Seoane, Alberto Sesana, Loren Hoffman, Matthew Benacquista, Christoph Eichhorn, Junichiro Makino, Rainer Spurzem
aa r X i v : . [ a s t r o - ph . C O ] N ov Mon. Not. R. Astron. Soc. , 1–15 (31 October 2018, submitted to MNRAS) Printed 31 October 2018 (MN L A TEX style file v2.2)
Triplets of supermassive black holes:Astrophysics, Gravitational Waves and Detection
Pau Amaro-Seoane ⋆ , Alberto Sesana , Loren Hoffman ,Matthew Benacquista , Christoph Eichhorn , , Junichiro Makino & Rainer Spurzem , , Max-Planck Institut f¨ur Gravitationsphysik (Albert-Einstein-Institut), Am M¨uhlenberg 1, D-14476 Potsdam, Germany andInstitut de Ci`encies de l’Espai (CSIC-IEEC), Campus UAB, Torre C-5, parells, na planta, ES-08193 Bellaterra, Barcelona, Spain Penn State University, 104 Davey Lab, Northwestern University, Dearborn Observatory, 2131 Tech Drive, Evanston, IL 60208-2900, USA Center for Gravitational Wave Astronomy, University of Texas at Brownsville, Brownsville, TX 78520, USA Institut f¨ur Raumfahrtsysteme, Universit¨at Stuttgart, Pfaffenwaldring 31, D-70550 Stuttgart, Germany Division of Theoretical Astronomy, National Astronomical Observatory, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan National Astronomical Observatories of China, Chinese Academy of Sciences, 20A Datun Lu, Chaoyang District, 100012, Beijing, China Kavli Institute for Astronomy and Astrophysics, Peking University, China Astronomisches Rechen-Institut, M¨onchhofstraße 12-14, 69120, Zentrum f¨ur Astronomie, Universit¨at Heidelberg, Germany
31 October 2018, submitted to MNRAS
ABSTRACT
Supermassive black holes (SMBHs) found in the centers of many galaxies are understoodto play a fundamental, active role in the cosmological structure formation process. Inhierarchical formation scenarios, SMBHs are expected to form binaries following themerger of their host galaxies. If these binaries do not coalesce before the merger with athird galaxy, the formation of a black hole triple system is possible. Numerical simula-tions of the dynamics of triples within galaxy cores exhibit phases of very high eccentric-ity (as high as e ∼ . LISA ) isestimated using several population models of SMBHs with masses & M ⊙ . Assuming10% or more of binaries are in triple systems, we find that up to a few dozen of thesebursts will produce residuals > − & M ⊙ ) with LISA isnegligible.
Key words: black hole dynamics gravitational waves cosmology: theory pulsars:general
It is well established that most galaxies host supermas-sive black holes (SMBHs) in their centers (Richstone et al.1998). In the past decade, compelling evidence of thecorrelation between the mass of the central SMBH and ⋆ e-mail: [email protected] (PAS) the bulge velocity dispersion and luminosity has been col-lected (Ferrarese & Merritt 2000; Gebhardt, et al. c (cid:13)
31 October 2018, submitted to MNRAS RAS
Pau Amaro-Seoane et al. tral SMBH evolution in which every merger of galaxies leadsquickly to coalescence of their central black holes can quan-titatively reproduce both the SMBH mass-bulge luminosityrelation (Kauffmann & Haehnelt 2000) and the SMBH mass-velocity dispersion relation (Haehnelt & Kauffmann 2000).In this general picture, if both of the galaxies involvedin a merger host a SMBH, then the formation of a SMBHbinary is an inevitable stage of the merging process. Fol-lowing the merger, the two black holes sink to the cen-ter of the merger remnant because of dynamical friction(Begelman, Blandford & Rees 1980). When the mass (eitherin gas or stars) enclosed in their orbit is of the order of theirown mass, they start to feel the gravitational pull of eachother, forming a bound binary. The subsequent binary evolu-tion is, however, still unclear. In order to coalesce, the binarymust shed its binding energy and angular momentum; a dy-namical process known in literature as ‘hardening’. A crucialpoint in assessing the fate of the binary is the efficiency withwhich it transfers energy and angular momentum to the sur-rounding gas and stars.The case of SMBH binaries in stellar environments hasreceived a lot of attention in the last decade. The system isusually modeled as a massive binary embedded in a stellarbackground with a given phase space distribution. The re-gion of phase space containing stars that can interact withthe SMBH binary in one orbital period is known as the losscone (Frank & Rees 1976; Amaro-Seoane & Spurzem 2001;Milosavljevi´c & Merritt 2003). As the binary evolves, it ejectsstars on intersecting orbits via the so called ‘slingshot mecha-nism’, causing a progressive emptying of the loss cone, whichultimately increases the hardening time scale. Without an ef-ficient physical mechanism for repopulating the loss cone, thebinary will never proceed to small separations where coales-cence induced by gravitational radiation takes place withina Hubble time. This is known as the stalling or ‘last parsec’problem (Milosavljevi´c & Merritt 2001).In the last decade, several solutions to the stallingissue have been proposed. Axisymmetric or triaxial stel-lar distributions may significantly shorten the coalescencetimescale (Yu 2002; Merritt & Poon 2004; Berczik, et al. et al. (2006). These orbits produce centrophilic stel-lar orbits and, therefore, replenish the loss-cone. However,more recent calculations by Amaro-Seoane & Santamaria(2009) of the outcome of the merger of two clusters ini-tially in parabolic orbits (Amaro-Seoane & Freitag 2006)have not been able to reproduce the rotation necessary tocreate the unstable bar structure. Other studies have in-voked eccentricities of the binary to refill the loss cone,since this effect could alter the cross section for super-elastic scatterings (thus altering the state of the loss cone)and shorten the gap to the onset of gravitational radia-tion effects (e.g.: Hemsendorf, Sigurdsson, & Spurzem 2002;Aarseth 2003a; Berczik, et al. > M ⊙ ),the inferred coalescence timescales in a stellar dominated en-vironment are of the order of few Gyrs, indicating that SMBHbinaries may be relatively long living systems. If the typicaltimescale between two subsequent mergers is comparable theSMBH binary lifetime, then a third black hole may reachthe nucleus when the binary is still in place, and the forma-tion of SMBH triplets might be a common step in the galaxyformation process. Recent studies of galaxy pairs lead to theconclusion that 30 −
70% of present day massive galaxies haveundergone a major merger since redshift one (Bell et al. 2006;Lin et al. 2008), where ’major’ means with baryonic mass ra-tio of the two components larger than 1 / / Nbody simulations,we study the dynamical evolution of the system, findingsurprisingly high eccentricities of the inner SMBH binary(up to e > . et al. ∼ yrs, Hoffman & Loeb (2007)), and final coales-cence is more common than ejection, confirming analyticalresults by Makino & Ebisuzaki (1996). We model at the lead-ing quadrupole order (Peters & Mathews 1963) the bursts ofgravitational radiation emitted in the highly eccentric phase,assessing detectability with future GW experiments. Adopt-ing cosmologically and astrophysically motivated models forSMBH formation and evolution, we estimate reliable eventrates.In order to cover the low frequencies generated bythe expected cosmological population of coalescing SMBHbinaries (e.g., Wyithe & Loeb 2003; Sesana et al. 2004,2005; Sesana, Volonteri & Haardt 2007) or plunges of com-pact objects such as stellar black holes on to supermas-sive ones (see e.g. Amaro-Seoane et al. 2007, for a reviewand references therein), the space-born observatory LISA (Bender et al. 1998) has been planned to be covering therange of frequencies of ∼ − − − Hz. Moving toeven lower frequencies, the Parkes Pulsar Timing Array(PPTA, Manchester 2006, 2008), the European Pulsar Tim-ing Array (EPTA, Janssen et al. 2008) and the NorthAmerican Nanohertz Observatory for Gravitational Waves(NANOGrav, Jenet et al. 2009) are already collecting dataand improving their sensitivity in the frequency range of c (cid:13)
31 October 2018, submitted to MNRAS RAS, MNRAS , 1–15 riplets of SMBHs: Astrophysics, GWs and detection ∼ − − − Hz, and in the next decade the planned SquareKilometer Array (SKA, Lazio 2009) will provide a majorleap in sensitivity.Throughout this paper we consider only very massivesystems, with total mass ∼ M ⊙ . Our goal is to in-vestigate if the high frequency nature of eccentric burstscan provide information about systems which would oth-erwise emit outside the frequency windows of the plannedGW experiments quoted above, by shifting wide (sepa-ration & . > M ⊙ ) sys-tems into the LISA domain. We note that the bursts an-alyzed here are different from the ‘bursts with memory’,which arise during the actual coalescence of SMBH binariesand are discussed in Pshirkov, Baskaran & Postnov (2009)and van Haasteren & Levin (2009).The structure of the paper is as follows. In Section 2, wedescribe our comprehensive study of the dynamics of triplesystems and investigate the eccentricity evolution of the innerbinary by using direct-summation N − body techniques and astatistical 3-body sample calibrated on the N − body results.In Section 3, we model the GW signal produced by eccentricbursts and we introduce observable quantities for PTAs and LISA . In Section 4 we construct detailed populations of emit-ting SMBH binaries and triplets, and we discuss our resultsin terms of signal observability and detection rates in Section5. Lastly, we briefly summarize our results in Section 6.
In modeling the dynamics of black hole triple systems withinthe centres of galaxy merger remnants, direct N -body inte-grations provide the most accuracy but are the most com-putationally expensive. We performed eight direct N -bodycalculations and used these to test the validity of an ap-proximation scheme involving three-body SMBH dynamicsembedded in a smoothed galactic potential with dynamicalfriction and gravitational radiation modeled by drag forces. N − body calculations The direct-summation
Nbody method we employed for allthe calculations includes the
KS regularisation . Thus, whentwo particles are tightly bound to each other or the separa-tion between them becomes very small during a hyperbolicencounter, the system becomes a candidate to be regularisedin order to avoid problematical small individual time steps(Kustaanheimo & Stiefel 1965). This procedure was later ex-ported to systems involving more than two particles. Inparticular, the
KS regularisation has been adapted to iso-lated and perturbed 3– and 4–body systems—the so-called triple (unperturbed 3-body subsystems), quad (unperturbed4-body subsystems) and the chain regularisation . The latteris invoked in our simulations whenever a regularised pair hasa close encounter with another single star or another pair(Aarseth 2003b).The basis of direct
Nbody codes relies on an improvedHermite integrator scheme (Aarseth 1999) for which we neednot only the accelerations but also their time derivative. Thecomputational effort translates into accuracy so that we canreliably keep track of the orbital evolution of every particle
Figure 1.
Initial conditions for the different direct N -body sim-ulations in the R P , V plane. For each simulation we choose theseparation between two SMBHs to be substantially smaller thanthe distance to the third SMBH. The initial parameters are se-lected in such a way that they are random but with an initialvelocity smaller than the escape velocity V esc and with a radius r < R P . We also show the circular velocity in the figure, V circ .Initially The set of 3 SMBHs of simulation A is represented withred solid bullets; for simulation B with green solid triangles; forsimulation C with cyan open triangles; for simulation D with blueopen squares; for simulation E with pink open circles; for simula-tion H with solid orange squares and for simulations F and G withblack crosses. In the case of simulations A, B, C, D, E, F and Gthe three SMBHs are set initially in a planar configuration. In thecase of simulation H they have a z component different from zeroin both the coordinates and velocities. In the cases of simulationsA, B, C, E and H we slightly modified the positions of the symbolsin order to avoid an overlap in our system. In order to make a highly accurate estimateof the eccentricity evolution of the SMBH system, we do notemploy a softening to the gravitational force (i.e. substitutingthe 1 /r factor with 1 / ( r + ǫ ), where r is the separationand ǫ the softening parameter) that weakens the interactionat small separations.The initial conditions for the set of three SMBHs, usedto conduct an exploration of the initial parameter space areshown in Table 2.1. For the stellar system, we use a Plum-mer model (Plummer 1911), which is an n = 5 polytropewith a compact core and an extended outer envelope. In thismodel the density is approximately constant in the centre anddrops to zero in the outskirts, φ = − GM ⋆ / p r + R , with M ⋆ the total stellar mass. This defines the Plummer radius R P . We depict the initial conditions in Figure 1, relative tothe circular and escape velocity of the Plummer potential. Wepresent results from eight direct numerical simulations, oneusing 512,000 stars using the special-purpose GRAPE6 sys-tem and the remaining simulations using Beowulf PC clustersand the AEI mini-PCI GRAPE cluster Tuffstein . c (cid:13)
31 October 2018, submitted to MNRAS RAS, MNRAS , 1–15
Pau Amaro-Seoane et al.
Model A B C D E F G H N ⋆ R/R P V/V esc
V/V circ
Table 1.
Initial conditions for the set of three SMBHs in each of the eight direct N -body simulations. N ⋆ is the total number of starsemployed in the simulation, V their velocity, R their position, R P the position of the SMBHs in terms of the Plummer radius, V esc theirescape velocity and V circ is their circular velocity. The mass of the SMBHs in N − body units is 1 · − , they are equal-mass and the massof a star 1 . · − While direct N -body simulations yield a very accurate result,they should be seen as a way to calibrate and test faster,more approximate simulations which can exhaustively coverthe parameter space and provide good statistics. We notethat the SMBHs in the N − body simulations are equal-massand (with the exception of simulation H), all of the systemsstudied with this method are coplanar. This was done becausesetting all SMBHs on a single plane accelerates the dynamics,shortening the integration time.In general, one wants to explore the whole parameterspace, including non coplanar systems with different SMBHmasses. For this purpose, we performed an ensemble of 1000three-body experiments, with the three Euler angles of theouter orbit sampled uniformly and a distribution of massratios motivated by Extended Press-Schechter theory (withtypical mass ratios m : m : m around 3.5:1). In each exper-iment we computed the Newtonian orbits of three SMBHsembedded in a smooth galactic potential and added dragforces to account for gravitational radiation and dynamicalfriction. We also included coalescence conditions when ei-ther of the two SMBHs pass within three Schwarzschild radiiof each other, or the gravitational radiation timescale be-comes short relative to the orbital period of the binary. Closetriple encounters were treated using a KS-regularised few-body code provided by Sverre Aarseth (Mikkola & Aarseth1990, 1993), while the two-body motion in between close en-counters was followed with a simple 4th-order Runge-Kuttaintegrator. See Hoffman & Loeb (2007) for further details onthe code. The initial conditions are those for the canonical ICs as in Hoffman & Loeb (2007). We performed each runtwice—once with gravitational radiation drag and the coa-lescence conditions, and once without.The 3-body experiments are divided into two computa-tional regimes based upon a dimensionless parameter, α , thatmeasures the relative tidal perturbation to the inner binaryby the interloper at apoapsis: α = 2 R M single M bin , smaller D , (1)where R apo is the apoapsis separation of the two inner binarymembers, M single is mass of the single SMBH, M bin , smaller is mass of the smaller inner binary member and D single is thedistance of the interloper (single SMBH) from the inner bi-nary centre-of-mass. In the limit α →
0, we know that theperiod of the inner binary is perfectly Keplerian (plus gravita-tional radiation), since the perturbation to the force from thethird body is negligible, and thus we can do orbit-averagedintegration instead of precisely integrating the trajectories ofall three bodies. The two regimes are defined as follows:(i)
The first regime corresponds to when a 3-body inter-action is taking place (defined by α > − ), an extremelyconservative criterion for when we need to do the full three-body integration.(ii) The second regime corresponds to when the singleSMBH and binary are wandering separately through thegalaxy ( α < − ), often on the order of a Hubble time.In regime (i) the 3-SMBH dynamics is integrated using SverreAarseth’s high-precision, regularised CHAIN code. Gravita-tional radiation and stellar-dynamical friction are treated asperturbing, velocity-dependent forces on the three separatebodies. In regime (ii) the separate orbits of the single andbinary centre-of-mass are followed using a simple 4th-orderRunge-Kutta integrator and the evolution of the binary semi-major axis and eccentricity are evolved using orbit-averagedequations.In computational regime (ii) (between 3-body encoun-ters), the stellar interactions are treated using the hard bi-nary prescription of Quinlan (1996). The eccentricity evo-lution of the inner binary under stellar interactions for nearequal-mass hard binaries is quite weak so it is neglected en-tirely and only the binary semi-major axis is evolved understellar interactions. The eccentricity is evolved under gravita-tional radiation as given by Peters (1964). Dynamical frictiontends to increase the eccentricity of a binary in the supersonicregime (where the orbital speed exceeds the stellar veloc-ity dispersion), and to circularize it in the subsonic regime.Since the triple SMBH system starts out supersonic in these3-body experiments, this effect produces a slight increase inthe outer binary eccentricity during the initial inspiral of thethird SMBH. Although we neglect the eccentricity evolutionof the inner binary during this phase, we find that its eccen-tricity is thermalized by the first resonant 3-body encounter. c (cid:13)
31 October 2018, submitted to MNRAS RAS, MNRAS , 1–15 riplets of SMBHs: Astrophysics, GWs and detection −3 −2 −1 −4 −3 −2 −1 C u m u l a t i ve f r ac t i on o f t i m e (a) all, GRi < 39 ° , GRall, Newti < 39 ° , Newtthermal 10 −3 −2 −1 f t h e r m a l / f ac t u a l (b) GRNewt10 −4 −3 −2 −1 −5 −4 −3 −2 −1 a (1 − e) (pc) C u m u l a t i ve f r ac t i on o f t i m e (c) GRNewt 10 −5 −4 −3 −2 −1 a (1 − e) (Schwarzschild radii) C u m u l a t i ve f r ac t i on o f t i m e (d) GRNewt Figure 2.
Cumulative fraction of time for the set of 1000 three-body simulations. The red, solid line corresponds to all simulations withthe 2.5 post-Newtonian correction term; the blue, solid line corresponds to the same simulations but without it (i.e. purely Newtonian);the dot-dashed red curve is like the first case but taking into account only systems in which the third SMBH had an inclination belowthe critical Kozai angle of 39 ◦ , to compare with the direct N -body simulations of Fig.(4); the dot-dashed blue curve is the same, but forthe Newtonian cases and the dot-dashed black curve corresponds to the thermal distribution, since the direct N -body simulations do nothave relativistic correction terms. The top left panel shows the eccentricity computed from the instantaneous positions and velocities ofthe two binary members, and the top right panel shows the ratio of the actual to the thermal distribution. For the Newtonian runs, thedistribution is within a factor of 2 of thermal down to 1 − e = 0 .
01 and within a factor of 3 of thermal down to ∼ − e = 0 . The timescale of the chaotic 3-body interactions ( ∼ yr) ismuch shorter than the stellar-dynamical timescale ( 10 − yr), so any effect of stellar interactions during these encoun-ters is completely negligible. Thus, stellar interactions playonly two roles in our 3-body simulations:a They bring the third BH in to interact with the innerbinary from an initial marginally stable hierarchical tripleconfiguration;b During phase (ii), stellar-dynamical interactions gradu-ally decrease the binary semi-major axis, so that it enters thenext three-body encounter harder than it left the last one ifthe time between encounters is long ( > yrs). The binary may even coalesce between encounters due to this gradualshrinking.Consequently, the distribution in eccentricities is bestestimated using the fraction of time that binaries spend ata given eccentricity while in computational regime (i) sincethere is minimal evolution of the eccentricity while in regime(ii). We note that the overwhelming majority of the timespent in regime (i) is still spent with α small enough that thesystem can be though of as a separate inner and outer binary,and so the instantaneous inner binary semi-major axis andeccentricity are well-defined. c (cid:13)
31 October 2018, submitted to MNRAS RAS, MNRAS , 1–15
Pau Amaro-Seoane et al. log ( 1−e ) l og [ t g r ( y r s ) ] −3.5 −3 −2.5 −2 −1.5 −1 −0.5 002468101214 −8−6−4−2log ( 1−e ) l og ( t g r / t o r b ) log P −3.5 −3 −2.5 −2 −1.5 −1 −0.5 0−2024681012 −8−6−4−2 Figure 3.
Time spent at different locations in the t gr - (1-e) plane in the Newtonian simulations (blue lines) of Figure 2. P = t/t close is the probability of finding the binary in a given bin at a randomly chosen time, where t close is the total simulation time spent in close3-body encounters. The simulation being presented here is the same as the Newtonian simulation (blue lines) in Figure 2. The purposeof this figure is to understand whether gravitational radiation is the main reason for the fall-off in the eccentricity distribution, since thehighest eccentricity systems have short coalescence times and quickly disappear. t gr is in years in the upper panel and in orbital periodsin the lower panel (the typical resonant encounter takes on order 10 yrs). Note that 1-e = 0.002 is about where t gr falls to less than anorbital period Figure 2a shows the fraction of the time that the binary (clos-est SMBH pair) spends above a given eccentricity during closeencounters, averaged over 1000 three-body experiments. Thered solid curves are the results of the standard runs includinggravitational radiation drag and coalescence conditions whilethe blue solid curves are the “Newtonian” case with theseeffects neglected. The dashed red and blue curves are aver-aged over only those experiments where the initial BH con-figuration has inclination < ◦ , the critical angle for Kozaioscillations, for comparison with the direct N -body simula-tions that use coplanar initial conditions. The black (dot-dashed) line shows the thermal distribution of eccentricitiesfor reference purposes. Three-body interactions result in athermal distribution of eccentricities, truncated at very higheccentricities by coalescence in collisions when gravitationalradiation is included. The similarity of the dashed and solidcurves shows that once the initial secular evolution is over, the system quickly thermalizes and little memory of the ini-tial configuration is maintained. Figure 2b shows the ratioof the thermal to the actual distribution as a function of ec-centricity. The first close encounter in each experiment hasbeen excluded from these plots, since it begins from a stablehierarchical triple configuration and includes a long period ofsecular evolution, whereas chaotic three-body encounters arethe focus of this work.The runs with and without gravitational radiationclosely follow each other and the thermal distribution up toeccentricities e ∼ .
99. At higher eccentricities, the Newto-nian distribution remains within a factor of 2 to 3 of thethermal distribution, but the gravitational radiation curvefalls off sharply, since these high-eccentricity systems coalescequickly through emission of gravitational waves. To verify thisinterpretation of Figure 2, we plot the time spent at differentlocations in the t gr -(1 − e ) plane in Figure 3, where the grav-itational radiation timescale is computed from the instanta- c (cid:13)
31 October 2018, submitted to MNRAS RAS, MNRAS , 1–15 riplets of SMBHs: Astrophysics, GWs and detection neous ( a, e ) of the binary. The eccentricity where the red andblue curves diverge in Figure 2b (1 − e ∼ .
01) is the valuewhere the typical t gr falls to just a few orbital periods, so thatthe binary can coalesce quickly by gravitational radiation be-fore the third body scatters it on to a lower-eccentricity or-bit. Figures 2c-d show the time spent by the binary at variouspericenter separations, in pc and in Schwarzschild radii. Notethat the gravitational radiation curve diverges sharply fromthe Newtonian one when the pericenter separation reaches ∼
100 Schwarzschild radii. We show the fraction of time spentat different eccentricities and the fraction of time spent at dif-ferent pericenter separations for the N − body simulations inFigure 4. The qualitative features are retained. In this section we make use of the leading Newtonian orderderivation of the GW radiation from eccentric binaries, asdescribed in Peters & Mathews (1963). We also use of ge-ometric units, with G = c = 1. Consider a system withmasses M < M orbiting with an orbital rest frame fre-quency f r = ω/ π , and with eccentricity e ; at the quadrupoleleading order, the luminosity ˙ E emitted by the system aver-aged over one complete orbit is:˙ E = 325 M / (2 πf r ) / F ( e ) = ˙ E c F ( e ) (2)where F ( e ) = ∞ X n =0 g ( n, e ) = 1 + e + e (1 − e ) / , (3)and M = M / M / / ( M + M ) / is the chirp mass of thesystem. ˙ E c , defined by the right-hand side in equation (2), isthe luminosity emitted by a circular binary orbiting at thesame frequency f r . The binary radiates GWs in the wholespectrum of harmonics f r,n = nf r ( n = 1 , , ... ), and therelative power radiated in each single harmonic is describedby the function g ( n, e ), defined as: g ( n, e ) = n „ B n + ` − e ´ A n + 43 n J n ( ne ) « , (4)where J n ( x ) are the Bessel functions and A n and B n are alsodefined in terms of the J n as: B n = J n − ( ne ) − eJ n − ( ne ) + 2 n J n ( ne )+2 eJ n +1 ( ne ) − J n +2 ( ne ) (5) A n = J n − ( ne ) − J n ( ne ) + J n +2 ( ne ) . (6)The total luminosity of the source can then be written, usingequations (2) and (3), as the sum of the component radiatedat each single harmonic:˙ E = ∞ X n =0 ˙ E n = ∞ X n =0 ˙ E c g ( n, e ) . (7)Given a general GW characterised by the two polarisedcomponent waves h + and h × , the rms amplitude of the waveis defined as h = q h h + h × i , where h i denotes the averageover directions and over time. The flux radiated in the GW field is related to the derivatives of its amplitude componentsby the relation (Thorne 1987) dEdtdA = 116 π “ ˙ h + ˙ h × ” . (8)The sinusoidal nature of the waves implies h ˙ h + ˙ h × i =4 π f r h h + h × i . So that integrating equation (8) over a spher-ical surface of radius d L (the luminosity distance from thesource) centered at the source and averaging over an orbitalperiod, directly relates ˙ E to the wave rms amplitude. We canthen infer that the rms amplitude and the energy radiated inthe n-th harmonic are related as (Finn & Thorne 2000) h n = 1 πd p ˙ E n f r,n = 2 r M / nd L (2 πf r ) / p g ( n, e ) , (9)where d = d L / (1 + z ). In the limit of a circular orbit (i.e., g ( n, e ) = δ n, in the Kronecker- δ notation), equation (9) re-turns the usual sky-polarization averaged amplitude (Thorne1987). Since we are interested in an estimate of the detectability ofextremely eccentric binaries (induced by triple interactions)by means of pulsar timing (and possibly
LISA ) observations,we first introduce an extension of the characteristic amplitude to include eccentric binaries. Eccentric binaries emit pulsesof GWs at their periapsis passages, and the rms amplitudeof each harmonic is given by equation (9). However, h n is anaverage amplitude related to the average luminosity along theorbit. The actual relevant time for the burst is the periapsispassage timescale T p = (1 − e ) / T orb ( T orb is the binaryorbital period), and if the burst is detected, almost all theenergy radiated along the whole orbit is seen on the timescale T p . This means that the relevant detectable amplitude of eachharmonic during the burst is h obs , n = h n p T f r,n , (10)where the factor T = max( T orb , T obs ) takes into account thefact that, if T orb < T obs , multiple bursts are visible during theobservation. Equation (10) is a crude approximation, never-theless it catches the basic features of the observed signal:this is given by the rms amplitude of each single n-th har-monic multiplied by the square root of the cycles completedby the harmonic in an orbital period, assuming that the bi-nary orbit is a fixed ellipse and GW emission does not changethe orbital parameters.The search for GWs using pulsar timing data exploitsthe effect of gravitational radiation on the propagation ofthe radio waves from one (or more) pulsar(s). A pass-ing GW would imprint a characteristic signature on thetime of arrival of radio pulses (e.g. Sazhin 1978; Detweiler1979; Bertotti, Carr & Rees 1983), producing a so called timing residual . We refer the reader to Jenet et al. (2004)and Sesana, Vecchio & Volonteri (2009), (hereinafter SVV09)for a detailed mathematical description of the GW inducedresiduals. The residuals are defined as integrals of the GWduring the observation time. For a collection of harmonics,the residuals are given by: R ( T ) = Z T ∞ X n =0 „ α − β γ ) h + ,n + αβ (1 + γ ) h × ,n « dt, (11) c (cid:13)
31 October 2018, submitted to MNRAS RAS, MNRAS , 1–15
Pau Amaro-Seoane et al.
Figure 4.
Left panel: Fraction of time for the direct N -body simulations at which the binary black hole has eccentricity in the range1 − e for the 512,000 stars simulation (thin green curve) and for the lower-resolution simulations. The thick red curve corresponds to theaverage of these lower-resolution computations and the dashed curve to the thermal distribution. Right panel: Cumulative fraction of timespent on a certain periapsis distance for all direct N -body simulations following the same colour labelling. One N − body unit of distanceis U| R = 1 . where α , β , and γ are the direction cosines of the pulsarrelative to a Cartesian coordinate system defined with the z -axis along the direction of propagation of the gravitationalwave and the x and y axes defining the + polarization. Theharmonics of the two polarizations, h + ,n and h × ,n , can befound in Section 3.2 of Pierro, et al. (2001). The rms residual δt gw is then formally defined as p h R ( T ) i .A simple derivation of the average timing residual δt gw generated by a circular binary is given by SVV09. With thenotations adopted above, their equation (20) reads: δt gw ( f ) = r h πf r p fT obs , (12)where the observed frequency f is related to f r as f = f r / (1+ z ) (being z the redshift of the source), the factor √ fT obs takesinto account for the signal ’build-up’ with the square rootof the number of cycles, and p /
15 comes from the angleaverage of the amplitude of the signal (cf. equation (17)-(21)of SVV09). We can generalise this derivation to the case ofbursts produced by eccentric binaries, relating the h obs , n ofeach harmonic to the induced residual residual at its peculiarfrequency via: δt gw ( f n ) = r h obs , n πf r,n , (13)The total residual can then be assumed to be of the order: δt gw = ∞ X n =0 δt ( f n ) ! / . (14)The estimation given in equations (13) and (14) is justifiedbecause the integral in equation (11) gives products of sines and cosines of different harmonics, that drop to zero when av-eraged over the observation, leaving only a sum of the squaresignals produced by each single harmonic (those terms includ-ing cos (2 πf n t ) and sin (2 πf n t )). We shall plot, in Section 4, R ( T ) for selected eccentric bursts, and we will see that δt gw as defined by equations (13) and (14) gives a good estimateof the amplitude of the induced residual.For inferring LISA detectability, given h obs , n , an estimateof the signal to noise ratio (SNR) in the LISA detector isstraightforwardly computed as:SNR = 4 ∞ X n =0 h , n fS f , (15)where S f is the one-side noise spectral density of the detector.We adopted the S f given in equation (48) of Barack & Cutler(2004), based on the LISA
Pre-Phase A Report. We extendedthe sensitivity down to 10 − Hz and we considered detec-tion with two independent TDI interferometers (which im-plies a gain of a factor of two in S f ). The SNR computed inthis way may seem a poor approximation. However, we havechecked the SNRs against those obtained following the proce-dure given in Section V-B of Barack & Cutler (2004), wherethe binary is consistently evolved with orbit averaged post-Newtonian equations, and found agreement at a 20-30% level,which is acceptable since we are interested in a preliminaryestimation of source detectability . The difference is mainly due to the fact that the orbital param-eters change during the strong GW emission burst, and this is nottaken into account in equations (10) and (15).c (cid:13)
31 October 2018, submitted to MNRAS RAS, MNRAS , 1–15 riplets of SMBHs: Astrophysics, GWs and detection The previous derivation can be use to achieve a heuristicunderstanding of what we may expect to actually detect. Letus consider two binaries ‘1’ and ‘2’ with the same masses, andsemimajor axes related as a = a (1 − e ) (suppose ‘2’ is incircular orbit and ‘1’ on a very eccentric orbit, i.e. 1 − e ≪ averaged over anorbital period. The eccentric binary ‘1’ has an orbital period T ∝ a / . But, it emits GWs in a short burst of durationof the order of its periapsis passage that is T p ∝ [ a (1 − e )] / . The mean luminosity of the eccentric binary duringthe periapsis burst is then˙ E ,p = ˙ E T T p ∝ F ( e ) a (1 − e ) − / , (16)where we used the Newtonian relation to switch from f to a in equation (2), and we ignored the source redshift. Accordingto equations (8) and (9), we can write h ≈ q ˙ E ,p /f p , wherewe make the assumption that f p is the ‘dominant frequencyof the burst’, which corresponds to the ‘periapsis frequency’, f p ∼ f (1 − e ) − / . A circular binary with semimajor a simply emits a periodic wave with amplitude h ≈ p ˙ E /f ,where ˙ E ∝ /a . Remembering that a = a (1 − e ) and,consequently, f p ≈ f , the h /h ratio reads: h h ∼ " e + e (1 + e ) / / = O (1) . (17)Since PTAs detect a timing residual that is δt ∼ h/f , itfollows that δt /δt ∼ h /h . The timing residual caused bya burst that happens to be at the right frequency for PTA( ∼ − Hz), generated by a very eccentric binary with anorbital frequency f ≪ − , is then of the same order of theresidual caused by a circular binary emitting at f = 10 − .The signal is, however, quite different and it is spread over abroad frequency band. This heuristic consideration suggeststhat PTA detection of such extreme events may be ratherdifficult, because their signal may be overwhelmed by GWemitted by ’conventional’ binaries with shorter periods. Onthe other hand, we might expect some interesting effect for LISA , since this mechanism can boost the GW frequencyby more than three order of magnitudes and signals fromsystems that would emit at much lower frequencies, may beshifted into the
LISA domain.
To draw sensible predictions about the number of expecteddetectable GW bursts, we need to model the population oftriple systems that form during the SMBH hierarchical buildup. We start by considering the SMBH binary , population.We are mainly interested here in probing massive systems M = M + M > M ⊙ , so that we can use catalogs ofsystems extracted from the Millennium Run (Springel et al.2005). We employ the very same catalogs used in SVV09;the reader is referred to Section 2 of that paper for details, here we merely summarise the basics of the procedure. Wecompile catalogs of galaxy mergers from the semi-analyticalmodel of Bertone, De Lucia & Thomas (2007) applied to theMillennium Run. We then associate a pair of merging SMBHsto each merging pair of spheroids (elliptical galaxies or bulgesof spirals) according to four different SMBH-host prescrip-tions (Section 2.2 of SVV09). Here we consider the threeTu models presented in SVV09, in which SMBHs correlatewith the spheroid masses according to the relation givenby Tundo et al. (2007), and differ from each other in theadopted accretion prescription: the Tu-SA model (accretiontriggered on to the more massive black hole before the final co-alescence), the Tu-DA model (accretion triggered before themerger on to both black holes) and the Tu-NA model (ac-cretion triggered after the coalescence). We also investigatethe dependence on the adopted SMBH binary population byconsidering the La-SA and Tr-SA models (see SVV09 for de-tails). The catalogs of coalescing binaries obtained in thisway are then properly weighted over the observable volumeshell at each redshift to obtain the differential distribution d N/d M dzdt r , i.e. the coalescence rate (the number of coa-lescences N per unit proper time dt r ) in the chirp mass andredshift interval [ M , M + d M ] and [ z, z + dz ], respectively. The GW signal can be divided into two contributions—onefrom the binaries, and one from the triplets. We will referto the latter as bursting sources , since we consider the GWbursts they emit at the periastron in their eccentric phase. Inthis study, we consider the binary population emitting in thePTA domain to be composed of circular systems dynamicallydriven by GW emission only. The GW signal is then givenby (Sesana, Vecchio & Colacino 2008): h c ( f ) = Z ∞ dz Z ∞ d M d Ndzd M d ln f r h ( f r ) , (18)where h is the sky-polarization average of each single source(Thorne 1987), and d N/dzd M d ln f r is the instantaneouspopulation of comoving systems emitting in a given logarith-mic frequency interval with chirp mass and redshift in therange [ M , M + d M ] and [ z, z + dz ], and is given by: d Ndzd M d ln f r = (1 − F t ) d Ndzd M dt r dt r d ln f r , (19)where dt r d ln f r = 564 π / M − / f − / r . (20)In equation (19), F t is the fraction of coalescing binaries thathave experienced a triple interaction. This can be estimatedsimply by knowing the likelihood of forming triple systemsbecause of two subsequent mergers. The galaxy merger ratedrops dramatically at low redshift, and the typical timescalebetween two subsequent major merger could be as long as ∼ yr. This means that massive galaxies may have ex-perienced, on average, just one major merger since z = 1(see, e.g. Bell et al. 2006). If we assume survival time of abinary is ∼ yrs, adopting the simplifying assumptionsof uncorrelated mergers with a Poissonian delay distributionwith a characteristic time of 10 yr, the probability of hav-ing two subsequent mergers in a 10 yr time interval is ∼ . c (cid:13)
31 October 2018, submitted to MNRAS RAS, MNRAS , 1–15 Pau Amaro-Seoane et al. log ( 1 − e ) l og [ a ( − e ) ] ( p c ) −3.5 −3 −2.5 −2 −1.5 −1 −0.5 0−4−3−2−10 l og P −9−8−7−6−5−4−3−2 Figure 5.
Two dimensional joint probability distribution P ( r p , e )for the inner binary in the [ r p , e ] plane (where r p = a (1 − e )).The distribution is obtained averaging over the 1000 3-body ex-periments described in Section 2.2 We will consider two different situations, choosing the frac-tion of SMBH binaries experiencing a triple interaction to be F t = 0 . F t = 0 . F t , we can write the coalescence rate of bi-naries that have experienced a triple interactions as F t × d N/dzd M dt r . From the 3-body scattering presented in Sec-tion 2.2, we derive the joint probability distribution for theinner binary of having a certain periastron and a certain ec-centricity, P ( r p , e ). This quantity is plotted in figure 5 for ourset of the 1000 3-body realizations. This probability distribu-tion refers to systems with mean total mass ∼ × M ⊙ .To extend it to a wider range of masses, we assume a tripletlifetime T = 10 yrs independently of the masses (which arein a narrow range peaked around 10 M ⊙ in our case) andwe rescale P ( r p , e ) so that, in the GW dominated regime,elements in the ( r p , e ) space having the same coalescencetimescale T gw ( r p , e ), have the same probability value P ( r p , e ).Since T gw ∝ a / [ M M ( M + M )], assuming an invariant bi-nary mass ratio distribution in the relevant mass range (whichis a good approximation given the narrow mass range we aredealing with), the y axis in figure 5 is rescaled for any giventotal mass of the binary M according to ( M/ × M ⊙ ) / .We then compute the distribution of eccentric binaries emit-ting an observable burst as: N ( M , z, r p , e ) = d Ndzd M dt r × F t × T × P ( r p , e ) ×× min[1 , ( T obs /T orb )] . (21)Where the factor min[1 , ( T obs /T orb )] takes into account thefact that if the binary period is longer than the observa-tion time, only a fraction T obs /T orb of the systems is actuallybursting during the observation. The relevant frequency band for pulsar timing observations isbetween 1 /T obs and the Nyquist frequency 1 / (2∆ t ) – where∆ t is the time between two adjacent observations–, corre-sponding to 3 × − Hz - 10 − Hz. The frequency resolution bin is ∆ f = 1 /T obs , and we assume T obs = 10 yr through-out the paper. Every realistic frequency–domain computationof the signal has to take into account the frequency resolu-tion bin ∆ f of the observation. The signal is therefore eval-uated for discrete frequency bins ∆ f j centered at discretevalues of the frequency f j , where f ( j +1) = f j + ∆ f . Whatwe actually collect in our code is the numerical distribution∆ N/ ∆ z ∆ M ∆ f r , where ∆ f r = (1 + z )∆ f . The integral inequation (18) is then replaced as a sum over redshift andchirp mass, and the value of the characteristic strain at eachdiscrete frequency f j is computed as h c ( f j ) = X z X M ∆ N ∆ z ∆ M ∆ f r,j f r h ( z, M , f r )∆ z ∆ M , (22)where ∆ f r,j = (1 + z )∆ f j is the j th frequency bin shifted ac-cording to the cosmological redshift of the sources. Equation(22) is simply read as the sum of the squares of the char-acteristic strains of all the sources emitting in the observed frequency bin ∆ f j . If we produce a family of α = 1 , ..., K sources by performing a Monte Carlo sampling of the nu-merical distribution ∆ N/ ∆ z ∆ M ∆ f r of the emitting binarypopulation, the characteristic strain is computed as h c ( f j ) = K X α =1 h c,α ( z, M , f α,r ) Θ[ f α,r , ∆ f j (1 + z )] (23)where Θ[ f α,r , ∆ f j (1 + z )] = 1 if f α,r ∈ ∆ f j (1 + z ) andis null elsewhere. To recover equation (22), the character-istic amplitude of the individual source is given by: h c,α = h α f α,r / ∆ f r,j ≈ h α f j / ∆ f j = h α f j T obs ; ı.e., the sky andpolarization averaged amplitude square, multiplied by thenumber of cycles completed in the observation time. The in-duced rms residual of each individual source is then given byequation (12). Note that in the limit of large K (formally, K → ∞ ), h c ( f j ) computed according to equation (23) isindependent of T obs (because the increment of the contribu-tion of each single source according to the number of cyclescompleted during T obs is balanced by the fact that we sumover a frequency bin that is proportional to 1 /T obs ), and itsvalue coincides with the one obtained from the standard en-ergy based definition of h c ( f ) (Sesana, Vecchio & Colacino2008). On the other hand, when K is small (i.e. we sum overa small number of sources), fluctuations become importantin the computation of the signal in each frequency bin. Nu-merical computation according to equation (23) allows us toaccount for signal fluctuations, which are missing in the an-alytical definition of the characteristic amplitude of the GWspectrum (e.g. Phinney 2001), but are important in the ac-tual computation of the observed signal. Given h c , the in-duced rms timing residual produced by the whole emittingpopulation is simply given by h c ( f i ) / (2 πf i ).We generate a population of emitting binaries accord-ing to the numerical distribution ∆ N/ ∆ z ∆ M ∆ f r , and wesum all the h c,α contributions in every frequency bin to ob-tain the characteristic strain of the signal. We then generate,again using a Monte–Carlo sampling, a population of emit-ting eccentric binaries in triple systems from the distributiongiven in equation (21) and we compute their GW bursts andthe induced rms residuals according to equations (9, 10, 13,14). For the few systems reaching 10 − Hz with their higherharmonics, we also compute the SNR produced in the
LISA detector using equation (15), adopting the S f given in equa-tion (48) of Barack & Cutler (2004), extended downward to c (cid:13)
31 October 2018, submitted to MNRAS RAS, MNRAS , 1–15 riplets of SMBHs: Astrophysics, GWs and detection Figure 6.
Representation of all the relevant features of a Monte–Carlo generated signal. The Tu-DA model with F = 0 . − Hz as described in Section 3.1. We consider five differ-ent SMBH binary populations presented in SSV09 (Tu-SA,Tu-DA, Tu-NA, La-SA, Tr-SA) with two different fractionsof triplets F t = 0 . , .
5, for a grand total of 10 different mod-els. We run 50 (when F t = 0 .
5; 100 if F t = 0 .
1) independentMonte–Carlo realizations of each single model, which allowsus to perform a statistical study of the properties of the burst-ing sources. We consider only systems with e > .
66, becausehighly eccentric systems are those expected to burst at highfrequencies, where the contribution of the overall circular bi-nary population declines. And also because high eccentricitiesresult in a well defined burst shape which may be essentialto distinguish it from periodic sources.
All the relevant features of the signal are plotted in figure6 for a realization of the Tu-DA model with F t = 0 .
5. AMonte–Carlo generated signal is depicted as a blue jaggedline. The magenta points represent all the binary systems producing a δt gw > . ∼ resolvable and we mark it with a superposed red triangle.The red jagged line is the resulting stochastic level of thesignal, after the contribution from the resolvable sources hasbeen subtracted. The arc-like black (red) tracks represent themore luminous (resolvable) bursting systems in the realiza-tion. In this particular case there were five resolvable burstswith rms residual δt gw = 3 . , . , . , . , .
002 ns. How-ever, considering realistic PTA sensitivities achievable in thenear future ( ∼ resolvablesource in the frequency domain, assuming that a source isresolvable if its strain is larger than the sum of the strainsof all the other sources in that frequency bin. This definitionis, however, only appropriate for monochromatic sources. Avery eccentric burst, emitting a whole spectrum of harmonics,may not be the brightest source in any of the frequency bins,however, it may produce a significantly larger rms residualswith respect to other individual circular binaries. Moreover,in the time domain, the signature of these bursts is quitedifferent with respect to periodic circular binaries, resultingin long bumps or narrow well localized bursts (see figure 10).Given these caveats, we will also present results in terms oftotal number of sources, independent of their resolvabilityaccording to our definition.A sample of different realizations of the signal is col-lected in figure 7, for individual realizations of the three dif-ferent Tu models. Only the brightest sources are plotted inthis case. Given the small number of systems involved, theirphenomenology is quite variable. For example, the realiza-tion illustrated in the left-middle panel, shows three resolv-able bursts with δt gw & To quantify the statistics of the bursting sources, we cast theresults in terms of the cumulative number of sources as afunction of the timing residuals: N ( δt gw ) = Z ∞ δt gw dNd ( δt ′ gw ) d ( δt ′ gw ) , (24)where the distribution dN/d ( δt ′ gw ) is the average over the50 (100) Monte–Carlo realizations of each model. We com-pute this average both considering all the sources emittingover a given δt gw threshold (obtaining the total distributionof bursting sources), as well as considering only resolvablesources as defined in the previous section (obtaining the dis-tribution of bursting resolvable sources). In figure 8, N ( δt gw )for all the sources is shown. Depending on the adopted model,and on the fraction of triplets assumed, there are few hun-dred to few thousand binaries contributing to the signal ata level & c (cid:13)
31 October 2018, submitted to MNRAS RAS, MNRAS , 1–15 Pau Amaro-Seoane et al.
Figure 7.
Sample of individual Monte–Carlo realizations of the signal generated using Tu-DA (left panels), Tu-SA (central panels) andTu-NA (right panels) models ( F = 0 . between 20 and 60 assuming F = 0 .
5, and, not surprisingly,a factor of 5 lower if we assume F = 0 .
1. If triple interactionsof SMBHs are common (say, F > . > . , resolvable sources only, the figures are not as promising. As shown in figure 9,a timing precision of 0.1 ns is needed to guarantee the de-tectability of a resolvable burst if F = 0 .
5. At a 1 ns level, we have less than one resolvable burst, we can then interpretthe results in terms of the probability of having such burstsin our observable Universe. This probability ranges from 2%to 50% depending on the adopted model, and the eccentric-ity distribution of these resolvable events is biased towardshigh values, peaking around e = 0 .
9. La-SA and Tr-SA givesimilar results both qualitatively and quantitatively, there-fore we don’t plot them in the figures in order to keep themclear. Again, we stress the fact that our definition of resolvablesource is rather arbitrary, and does not take into account for c (cid:13)
31 October 2018, submitted to MNRAS RAS, MNRAS , 1–15 riplets of SMBHs: Astrophysics, GWs and detection Figure 8.
Cumulative number N ( δt gw ) of circular binaries (thinlines) and bursting triplets (thick lines) emitting over a given δt GW threshold as a function of δt GW . In each panels the differ-ent linestyles refer to the Tu-SA (solid), Tu-DA (long–dashed) andTu-NA (short–dashed) models. The fraction of triplets assumed islabeled in each panel. the peculiar shape of the burst, we then consider these figuresas lower limits to the actual detectability of these bursts. To give a feeling of how the actual signals would appear,we also computed residuals in the time domain for selectedsources. To this purpose, we evolved the system using equa-tions (27)-(31) of Barack & Cutler (2004) assuming non spin-ning SMBHs. We then computed all the components h + ,n and h × ,n (following Pierro, et al. R ( T ) integrating equation (11). The actualshape of the residuals is rather complex and depends on thegeometry of the system: the relative orientation of the sourceto the pulsars (encoded in the direction cosines α , β , and γ in equation (11)); the polarization angle of the source Ψ; theinclination i ; the initial phase of the orbit Φ ; and an angle φ p describing the orientation of the periastron in the orbitalplane (see, e.g., Barack & Cutler 2004 for a definition of allthese quantities).Examples of the phenomenology of bursting sources aregiven in figure 10 for a sample of eccentric systems found inone selected realisation of the model Tu-DA. In the left pan-els we show the three brightest resolvable sources, while inthe right panels we show three of the brightest bursts whichwould be unresolvable according to our definition, becausetheir power spectra would be overwhelmed by the signal pro-duced by the standard circular binaries found in the realiza-tion. Parameters of the binaries are given in table 2. Burstscan be generated by very eccentric-long period binaries (as inthe two lower panels), or by relatively short-period systems(e.g., central left panel), in which case multiple bursts are vis- Figure 9.
Same as figure 8, but considering only the resolvablesources in the computation of N ( δt gw ) according to equation (24).Linestyle as in figure 8. M [M ⊙ ] M [M ⊙ ] f r [Hz] z e δt gw [ns]1 . × . × . × − . × . × . × − . × . × . × − . × . × . × − . × . × . × − . × . × . × − Table 2.
Parameters of the sources plotted in figure 10. Rows inthe table (from the top to the bottom), correspond to panels offigure 10 considered counterclockwise, starting from the upper leftpanel. ible in the observation times. The width of the burst dependson the periastron passage timescale: systems with T p ≪ T obs produce narrow features in the data stream (e.g., lower leftpanel), while systems with T p ≈ T obs give a characteristicbump shaping all over the data span (e.g. upper right panel).Given the integral nature of the signal (equation (11)), itsshape is also heavily dependent on Φ (e.g. the cumulativeresidual can be positive or negative depending on the binaryorbital phase at the beginning of the detection), on φ p andon Ψ. The inclination of the source i and its aperture angleto the pulsar, determine the amplitude of the signal. LISA
We also collected catalogs of systems bursting in the
LISA window, to check for detectability. Unfortunately, prospectsfor detection with
LISA are not as promising as for PTAs.In a total of 750 realization of the ten different models, wefound ∼
50 sources bursting in the
LISA window produc-ing an SNR > .
1. Unfortunately, none of them produced c (cid:13)
31 October 2018, submitted to MNRAS RAS, MNRAS , 1–15 Pau Amaro-Seoane et al.
Figure 10.
Examples of timing residual, found in a particularrealization of the model Tu-DA, computed according to equation(11). Different linestyles correspond to different aperture angle θ between the pulsar and the source; θ = π/ π/ π/ π/ i = π/ φ p is random and Φ is chosen so that the burst occurs duringthe observation. The rms residual computed according to equation(14) are also shown. Parameters of the sources are listed in table2. a SNR >
8, necessary for a confident detection. We thenconclude, that even with a consistent population of SMBHtriplets forming during the cosmic history, burst from mas-sive (say,
M ∼ M ⊙ ) eccentric binaries are unlikely to beproduced at a significant rate for LISA . On the other hand,if formation of triple systems was common in the past, forsystem in the
LISA mass range ( ∼ − M ⊙ ), very pe-culiar signals from coalescing eccentric binaries may be com-mon in the data stream. However, this is beyond the scopeof the present paper, where we focused on massive binaries( M > M ⊙ ) only. We have addressed in this work three different points in theevolution of triplets of SMBHs in the Universe: The Astro-dynamics of the system, the potential GW signature and thedetectability.We have performed eight different direct-summation N − body simulations, one including more than half a millionof particles, to calibrate 1,000 3-body scattering experiments,which include post-Newtonian corrections, in order to havea statistical description of the system. Both numerical toolsagree that the inner binary of SMBHs will go through a phaseof extremely high eccentricity, which is the motivation for therest of the work.These three-body excitations of episodic high eccentric-ity configurations of the close SMBH binary produce inter-esting GW bursts that may be detectable with forthcoming experiments such as PTAs and LISA . The extreme eccentrici-ties of such bursts on one hand would leave a very distinctivesignature, but on the other require the development of ap-propriate analysis techniques.To compute likely event rates, we extracted cataloguesof merging galaxies from the Millennium Run, and we pop-ulated them with SMBHs following the known MBH-bulgesrelations. We then estimated the fractions of triplets and theireccentricity distribution and we computed the induced signalsin both PTAs and the
LISA detector.We found that, depending on the details of the SMBHpopulation model, if the fraction of triplets is > .
1, few to ahundred of GW bursts would be produced at a > − LISA is essentially nil. However, we stress the fact that wefocused on systems with M > M ⊙ ; our results then sim-ply imply that it is extremely unlikely that a system whichwould normally emit outside the LISA range will produce aburst in the
LISA window because of resonant three body in-teractions. On the other hand, if a consistent fraction of lightbinaries ( M < M ⊙ ) is involved in triple systems, we mayexpect several eccentricity-driven coalescences to be observedby LISA . This eventuality would call for the development ofextremely eccentric templates ( e > .
9) for merging SMBHs,and of adequate analysis techniques to extract the signal.
ACKNOWLEDGMENTS
The work of PAS has been supported by the Deutsches Zen-trum f¨ur Luft- und Raumfahrt. PAS is indebted with SterlPhinney for discussions which motivated the work and withMarc Freitag for comments on the article and help with someof the diagrams. PAS, AS and MB acknowledge the sup-port of the Aspen Center for Physics. Some of the numer-ical simulations were done with the
Tuffstein cluster lo-cated at the Max-Planck Institut f¨ur Gravitationsphysik (Al-bert Einstein-Institut). RS and PAS acknowledge comput-ing time on the GRACE cluster in Heidelberg (Grants I/80041-043 of the Volkswagen Foundation and 823.219-439/30and /36 of the Ministry of Science, Research and the Artsof Baden-W¨uurttemberg). MB acknowledges the support ofNASA grant NNX08AB74G and the Center for GravitationalWave Astronomy, supported by NSF award 0734800.
REFERENCES
Aarseth, S. J., PASP 111, 1333, 1999Aarseth, S. J. 2003a, AAS, 285, 367 c (cid:13)
31 October 2018, submitted to MNRAS RAS, MNRAS , 1–15 riplets of SMBHs: Astrophysics, GWs and detection Aarseth, S. J. 2003b, Gravitational N-Body Simulations,pp. 430. ISBN 0521432723. Cambridge, UK: CambridgeUniversity Press, November 2003Amaro-Seoane, P., & Spurzem, R. 2001, MNRAS, 327Amaro-Seoane, P.; Freitag, M., ApJ, 653, Issue 1, pp. L53-L56 (2006)Amaro-Seoane P., Gair J. R., Freitag M., Miller M. C., Man-del I., Cutler C. J., Babak S., 2007, Classical and QuantumGravity, 24, 113Amaro-Seoane, P.; Miller, M. C.; Freitag, M., ApJ, 692, Is-sue 1, pp. L50-L53 (2009)Amaro-Seoane, P. and Eichhorn, C. and Porter, E. andSpurzem, R., 2009, accepted for publication in MNRAS,2009Amaro-Seoane, Pau; Santamaria, Lucia, ArXiv e-printsarXiv:0910.0254Barack L. & Cutler C., 2004, PhRvD, 082005Begelman, M.C., Blandford, R.D., Rees, M.J. 1980, Nature287, 307Bell E. F., Phleps S., Somerville R. S., Wolf C., Borch A. &Meisenheimer K., 2006, ApJ, 652, 270Bender P. L. et al., 1998
LISA Pre-Phase A Report; SecondEdition , MPQ 233Berczik, P.; Merritt, D.; Spurzem, R.; Bischof, H.-P., 2006,ApJ, 642, L21Bertone S., De Lucia G. & Thomas P. A., 2007, MNRAS,379, 1143Bertotti B., Carr B. J. & Rees M. J., 1983, MNRAS, 203,945Detweiler, S.L. 1979, ApJ, 234, 1100–1104Dotti, M., Colpi, M., & Haardt, F., 2006, MNRAS, 367, 103Escala, A., Larson, R. B., Coppi, P. S., & Mardones, D.,2005, ApJ, 630, 152Ferrarese, L., Merritt, D. 2000, ApJ 539, L9Finn L. S. & Thorne K. S., 2000, PhRvD, 124021Frank, J. & Rees, M. J. 1976, MNRAS, 176, 633Gebhardt, K., Richstone, D., Kormendy, J., Lauer, T.R.,Ajhar, E.A., Bender, R., Dressler, A., Faber, S.M., Grill-mair, C., Magorrian, J., Tremaine, S. 2000b, AJ 119, 1157Haehnelt, M.G.; Kauffmann, G.: 2000, MNRAS, 318, L35.Hemsendorf, M., Sigurdsson, S., & Spurzem, R. 2002, ApJ,581, 1256Hoffman, L. & Loeb, A., 2007, MNRAS, 377, 957Janssen G.H. et al., 2008, in 40 YEARS OF PULSARS:Millisecond Pulsars, Magnetars and More. AIP, ConferenceProceedings, Volume 983, 633Jenet, F. A., Lommen, A., Larson, S. L., & Wen, L. 2004,ApJ, 606, 799Jenet, F. A. et al., 2009, preprint (arXiv:0909.1058)Kauffmann, G. & Haehnelt, M.G. 2000, MNRAS, 311, 576Kustaanheimo, P.E. & Stiefel, E. L. 1965, Publ. Astron.Obs. Helsinki 110, 204-219, J. Reine Angew. Math. 218,204Lazio, J., 2009, arXiv:0910.0632Lin, L. et al., 2008, ApJ, 681, 232Makino, J. & Ebisuzaki, T. 1996, ApJ, 436, 607Manchester, R.N., 2006, Chin. J. Astron. Astrophys., 6, 139Manchester R.N., 2008, in 40 YEARS OF PULSARS: Mil-lisecond Pulsars, Magnetars and More. AIP, ConferenceProceedings, Volume 983, 584Merritt, D. & Ferrarese, L. 2001, MNRAS, 320, L30Merritt, D. & Poon, M. Y., 2004, ApJ, 606, 788 Mikkola S. & Aarseth S., 1990, Celest. Mech. Dyn. Astron.,47, 375Mikkola S. & Aarseth S., 1993, Celest. Mech. Dyn. Astron.,57, 439Milosavljevi´c, M. & Merritt, D. 2001, ApJ, 563, 34Milosavljevi´c, M. & Merritt, D. 2003, ApJ, 596, 860Peters, P.C. and Mathews, J. 1963, Phys. Rev., 131, 434–440Peters, P.C., 1964, Phys. Rev., 136, 1224-1232Perets, H. B., Hopman, C., & Alexander, T., 2007, ApJ,656, 709Phinney E. S., 2001, arXiv:astro-ph/0108028Pierro, V., Pinto, I.M., Spallicci, A.D., Laserra, E., & Re-cano, F. 2001, MNRAS, 325, 358–372Plummer H. C., 1911, MNRAS, 71, 460Pshirkov M.S., Baskaran, D. & Postnov, K.A., 2009,arXiv:0909.0742Quinlan, G. 1996, NewA, 1, 35.Richstone, D. et al., 1998, Nature, 395, 14Sazhin M. V., 1978, Soviet Astron., 22, 36Sesana, A., Haardt, F., Madau, P., & Volonteri, M. 2004,ApJ, 611, 623.Sesana A., Haardt F., Madau P. & Volonteri M., 2005, ApJ,623, 23Sesana A., Volonteri M. & Haardt F., 2007, MNRAS, 377,1711Sesana A., Vecchio A. & Colacino C. N., 2008, MNRAS,390, 192Sesana A., Vecchio A. & Volonteri M., 2009, MNRAS, 394,2255Springel V. et al., 2005, Nature, 435, 629Thorne, K. S. 1987, in 300 Years of Gravitation, ed. S. Hawk-ing & W. Israel (Cambridge: Cambridge Univ. Press), 330Tremaine, S., et al. 2002, ApJ, 574, 740Tundo E., Bernardi M., Hyde J. B., Sheth R. K. & PizzellaA., 2007, ApJ, 663, 57Valtonen, M.J., Mikkola, S., Heinamaki, P., Valtonen, H.1994, ApJS, 95, 69van Haasteren, R. & Levin, Y., 2009, arXiv:0909.0954Wyithe J. S. B. & Loeb A., 2003, ApJ, 590, 691Yu, Q., 2002, MNRAS, 331, 935 c (cid:13)
31 October 2018, submitted to MNRAS RAS, MNRAS000