Gravitational Waves from the Inspiral of Supermassive Black Holes in Galactic-scale Simulations
Matias Mannerkoski, Peter H. Johansson, Pauli Pihajoki, Antti Rantala, Thorsten Naab
DDraft version December 11, 2019
Typeset using L A TEX twocolumn style in AASTeX63
Gravitational Waves from the Inspiral of Supermassive Black Holes in Galactic-scale Simulations
Matias Mannerkoski, Peter H. Johansson, Pauli Pihajoki, Antti Rantala,
1, 2 and Thorsten Naab Department of Physics, Gustaf H¨allstr¨omin katu 2, FI-00014, University of Helsinki, Finland Max-Planck-Institut f¨ur Astrophysik, Karl-Schwarzchild-Str. 1, D-85748, Garching, Germany
ABSTRACTWe study the orbital evolution and gravitational wave (GW) emission of supermassive black hole(SMBH) binaries formed in gas-free mergers of massive early-type galaxies using the hybrid tree–regularized N-body code
Ketju . The evolution of the SMBHs and the surrounding galaxies is followedself-consistently from the large-scale merger down to the final few orbits before the black holes coalesce.Post-Newtonian corrections are included up to PN3.5-level for the binary dynamics, and the GWcalculations include the corresponding corrections up to PN1.0-level. We analyze the significance ofthe stellar environment on the evolution of the binary and the emitted GW signal during the final GWemission dominated phase of the binary hardening and inspiral. Our simulations are compared to semi-analytic models that have often been used for making predictions for the stochastic GW backgroundemitted by SMBHs. We find that the commonly used semi-analytic parameter values produce largedifferences in merger timescales and eccentricity evolution, but result in only ∼
10% differences in theGW spectrum emitted by a single binary at frequencies f (cid:38) − yr − , which are accessible by currentpulsar timing arrays. These differences are in part caused by the strong effects of the SMBH binarieson the surrounding stellar population, which are not included in the semi-analytic models. INTRODUCTIONSupermassive black holes (SMBHs) with masses in therange M = 10 –10 M (cid:12) are believed to reside in thecenters of most, if not all, massive galaxies (for a review,see Kormendy & Ho 2013). There is also increasing ob-servational evidence for systems with multiple SMBHs,such as an SMBH binary system at z = 0 .
055 with aprojected separation of just ∼ z = 0 . ∼
140 pc (Deane et al.2014). The observed quasi-periodic outbursts of thequasar OJ287 have also plausibly been modeled as abinary system with a separation of only ∼ .
05 pc (Val-tonen et al. 2008; Dey et al. 2018).In the ΛCDM hierarchical picture of structure forma-tion, galaxies grow through mergers and gas accretion,resulting in situations with multiple black holes in thesame galaxy (e.g. Begelman et al. 1980; Volonteri et al.2003; Tremmel et al. 2018). Galaxy mergers are par-ticularly relevant for the most massive, slowly-rotatingearly-type galaxy population hosting the largest SMBHsin the Universe. These galaxies are believed to have as-
Corresponding author: Matias Mannerkoskimatias.mannerkoski@helsinki.fi sembled through a two-stage process in which the earlyassembly is dominated by rapid in situ star formationfueled by cold gas flows and hierarchical merging of mul-tiple star-bursting progenitors, whereas the later growthbelow redshifts of z (cid:46) a ∼
10 pc is formed in the center of the merging galaxypair. As the SMBH binary continues hardening starsbecome the primary scatterers, experiencing complexthree-body interactions that carry away energy and an-gular momentum from the SMBH binary system (e.g.Hills & Fullerton 1980). The largest uncertainty in thisprocess is the rate at which the ‘loss cone’ is filled, i.e.the region of parameter space where the stars have suffi-ciently low angular momenta to interact with the SMBHbinary. If the SMBH ‘loss cone’ is depleted, the binaryhardening is halted and we are faced by the so-calledfinal-parsec problem (Milosavljevi´c & Merritt 2001; Mer-ritt 2013). However, if the SMBH binaries are able to a r X i v : . [ a s t r o - ph . GA ] D ec reach smaller separations ( (cid:46) . f = 10 − –10 − Hz) and thusable to detect the final stages of the inspiral and coales-cence of SMBHs with masses of M • (cid:46) M (cid:12) (Amaro-Seoane et al. 2017).In the meantime before LISA becomes operational,pulsar timing arrays (PTAs) constitute the most promis-ing method for detecting gravitational waves emittedby SMBH binaries. PTAs attempt to detect GWs atnanohertz frequencies by measuring correlated offsets inthe arrival times of the highly regular pulses emittedby pulsars in the Milky Way (see e.g. Tiburzi 2018for an overview). It is expected that the large popu-lation of GW sources forms a stochastic gravitationalwave background (GWB), which has a nearly power-law spectrum (Phinney 2001) with a characteristic am-plitude of h c ∼ − at the reference frequency of f = 1 yr − ≈ × − Hz. The exact properties of theGW spectrum are affected by the eccentricity distribu-tion of the binaries and different environmental effects,such as stellar ‘loss-cone’ scattering and the viscous dragfrom circumbinary gas discs that will influence the bi-nary population (e.g. Sesana 2013; Kelley et al. 2017b;Burke-Spolaor et al. 2019). The PTA observations holdgreat promise for constraining the properties and forma-tion mechanisms of the SMBH binary population (e.g.Arzoumanian et al. 2018; Middleton et al. 2018; Chenet al. 2019).Previous work on the GW emission from mergingSMBH binaries can be classified into two categoriesbased on whether the final sub-parsec dynamics of thebinary is directly resolved or instead modeled usingsemi-analytic formulae. Studies of the GWB amplitudeare typically in the latter category, with the evolution due to gravitational wave emission implemented oftenusing the classic leading order formulae of Peters (1964).Effects of the environment may be included using sim-plified models which cannot accurately account for e.g.highly anisotropic stellar populations. The main differ-ences between the different semi-analytic studies lie inthe treatment of the population of binaries, using eitheranalytic formulae constrained by observations of galax-ies (e.g. McWilliams et al. 2014; Huerta et al. 2015; In-ayoshi et al. 2018) and observed candidate binary sys-tems (Sesana et al. 2018), or the results from cosmolog-ical simulations (Sesana et al. 2008; Salcido et al. 2016;Kelley et al. 2017b; Ryu et al. 2018).Studies where the dynamics of the SMBH binary isresolved down to the scales relevant for GW emissioneither model the environment using semi-analytic meth-ods (e.g. Bonetti et al. 2016, 2018), or they utilize N-body codes to directly simulate the interaction betweenthe binary and the surrounding field of stars, using ei-ther Newtonian gravity (Preto et al. 2011) or includingalso post-Newtonian (PN) corrections to the motion ofthe SMBHs to account for relativistic effects (Berentzenet al. 2009; Khan et al. 2018a). Typically, due to thesteep O ( N ) scaling of computational time with parti-cle number, these simulations are limited to a relativelylow number of particles, which typically restricts themto studying only the core regions of the galaxies sur-rounding the SMBHs. Some of these resolution issuescan be mitigated by using a combination of codes fordifferent phases of the galaxy merger and SMBH inspi-ral (Khan et al. 2016, 2018b).In this paper we utilize the recently developed hybridtree–N-body code Ketju (Rantala et al. 2017, 2018) todirectly compute the trajectories of the SMBHs and theresulting GW signal in gas-free mergers of massive early-type galaxies. The
Ketju code uses an algorithmicchain regularization (AR-CHAIN); (Mikkola & Merritt2006, 2008) method to efficiently and accurately com-pute the dynamics close to SMBHs, and combines it withthe fast and widely used tree code GADGET-3 (Springelet al. 2005; see also Karl et al. 2015 for a similar codecombination). The hybrid nature of
Ketju enables usto follow the SMBHs from before the galaxies mergeuntil the final few orbits before the coalescence of theSMBH binary that forms after the merger. In particu-lar, the distribution of stars is followed self-consistentlythroughout the simulation, correctly accounting for thechanging properties of the surrounding stellar popula-tion during both the dynamical friction dominated phaseand the stellar scattering driven phase.The main goal of the present work is to compare thegravitational wave emission from dynamically resolved
Ws from SMBHs in Galactic Simulations
Ketju simulations to semi-analyticmodels that have been used to compute the GW emis-sion for SMBH binaries in e.g. cosmological simulations,where the spatial resolution is insufficient to model theSMBHs directly (Kelley et al. 2017a,b). Our main focusfor the GW calculations is in the PTA frequency bandand the comparison with semi-analytic models allowsus to validate the semi-analytic models and directly ad-dress how important the additional accuracy gained byusing
Ketju is in the context of making predictions forPTA observations.The remainder of this paper is organized as follows:in section 2 we give a brief overview of the
Ketju codeand describe our merger simulations. In section 3 wepresent the methods used for analyzing the evolution ofthe SMBH binary, including the use of quasi-Keplerianorbital elements that account for post-Newtonian correc-tions as well as different simple binary evolution modelsused for comparison. Then, in section 4, we present thecombination of two methods used to compute the emit-ted gravitational waves in different phases of the binaryevolution. These methods are applied to the simulationsin section 5, and the implications of the obtained resultsare discussed in section 6. Finally, we present our con-clusions in section 7. NUMERICAL SIMULATIONS2.1.
The KETJU Code
The simulations analyzed in this paper have been runusing the recently developed
Ketju code (see Rantalaet al. 2017, 2018 for full code details), which is an ex-tension of the tree-SPH simulation code GADGET-3(Springel et al. 2005). The central idea of the codeis the inclusion of a regularized region around everySMBH particle, in which the non-softened gravitationaldynamics is computed using the regularized AR-CHAIN(Mikkola & Merritt 2008) integrator while the dynam-ics of the remaining particles is computed with theGADGET-3 leapfrog using the tree force calculationmethod.In practice the code operates by dividing simulationparticles into three categories. The SMBH and all stel-lar particles, which lie within a user defined chain ra-dius ( r chain ) are marked as chain subsystem particles.Particles that lie just outside the chain radius, but in-duce a strong tidal perturbation on the chain system aremarked as perturber particles. Finally, all the remain-ing particles that are far from any SMBHs are treated asordinary GADGET-3 particles with respect to the forcecalculation. The Ketju code also allows for both multi-ple simultaneous chain subsystems and several SMBHsin a given single subsystem. During each global GADGET-3 timestep the particlesin the chain subsystems are propagated using the AR-CHAIN algorithm (Mikkola & Merritt 2008; Rantalaet al. 2017). This algorithm has three main aspects:algorithmic regularization, the use of relative distancesto reduce round-off errors by organizing the particlesin a chain (e.g. Mikkola & Aarseth 1993) and finallythe use of a Gragg–Bulirsch–Stoer (GBS); (Gragg 1965;Bulirsch & Stoer 1966) extrapolation method that yieldshigh numerical accuracy in orbit integrations at a presetuser-given error tolerance level ( η GBS ). In essence, algo-rithmic regularization works by transforming the equa-tions of motion by introducing a fictitious time variablesuch that integration by the common leapfrog methodyields exact orbits for a Newtonian two-body problemincluding two-body collisions (e.g. Mikkola & Tanikawa1999; Preto & Tremaine 1999).2.2.
Post-Newtonian Corrections
The AR-CHAIN algorithm within the
Ketju codefeatures an extension of phase space with the help of anauxiliary velocity variable (Hellstr¨om & Mikkola 2010;Pihajoki 2015). This allows for the efficient implemen-tation of the velocity-dependent post-Newtonian correc-tions in the motion of the SMBH particles (e.g. Will2006) using an explicit integrator.Schematically the PN-corrected acceleration can bewritten as a = a N + a PN1 + a PN2 + a PN3 + a PN2.5 + a PN3.5 , (1)where the Newtonian acceleration a N is computed in-cluding the surrounding stellar particles, while the PN-terms only include contributions from other SMBHs.The PN-correction terms are labeled so that they areproportional to the corresponding power of the formalPN expansion parameter (cid:15) PN , i.e. | a i PN | ∝ (cid:15) i PN ∼ (cid:16) vc (cid:17) i ∼ (cid:18) R s R (cid:19) i , (2)where v and R are the relative velocity and separationof a pair of SMBHs, R s = 2 GM/c is the Schwarzschildradius corresponding to the total binary mass M , c isthe speed of light and G the gravitational constant. ThePN terms of integer order are conservative and are asso-ciated with conserved energy and angular momentum,while the half integer order terms are dissipative radi-ation reaction terms caused by the emission of gravita-tional radiation.In the Ketju simulations studied in this paper weuse PN-correction terms up to order PN3.5 derived fora binary system of arbitrary eccentricity in the modi-fied harmonic gauge as given in Mora & Will (2004).
433 Myr 100 pc 1104 Myr 100 pc 1286 Myr 100 pc5 pc 5 pc 0.5 pc
Figure 1.
A sequence of illustrative example snapshots of simulation run A (5:1 mass ratio merger). The main images showthe projected stellar density, while the insets show the stellar particles (blue), SMBHs (black), and sections of their trajectoriesin the regularized region (grey). The snapshots illustrate different characteristic phases of the SMBH binary evolution: theSMBHs first sink to the center of the merging galaxy due to dynamical friction (left panel) and form an eccentric bound binary(middle panel) that shrinks due to stellar scattering, finally entering the strongly relativistic regime (right panel) where thebinary shrinks due to GW emission, and other relativistic effects such as precession of the orbit are apparent as well. Thecorresponding timescales and spatial scales are indicated on the figure.
Terms depending on the spin of the black holes as wellas the lowest order cross term corrections for a systemconsisting of more than two particles have also been im-plemented in the
Ketju code, but are not used in thispresent work. The spin terms are ignored due to the factthat the spin of the SMBHs is largely determined by in-teractions with gas, which is not included in this study.Similarly, the cross terms as well as other PN-correctionsfor stellar particles are ignored due to the unphysicallystrong corrections resulting from stellar particles withvery large masses ( m (cid:63) (cid:29) M (cid:12) ). Thus post-Newtoniancorrections are only included in the interactions betweenSMBHs and only when they are within the same chainregion.We note that the missing higher order PN-correctionsmay potentially lead to significant effects over long peri-ods of time. Corrections up to the PN4-level have beenrecently derived (e.g. Bernard et al. 2018), but these arehighly impractical to implement due to the appearanceof the so called ‘tail’-terms that depend on the entirehistory of the binary in addition to its instantaneousstate. The lack of these higher order terms also leadsto the simulated motion of the SMBHs becoming unre-liable when (cid:15) PN ∼ . R = 6 R s , where the equations of motion arestill well behaved. However, as will be discussed in sec-tion 5, the actual limit of reliability for our purposes iscloser to R ∼ R s .2.3. Progenitor galaxies and merger orbits
Our merger progenitor galaxies consist of sphericallysymmetric, isotropic multi-component systems consist-ing of a stellar bulge, a dark matter halo and a centralSMBH. The stellar component follows a Dehnen profile(Dehnen 1993) with γ = 3 /
2, whereas the dark matterhalo is modeled by a Hernquist ( γ = 1) profile (Hern-quist 1990). The multi-component initial conditionsare generated using the distribution function method(e.g. Merritt 1985) following the approach of Hilz et al.(2012). In this method the distribution functions f i forthe different mass components are obtained from thecorresponding density profile ρ i and the total gravita-tional potential Φ T is then derived using Eddington’sformula (Binney & Tremaine 2008).In this paper we primarily study the orbits and merg-ing of SMBHs in unequal-mass mergers. The primarygalaxy is identical to the ‘ γ -1.5-BH-6’ initial conditionof Rantala et al. (2018, 2019), which was used to modelthe major merger progenitor of NGC 1600 (Thomaset al. 2016). This model corresponds to a massive gas-free early-type galaxy with a stellar mass of M (cid:63) =4 . × M (cid:12) with an effective radius of R e = 7 kpc Ws from SMBHs in Galactic Simulations Table 1.
Summary of the studied merger runs.Run m • m • M • M DM M (cid:63) N DM N (cid:63) q a e t t merger (10 M (cid:12) ) (10 M (cid:12) ) (10 M (cid:12) ) (10 M (cid:12) ) (10 M (cid:12) ) ( × ) ( × ) (pc) (Gyr) (Gyr)A 8.5 1.7 1.02 9.0 49.8 12.0 4.98 5:1 4.93 0.958 0.979 1.30B 10.2 1.7 1.19 10.5 58.1 14.0 5.81 6:1 5.45 0.971 1.09 1.22C 11.9 1.7 1.36 12.0 66.4 16.0 6.64 7:1 4.68 0.961 1.37 1.50D 13.6 1.7 1.53 13.5 74.7 18.0 7.47 8:1 5.32 0.954 2.10 2.38X 0.4 0.4 0.08 1.942 16.82 19.42 1.682 1:1 0.520 0.925 0.839 1.05 Note — Listed are the masses of the individual SMBHs ( m • , m • ), the total binary mass ( M • ), the total dark matterand stellar masses and particle numbers of the merger remnants ( M DM ,(cid:63) , N DM ,(cid:63) ), the merger mass ratio q , the SMBHbinary orbital parameters a and e at the time t , as well as the time of the SMBH binary merger t merger . The times aremeasured from the start of the simulation, and the time period considered in this work spans the interval ( t , t merger ). and a dark matter halo mass of M DM = 7 . × M (cid:12) .The dark matter fraction within the effective radius isset to f DM ( R e ) = 0 .
25 and the central SMBH has amass of M • = 8 . × M (cid:12) . The secondary galaxyis a scaled down version of the primary galaxy, withthe masses of all components divided by a factor of5, i.e. M (cid:63) = 8 . × M (cid:12) , M DM = 1 . × M (cid:12) , M • = 1 . × M (cid:12) and a resulting effective radius of R e = 3 . m (cid:63) = 10 M (cid:12) and m DM = 7 . × M (cid:12) for the dark matter particles resulting in N (cid:63) = 4 . × stellar and N DM = 1 . × dark matter particles for theprimary galaxy, the number of particles being a factorof five lower for the secondary galaxy.We first simulate a 5:1 minor merger between the pri-mary galaxy (Simulation A in Table 1). Next we con-tinue the 5:1 merger run with subsequent merger gen-erations until a total of 4 minor mergers are completed(Simulations B–D in Table 1). The fifth merger genera-tion simulation could not be used for the present study,as the black holes did not merge within a reasonablesimulation time due to the formation of a very low den-sity core in this system, with the missing stellar masseffectively giving rise to a final-parsec problem.Finally, in addition to the unequal-mass mergers wesimulated an equal-mass major merger between twoidentical galaxies with significantly lower masses of M (cid:63) = 8 . × M (cid:12) , M DM = 9 . × M (cid:12) , M • =4 × M (cid:12) and an effective radius of R e = 4 . N (cid:63) = 8 . × stellar and N DM = 9 . × dark matter particles foreach galaxy, resulting in particle masses of m (cid:63) = 10 M (cid:12) and m DM = 10 M (cid:12) for the stellar and dark matter com-ponents, respectively and a final merged SMBH mass of M • = 8 × M (cid:12) (simulation X in Table 1). The moti-vation for this run was to study the effect of the stellarenvironment on the SMBH merging process in a settingwhere the black hole masses are lower by a factor of ∼ r p ∼ . R e of the primary galaxy.After each minor merger (runs A–D), the merger rem-nant is reoriented so that the satellite galaxies fall infrom random directions with respect to the principal axisof the primary galaxy. In Figure 1 we show a sequence ofillustrative snapshots of run A, depicting the evolutionof the SMBH binary from the dynamical friction phase(left panel), through the stellar scattering stage (mid-dle panel) and into the final gravitational wave driveninspiral phase (right panel).2.4. Modeling the Final Phases of SMBH Inspirals
Following Rantala et al. (2018, 2019) we set theGADGET-3 integrator error tolerance to η = 0 .
002 andthe force accuracy to α = 0 . r chain = 10 pc and the perturber ra-dius to twice this value, i.e. r pert = 2 r chain , for allthe simulation runs. The gravitational softening lengthsare set to (cid:15) (cid:63) = 3 . (cid:15) DM = 100 pc, respectively.The softening lengths are chosen to fulfill the crite-rion r chain > . × (cid:15) (Rantala et al. 2017). In interac-tions between particles with different softening lengthsGADGET-3 uses the larger softening length.The mergers are initially evolved using Ketju with achosen GBS tolerance of η GBS = 10 − until the binaryblack hole semi-major axis is around a ∼ R s . Thiscorresponds to about a ∼ a ∼ . e (cid:38) .
9. Theseseparations are typically reached within t ∼ Ketju runs witha stricter GBS tolerance value of η GBS = 10 − . Theincreased accuracy is required for calculating the re-sulting GW spectrum and resolving in detail the post-Newtonian orbital motions of the individual SMBHs allthe way down to their final coalescence, which typicallyoccurs ∆ t merger ∼ Ketju runs that include the influence of thestellar background for our GW calculations until theseparation of the SMBHs is (cid:46) R s . This distance cor-responds to ∼ . ∼ × − pcfor run X. The final part of the merger is then runseparately without the surrounding stellar particles at amuch higher time resolution to facilitate direct waveformcomputations as described in section 4. The transitionseparation of ∼ R s was chosen as it provides a goodbalance between computational cost and the accuracyof these calculations. While in principle stellar parti-cles could also be included in this phase, we find thatignoring the environment in this very final stage doesnot cause significant changes in the SMBH evolution, asonly very few stellar particles are close enough to in-fluence the SMBH binary at our mass resolution. Thishas also been checked by comparing the evolution of theisolated binaries to the results from the full Ketju rundown to separations of ∼ R s , where the lower timeresolution Ketju data are still available. SMBH BINARY ORBIT ANALYSIS METHODS3.1.
Post-Newtonian Orbital Elements
For an analysis of the evolution of the SMBH binaryorbit, it is necessary to have a parametrization of theorbit that changes slowly and can be recovered from theinstantaneous position and velocity of the binary. Such aparametrization allows also computing the gravitationalwave emission of the system from only sparsely storeddata points, significantly reducing the computational ef-fort.In Newtonian gravity the orbit of a binary system canbe described by the Keplerian orbital elements includingthe semi-major axis a and the eccentricity e , which areconstant or evolve only slowly due to external perturba-tions. However, due to the inclusion of PN-correctionsto the motion of the SMBH binary, the standard Kep- lerian elements are no longer even approximately con-stant over an orbit, but rather oscillate quite stronglyespecially near the pericenter of an eccentric orbit, as isillustrated in Figure 2.As a result of this, the Keplerian parameters com-puted from sparsely stored data suffer from essentiallyrandom fluctuations due to the varying phase of the or-bit at the data points, which is problematic both foranalyzing the evolution of the orbit as well as for GWcomputations. For small binary separations, the Keple-rian orbital parameters are also no longer related to thegeometry of the orbit, as is apparent from the increaseof the Keplerian eccentricity for a nearly circular binaryin the right panel of Figure 2.The problems associated with the Keplerian parame-ters can be mostly remedied by using a quasi-Keplerianparametrization which accounts for the effects of thePN-corrections. Several such parametrizations accurateto different orders have been developed, and in thiswork we utilize the PN3 accurate quasi-Keplerian or-bital parametrization of Memmesheimer et al. (2004).This parametrization gives an approximate solution tothe conservative part of the PN3 accurate equations ofmotion of a non-spinning binary in the form R = a (1 − e R cos u ) , (3) n ( t − t ) = u − e t sin u + (cid:16) g t c + g t c (cid:17) ( v − u )+ (cid:18) f t c + f t c (cid:19) sin v + i t c sin 2 v + h t c sin 3 v, (4)2 π Φ ( φ − φ ) = v + (cid:18) f φ c + f φ c (cid:19) sin 2 v + (cid:16) g φ c + g φ c (cid:17) sin 3 v + i φ c sin 4 v + h φ c sin 5 v, (5)tan v (cid:115) e φ − e φ tan u . (6)The coordinates R and φ are the relative separationand angle in the orbital plane, with φ fixing the direc-tion of pericenter at time t = t . As in the Kepleriancase, the orbit is an ellipse with semi-major axis a , butnow the frequency of radial oscillations f r = n/ π is notthe same as the frequency of angular motion. Rather,during a single radial period the direction angle φ tra-verses an angle Φ > π leading to the precession of the Ws from SMBHs in Galactic Simulations t/T mean . . . . . . . . e e K e t e R e φ e g . . . t/T mean . . . . . . . . . e . . . . Figure 2.
Comparison of the different eccentricity definitions for an isolated highly eccentric binary with a semi-major axis of a ≈ R s evolved for four orbital periods using Ketju (left panel) and for a low eccentricity binary evolved from a ≈ R s down to a ≈ R s over ∼ T mean over the considered part of the orbital evolution. The eccentricities shown are the Keplerian eccentricity e K computed fromthe position and velocity using standard non-relativistic formulae; the quasi-Keplerian eccentricities e t , e R , and e φ ; and thegeometric eccentricity e g of equation (9). All the line widths are equal, with the thicker appearance of some lines caused byoscillations occurring during a single orbit. In the left panel the geometric eccentricity overlaps with e R and e φ , which differslightly from e t . The poor behavior of the Keplerian eccentricity is evident in both cases. orbit. The angles u and v are generalizations of the ec-centric and true anomaly, although their interpretationin terms of the geometry of the orbit is not as straight-forward as in the Keplerian case due to the appearanceof the additional series expansion factors g, f, i and h in the generalization of Kepler’s equation (4) and inequation (5). In addition the single eccentricity of theKeplerian parametrization is replaced by three differenteccentricities e R , e t and e φ . In the limit where the post-Newtonian corrections are negligible this parametriza-tion reduces to the standard Keplerian parametrization.The orbital parameters appearing in the aboveparametrization can be computed from the binary sep-aration vector R and relative velocity V in the modi-fied harmonic gauge used by Ketju via PN3 accurateapproximations of the conserved energy E and angularmomentum J in the form E = µ (cid:0) E + c − E + c − E + c − E (cid:1) (7) J = µ | R × V | (cid:0) J + c − J + c − J + c − J (cid:1) , (8)where µ = m m /M is the reduced mass and M = m + m the total mass of the binary. The expressionsfor the E i , J i in terms of R and V are quite lengthy, andcan be found in Memmesheimer et al. (2004). These ap-proximations of the conserved energy and angular mo-mentum are not constant along an orbit even when theradiation reaction terms are not included. Instead they oscillate slightly, and this oscillatory behavior carriesover to the orbital parameters, which are themselvesgiven as PN3 accurate expressions in terms of the ap-proximated E and J . The inclusion of the radiationreaction terms also introduces additional periodic oscil-lations (K¨onigsd¨orffer & Gopakumar 2006), but all ofthe oscillations are typically much smaller than the re-sulting secular evolution, and thus of no significant con-sequence.The appearance of three different eccentricities causessome ambiguity in their interpretation, as they do notcorrespond to a single orbital ellipse as in the Kep-lerian case. The radial eccentricity e R appears to bemost directly related to the dimensions of the orbit, butunfortunately the PN3 formulae for it and e φ fail fornearly circular relativistic orbits, resulting in e R < e φ <
0. The time eccentricity e t remains real, butit shows a small increase in this regime which does notcorrespond to an actual increased eccentricity of the or-bit. A measure of eccentricity that describes the actualgeometry of the orbit even in the strongly relativisticregime can be defined as (e.g. Csizmadia et al. 2012) e g = R max − R min R max + R min , (9)where R min , R max are the peri- and apocenter separa-tions of the binary. This definition is impractical sinceit cannot be computed from the instantaneous relativeposition and velocity of the binary, but it serves as auseful comparison.The right panel of Figure 2 shows the near-mergerbehavior of the different eccentricities for an isolated bi-nary with an eccentricity that is comparable to the onesseen in the main Ketju runs at this separation. Thetime eccentricity e t can be seen to be quite close to thegeometric eccentricity e g for the most part, while both e R and e φ go to zero too early, which is caused by thembecoming imaginary at that point. In regimes where thePN expansion parameter (cid:15) PN is small the differences be-tween the different eccentricities are negligible, so in thefollowing we elect to use e ≡ e t as the single parameterdescribing the eccentricity of the orbit. However, oneshould keep in mind that for the final stages of the bi-nary inspiral the eccentricity defined this way does notexactly describe the shape of the orbit. The other or-bital elements likely suffer from similar deterioration inthe very final stages of the inspiral, but for the vast ma-jority of the binary orbital evolution they work reliably.3.2. Extracting the Effects of the Environment
The stellar environment affects the evolution of thebinary orbit in addition to the evolution caused by GWemission. In the final stages of the binary evolutionconsidered in this work these effects are however fairlysmall, and to analyze them it is necessary to separatethem from the main GW-driven orbital evolution. Thetime derivative of a given quantity X , which can be e.g.the energy or the eccentricity, can be split into a partcorresponding to the effect from the PN forces includingthe GW reaction and a part corresponding to additionaleffects caused by interactions with the stellar environ-ment: d X d t = d X d t (cid:12)(cid:12)(cid:12)(cid:12) PN + d X d t (cid:12)(cid:12)(cid:12)(cid:12) env . (10)From this we find the cumulative effect of the environ-ment as∆ env X ( t ) = (cid:90) tt d X d t (cid:12)(cid:12)(cid:12)(cid:12) env d t = X ( t ) − (cid:90) tt d X d t (cid:12)(cid:12)(cid:12)(cid:12) PN d t. (11)This quantity is more robust than the instantaneousderivative d X/ d t | env , as numerically differentiating the X ( t ) derived from a Ketju run to compute the instan-taneous derivative amplifies any numerical noise presentin the data, which can easily hide the small effects ofthe stellar environment.To compute the effect of the PN forces at a giventime t (cid:48) , i.e. d X/ d t | PN ( t (cid:48) ), we perform a short integra-tion of an isolated binary with initial conditions taken from the Ketju computed SMBH positions and veloc-ities. The integration is performed for about 30 orbitsof the binary, and the derivative at the middle of theintegration interval is estimated by fitting a line to thedensely sampled evolution of X over this period. Theinitial conditions are chosen so that the middle of theintegration interval coincides with the desired time t (cid:48) .We evaluate the derivatives at about 50000 pointsspaced evenly in the semi-major axis a over the consid-ered section of the Ketju results. The cumulative ef-fect (cid:82) tt d X/ d t | PN d t is then found by integrating the in-stantaneous derivatives numerically. With these choicesthis method of computing the instantaneous derivativescaused by the PN terms appears to work well until thebinary separation is under ∼ R s , after which it failsas the evolution of the parameters starts to become non-linear over the integration period.3.3. Semi-Analytic Isolated Evolution
For comparison against the SMBH binary evolutioncomputed with
Ketju we consider two semi-analyticmodels used to describe binaries in simulations that can-not resolve them.The first model is the widely used Peters (1964) result,which considers the secular evolution of the Keplerianorbital parameters of an isolated binary due to the lead-ing radiation reaction term at PN2.5 level. The binaryis assumed to be otherwise Keplerian, and the orbitalparameters a and e are assumed to change only slowly.The evolution of the semi-major axis a is then given byd a d t = − G m m Mc a (1 − e ) / (cid:18) e + 3796 e (cid:19) (12)and the eccentricity e evolution byd e d t = − G m m Mc a (1 − e ) / e (cid:18) e (cid:19) . (13)3.4. Semi-Analytic Scattering Model
The main effect of the stellar environment on the evo-lution of the SMBH binary is through the scattering ofindividual stars during close encounters. These scatter-ing events remove energy and angular momentum fromthe binary, decreasing the semi-major axis a and increas-ing the eccentricity e . The second semi-analytic modelconsidered here is a commonly used (e.g. Sesana 2010;Kelley et al. 2017b) model describing this binary hard-ening process (Hills & Fullerton 1980; Quinlan 1996),which gives the evolution of the semi-major axis andeccentricity as dd t (cid:18) a (cid:19) = Gρσ H (14)
Ws from SMBHs in Galactic Simulations e d t = − Ka − d a d t = K Gρσ Ha, (15)where ρ is the stellar density and σ the velocity disper-sion. Equation (14) can also be written in terms of theenergy of the binary asd E d t = − GM µ Gρσ H, (16)where M = m + m is the total and µ = m m /M thereduced mass of the binary. This form is more suitedfor our applications as the post-Newtonian correctionsmodify the relation between the energy E and semi-major axis a compared to the non-relativistic Kepleriancase.The constants H and K are usually determinedby three-body scattering experiments (Quinlan 1996;Sesana et al. 2006), but Sesana & Khan (2015) find thatsimilar values for the constants can also be recovered inN-body simulations when the stellar density and veloc-ity dispersion are evaluated at the influence radius of theSMBH binary. Traditionally, the influence radius r inf isdefined by the condition that the enclosed stellar massis twice the binary mass, i.e. that M ∗ ( r < r inf ) = 2 M • . GRAVITATIONAL WAVE COMPUTATIONS4.1.
The Gravitational Wave Background
We focus our gravitational wave computations onquantities relevant for observations with PTAs, whichcan detect the low frequency GWs emitted by the mas-sive systems studied in this work. The main observablequantity describing the stochastic GWB is the charac-teristic strain h c ( f ), which can be computed using eithersemi-analytic or Monte Carlo (MC) methods when thedistribution of sources and their GW emission is known.In the semi-analytic method, the characteristic strain h c ( f ) of the gravitational wave background at frequency f is given by (Phinney 2001; Enoki & Nagashima 2007) h c ( f o ) = 4 Gπc f (cid:90) d z n c d E GW d f (cid:12)(cid:12)(cid:12)(cid:12) f =(1+ z ) f o , (17)where n c is the comoving number density of sources andthe source distribution is assumed to be spatially uni-form. As the source density n c is determined by thelarge-scale evolution of the Universe, here we computeonly d E GW / d f , the GW spectral energy density (SED)of a source integrated over its lifetime, which is thusaffected by any processes speeding up or delaying themerger of the SMBH binaries. In the MC methods ofcomputing the characteristic strain (e.g. Sesana et al.2008), the characteristic strain is instead computed bysimply summing the instantaneous contributions from each source. Thus all the environmental effects are in-cluded in the source distribution which we do not con-sider here.For computing the instantaneous GW flux and theSED of an SMBH binary we use two different meth-ods depending on the regime: in the weakly relativisticregime we save computational effort by applying analyt-ically derived formulae to compute the GW spectrumfrom the orbital elements, while in the regime wherethe relativistic effects become significant we compute theGW signal directly from the motion of the BHs. Thissplit between methods allows us to take advantage ofthe full PN3.5 accurate dynamics used in Ketju whilesaving computational effort where possible without sac-rificing numerical accuracy.4.2.
Semi-Analytic Spectrum
When the post-Newtonian effects are small, it is pos-sible to compute the contribution to the SED to goodaccuracy using analytic formulae derived assuming aslowly varying Keplerian orbit.A classic and widely used result derived by Peters &Mathews (1963) is the power emitted in gravitationalwaves by a binary in the lowest order quadrupole ap-proximation. For an eccentric orbit with eccentricity e and orbital frequency f p , the GWs are emitted in all har-monics of the orbital frequency, with the n th harmonichaving the power P n = 32 G / m m c a (2 πf p ) g ( n, e ) . (18)The function g ( n, e ) is given by g ( n, e ) = n (cid:26)(cid:104) J n − ( ne ) − eJ n − ( ne )+ 2 n J n ( ne ) + 2 eJ n +1 ( ne ) − J n +2 ( ne ) (cid:105) + (1 − e ) [ J n − ( ne ) − J n ( ne ) + J n +2 ( ne )] + 43 n J n ( ne ) (cid:27) , (19)where J n is the Bessel function of the first kind of order n . Since in this regime the PN effects are small, we maysimply use the PN quasi-Keplerian parameters in placeof the Keplerian parameters, thus taking f p = f r .From the instantaneous GW power and the knownevolution of the orbital elements the SED can be com-puted as (Enoki & Nagashima 2007)d E GW d f = ∞ (cid:88) n =1 P n ( f p , e ) τ GW ( f p , e ) nf p (cid:12)(cid:12)(cid:12)(cid:12) f p = f/n , (20)where τ GW = f p d t d f p = d t d(ln f p ) (21)0is the gravitational wave timescale. Implementing thisformula for computations from numerical simulationoutputs is mostly straightforward: for each frequency f we truncate the sum at some sufficiently high num-ber of harmonics, and interpolate the eccentricity e andgravitational wave timescale τ GW to the required orbitalfrequencies f p . Here e and τ GW are treated as functionsof the orbital frequency f p instead of time t , which ispossible as generally f p ( t ) is monotonic in the relevantphase of the binary evolution. The number of harmon-ics used for the computation is taken to be 750 in thiswork, which gives converged results for the eccentricitiesencountered.The computation of τ GW from the numerical data re-quires some care, as naively numerically differentiatingthe time-frequency data will amplify any noise causedfor example by interactions with the environment to un-usable levels. To avoid this, we instead compute thetimescale by integrating the time spent in logarithmicbins of orbital frequency. If two adjacent data pointsbelong to the same bin the whole interval is assigned tothat bin, while for points in different bins the time inter-val is divided over the spanned frequency interval. Wefind that using about 50 bins per decade in frequencygives results that match well with the differentiation ofsmooth data.The error in applying these leading order or PN0 re-sults is expected to be of order PN1, and in subsec-tion 4.4 we find that the error between this method andthe more accurate method described in the next sub-section is of the order of a few percent when the PNexpansion parameter (cid:15) PN ≈ .
01, in line with expecta-tions.Similar methods are commonly extended also beyondthe regime where PN effects are minor all the way downto the final orbits before the BHs merge (e.g. Berentzenet al. 2009; Kelley et al. 2017b). This may work wellif the orbital evolution is consistently performed in thesame approximation, but coupled with PN3.5 accurateorbits this leads to clearly erroneous results. There arealso analytic results for PN accurate spectra available inthe literature (e.g. Tessmer & Sch¨afer 2010; Klein et al.2018), however these are quite cumbersome to apply,and we find the methods of the next subsection to bemore convenient.4.3.
Direct Waveform Calculation
When the effects of the PN-corrections on the nownon-Keplerian binary orbit become significant, themethods of the previous subsection are no longer reli-able. In this regime we instead directly compute the emitted waveform, and use a Fourier transform to com-pute the spectrum.In principle this method could be used over the entirebinary evolution, but for most part the semi-analyticmethod is sufficiently accurate. Further, if the binaryeccentricity is high, the timestep required to resolve thehighest harmonics of the orbital frequency becomes verysmall resulting in considerable computational cost. Forthese reasons we limit this method to the short periodbefore the BH binary merger where it is necessary foraccurate results. In this regime the effect of the envi-ronment is also negligible, allowing us to easily computea densely outputted evolution of the SMBH binary inisolation using the
Ketju integrator. Inclusion of theenvironment at this stage would also be possible, butthis would require suitable schemes to resample the stel-lar population to lower mass particles in order to avoidartefacts in the spectrum caused by spurious unphysi-cally strong interactions from massive stellar particles.In computing the GW waveform emitted by theSMBH binary we follow Poisson & Will (2014). Theleading contribution to the strain tensor h ij is thequadrupole term h ij = 4 Gm m c dM (cid:16) V i V j − GM R i R j R (cid:17) , (22)where the relative velocity V i and separation R i are com-puted directly from the numerical integration of SMBHtrajectories. The observer distance d only scales theresult, and can in practice be ignored in these computa-tions. The radiation reaction force corresponding to thispart appears as the a PN2.5 term in the equations of mo-tion, and to maintain consistency with the PN3.5 accu-rate equations of motion we also include PN-correctionsto the waveform up to PN1 order beyond the leadingPN0 contribution of equation (22). The formulae forthese corrections for a binary system are lengthy, andare not reproduced here but can be found in Poisson &Will (2014). Note that by convention the waveform PNorders are offset from the corresponding effects in theequations of motion by 2 . h obtained as a function of time isthen projected into the physically observed polarizations h + and h × , which are related to emitted GW power asd E GW d t = d c πG (cid:90) (cid:0) ( ∂ τ h + ) + ( ∂ τ h × ) (cid:1) dΩ , (23)where τ = t − d/c is the retarded time, which can here bereplaced with the coordinate time t due to the arbitraryobserver distance. In terms of the Fourier coefficients˜ h + , × ( f ) the time dependent signals can be written as h + , × ( t ) = (cid:88) k √ (cid:12)(cid:12)(cid:12) ˜ h + , × ( f k ) (cid:12)(cid:12)(cid:12) sin(2 πf k t + φ k ) , (24) Ws from SMBHs in Galactic Simulations .
00 2 .
05 2 .
10 2 .
15 2 . f/f . . . . . . . d E G W / d f
102 3 5 7 20 30 f/f . . . . . . . . d E G W / d f PN3.5 + PN0PN3.5 + PN1PN3.5 + SAPN3.0 + PN0PN3.0 + SA
Figure 3.
The emitted GW energy spectrum (in code units) computed for PN3.5 accurate binary motion using the semi-analytic method (PN3.5 + SA, dashed green line) and directly from the waveform using either only the leading quadrupoleterms (PN3.5 + PN0, solid blue line) or including also PN1 corrections to the waveform (PN3.5 + PN1, solid orange line).Only the smoothed result is shown for the direct waveform calculations, but some artefacts caused by the finite signal lengthare still visible near the edges. The left panel shows the main n = 2 harmonic of the signal from a SMBH binary with aninitial semi-major axis of a = 250 R s and eccentricity e = 0 . × orbits. The right panel shows the fullspectrum from a binary with initial a = 40 R s and e = 0 . a ≈ R s . The frequency is given in units of the initialradial frequency f = f r ( t = 0). The right panel also shows in addition the smoothed quadrupole spectrum computed fromthe waveform (PN3.0 + PN0, solid red line) and using the semi-analytic method (PN3.0 + SA, dot-dashed cyan line) from anotherwise identical binary evolution but using only PN-corrections up to PN3 level. Note the logarithmic frequency scale in theright panel. so the total energy emitted over a time period of length T can be written as∆ E GW = (cid:90) t + Tt d t d E GW d t = d c T πG (cid:88) kp =+ , × (2 πf k ) (cid:90) (cid:12)(cid:12)(cid:12) ˜ h p ( f k ) (cid:12)(cid:12)(cid:12) dΩ . (25)As this can also be written as∆ E GW = (cid:88) k ∆ E GW ( f k )∆ f ∆ f, (26)the SED at the frequency f k is given by the correspond-ing term in the sumd E GW ( f k )d f ≈ ∆ E GW ( f k )∆ f = c (2 πf k T d ) G (cid:88) p =+ , × (cid:90) (cid:12)(cid:12)(cid:12) ˜ h p ( f k ) (cid:12)(cid:12)(cid:12) dΩ , (27)where ∆ f = f k +1 − f k = 1 /T is the frequency resolutionof a discrete Fourier transform of a signal with evenlyspaced samples.To numerically compute these quantities we use thefast Fourier transform (FFT) algorithm, which allows computing the coefficients ˜ h + , × efficiently from evenlyspaced samples of the waveform. The integral over thesolid angle Ω is conveniently performed numerically, al-though at least when ignoring the direction dependentPN corrections to the waveform it can also be treatedanalytically. The changing frequency of the signal leadsto rapid oscillations in the FFT results, which we sup-press by averaging them over a small window. Addi-tionally, the abrupt changes at the ends of the signallead to some artefacts. These could be reduced by win-dowing, but this would also reduce the accuracy of thecomputed energy of the GW signal. From the Fouriertransformed signals we may also evaluate the instanta-neous power from equation (23), using the relation be-tween derivatives and Fourier transforms. Evaluatingthe time derivatives via Fourier transforms is also moreaccurate for periodic signals than using finite differencetechniques.4.4. Accuracy of the Semi-Analytic Method
In order to confirm that the semi-analytic methodgives good results in the mildly relativistic regime and tofind the limits of its applicability we compare it to thedirect waveform calculation over short sections of the2evolution of an isolated binary. This allows also evaluat-ing the effect of including the additional PN correctionsto the waveform instead of simply using the quadrupoleformula (22).For the mildly relativistic regime we compute the orbitof an isolated SMBH binary with a mass ratio of q = 5 :1 and initial semi-major axis a = 250 R s and eccentricity e = 0 . × orbital periods. This fairly largenumber of orbits allows the orbital frequency to evolveenough for the GW spectrum to be resolved instead ofeffectively consisting of δ -spikes which are problematicnumerically. The post-Newtonian expansion parameteris approximately (cid:15) PN ≈ /
200 in this case, so the effectof the additional PN effects not taken into account inthe semi-analytic model can be expected to be of theorder of half a percent.The left panel of Figure 3 shows the resulting spec-trum around the strongest n = 2 harmonic. Apart fromthe small shift in frequency, which is of the order of (cid:15) PN , the semi-analytic result matches well with the di-rect waveform computation using the full PN1 accuratewaveform. This agreement with the PN1 waveform mayseem surprising, as the semi-analytic method only in-cludes the PN0 quadrupole radiation, but this can besimply explained based on energy conservation. In thisregime the binary is still close to Newtonian and themain effect of the PN1 waveform corrections is to main-tain energy balance with the PN3.5 radiation reactionterms. As a result, the semi-analytic method, which isalso constructed around energy conservation of a Newto-nian binary, gives matching results. The frequency shiftis caused by the precession of the binary, which leads tothe GW frequencies differing from integer multiples ofthe radial frequency used for the calculation. While thisshift could in principle be corrected for, we find that itis insignificant when the spectrum is computed over theentire binary evolution.To evaluate the methods in the highly relativisticregime we compute the evolution of the binary from aninitial semi-major axis of a = 40 R s down to a = 6 R s ,where the PN equations of motion become unreliable.The initial eccentricity is set to e ≈ .
1, which is com-parable to the eccentricities found in our
Ketju runs atthis binary separation. The resulting spectra in the rightpanel of Figure 3 clearly show the increasing significanceof consistently computing the waveform from the actualbinary motion as well as including the PN-corrections tothe waveform, with the differences between the methodsgrowing to tens of percents towards the end of the sim-ulation.The error of the semi-analytic method compared tothe full PN waveform is about 5% at the point where all the major harmonics are visible in the spectrum, butsomewhat lower when compared to the PN0 waveform.In this regime it is also of interest to look at the sig-nificance of the higher order radiation reaction terms atPN3.5 level. The right panel of Figure 3 includes also thequadrupole spectrum from an otherwise identical binaryevolved using PN-corrections only up to PN3 level. Theamplitude of the spectrum is clearly different from thefull PN3.5 result, showing the importance of includingthe higher order radiation reaction term. Interestingly,the semi-analytic method yields essentially exact resultscompared to the PN0 waveform in this case. This in-dicates that the differences between the semi-analyticmethod and the direct waveform methods are at leastpartly due to the PN3.5 radiation reaction term causingsignificant changes in the motion of the binary. RESULTS5.1.
Orbital Evolution
We begin by studying the evolution of the orbital ele-ments of the SMBH binaries in the
Ketju runs and com-pare them to the parameters predicted by semi-analyticmodels describing both isolated binaries and binaries in-teracting with a population of stars. The semi-analyticcomparison models are started from the initial condi-tions at the start of the analyzed section, i.e. using thevalues of a and e listed in Table 1. This is similar towhat would be done in a standard softened simulationwhen the SMBH binary separation has decreased to theorder of a gravitational softening length.Isolated binaries are modeled using the model de-scribed in subsection 3.3, hereafter the Peters model,which describes a Keplerian binary perturbed by theleading gravitational radiation reaction term.The interaction with stars is included by combiningthe Peters model with the scattering model of subsec-tion 3.4, to produce what we call the
Peters-Quinlan model. This model requires as inputs the stellar den-sity and velocity dispersion, which are computed fromthe
Ketju run results in a sphere within the influenceradius of the binary. For these calculations we use thevalues of H and K parameters determined by the fit-ting formulae of Sesana et al. (2006), which allows usto see how much the behavior of the Ketju runs differsfor typically used literature values. The values of H arecomputed using the fit for a circular binary, as Sesanaet al. (2006) do not tabulate the coefficients for the H fit for eccentric binaries. Their results however indicatethat at higher eccentricities the value of H should tendto increase slightly.The time evolution of the semi-major axis of the orbitsin the Ketju runs and the semi-analytic comparison
Ws from SMBHs in Galactic Simulations − . − . − . . . . t (Gyr)10 − − a ( p c ) ABCDX
Ketju
Peters − . − . − . . . . t (Gyr) Ketju
Peters-Quinlan
Figure 4.
Evolution of the semi-major axis a as a function of time. The solid lines show the Ketju runs while the dashed linesshow the Peters (left panel) and Peters-Quinlan (right panel) semi-analytic comparison models. The times have been shifted sothat the merger of the BHs in the
Ketju runs takes place at t = 0, corresponding to a shift of 1 Gyr to 2 . − − − − f r (yr − )0 . . . . . . e ABCDX − . − . . ∆ e Peters10 − − − − f r (yr − )0 . . . . ∆ e Peters-Quinlan
Figure 5.
Eccentricity evolution as function of radial orbital frequency f r . The left panel shows the results from the Ketju runs, and the right panel shows the difference to the
Ketju runs for the Peters model (top right) and the Peters-Quinlan model(bottom right). models is shown in Figure 4. The evolution of the bi-nary is mainly driven by GW emission at this stage, butthe environment still has a non-negligible effect. Thiscan be seen by studying the isolated Peters models, inwhich the SMBHs merge significantly later than in the
Ketju runs. On the other hand, the Peters-Quinlanmodels merge significantly earlier than the
Ketju runs,which is caused by the fact that the literature valuesof the H and K parameters from Sesana et al. (2006) significantly overestimates the effects of stellar scatter-ing for our Ketju runs with depleted stellar cores (runsA–D). This is most obvious in the case of run A, whilein the case of run X the match between the
Ketju runand the Peters-Quinlan model is reasonably good.The eccentricity evolution of the runs is shown inFigure 5 as a function of the orbital frequency f r ≈ (cid:112) GM/a / (2 π ), as this allows a comparison between themodels with different merger times and also relates the4 Table 2.
Fitted and literature values of the hardening pa-rametersRun H Ketju K Ketju H S K S Gρ/σ (10 − pc − yr − )A 2.18 0.032 16.2 0.056 5.49B 4.38 0.049 16.3 0.059 2.92C 5.98 0.052 16.4 0.068 1.68D 8.44 0.027 16.6 0.068 1.03X 11.5 0.098 14.4 0.054 24.4 Note — The values of the hardening parameters H Ketju and K Ketju of subsection 3.4 have been computed by linear fitsto the
Ketju data with the effects of the PN3.5 level cor-rections subtracted (see subsection 3.2). The values usedfor the Peters-Quinlan comparison models ( H S , K S ) com-puted from the fitting formulae of Sesana et al. (2006) arealso listed, as well as the value of Gρ/σ measured from thesimulations at the influence radii of the SMBH binaries. eccentricity to the frequency of emitted GWs. RunsA–D all show fairly similar behavior with high eccen-tricities in the range e ≈ . .
98 at the initial semi-major axis a ∼ f r ∼ − yr − where we begin consideringthese runs. Run X has a slightly lower eccentricity atits initial a ≈ . . Ketju runs.As can be seen from equation (12), the dependence of theevolution on the eccentricity is highly non-linear, and es-pecially at high eccentricities even small differences leadto large effects in the merger time. Notably the differ-ence in the evolution of the eccentricity is smallest inthe case of the Peters-Quinlan model run X and largestin the case of run A, which show also correspondinglythe smallest and largest difference in the SMBH mergertimes, respectively.5.2.
Strength of Environmental Effects
To quantify the strength of the environmental effects,we use the method of subsection 3.2 to subtract thedominant effects of the PN-corrections and recover theresidual effect caused directly by interactions with thestellar environment. The results for the energy and ec-centricity of the SMBH binaries are shown in Figure 6.The energies show a clear linear evolution with time, in agreement with the semi-analytic model equation (16).The eccentricities also show a secular evolution with thecorrect sign, but the effect is fairly weak.We perform linear fits to the data in Figure 6 to de-termine the values of the H and K parameters of thescattering model that would correctly describe the bi-nary evolution when the stellar densities and velocitydispersions are derived from the Ketju runs. Note thatwhile equation (16) implies that the change in energy islinear in time, equation (15) implies that the change iseccentricity is linear in the integrated semi-major axis:∆ env e ( t ) = KH Gρσ (cid:90) tt a d t. (28)The results are given in Table 2, which also shows the H S - and K S -values used for the Peters-Quinlan compar-ison models, derived from Sesana et al. (2006). In thecase of run X the fitted value of H is close to the valuedetermined from the Sesana et al. (2006) fitting formu-lae, being about 25% lower. This is comparable to theresults of Sesana & Khan (2015), who found that N-body simulations tend to produce approximately 30%lower values of H than those found from idealized 3-body scattering experiments, which were used by Sesanaet al. (2006).For runs A–D the agreement between the fitted H pa-rameter and the one used for the Peters-Quinlan modelis much poorer, with up to a factor of eight differencein the case of run A. This is not entirely unexpected asthese runs describe massive early-type galaxies with verylarge low density stellar cores. Thus, the derived smallvalues of H are caused by the fairly empty loss conesin these runs, as most of the stars on radial trajectorieshave already been ejected at this stage by the strongcore scouring effects of the SMBH binary (Rantala et al.2019).The fitted values of K are slightly lower for runs A–D,as can be expected due to the low density cores. How-ever, the differences are not quite as large as in the caseof the H parameter. On the other hand, the fitted valueof K for run X is almost twice as large as the literaturevalue used for comparison. This is likely partly due toeffect of eccentricity on the parameters, as the compar-ison value was computed for a circular binary. For runX the eccentricity evolution between the Peters-Quinlanmodel and the Ketju run agrees very well, despite slightdifferences in the scattering parameters (see Table 2).This is made possible since there is a certain degree ofdegeneracy between the H and K parameters in thePeters-Quinlan model, allowing very similar evolutionsto be produced for a range of parameter values. Ws from SMBHs in Galactic Simulations − . − . − . . t (Gyr) − . − . − . − . − . − . − . . ∆ e n v E / µ c × − ABCDX − . − . − . . t (Gyr)0 . . . . . . . . ∆ e n v e × − Ketju
Residual error
Figure 6.
The evolution of the energy per unit reduced mass
E/µc (left) and the eccentricity e (right) of the SMBH binarydue to the interaction with the stellar environment. The contribution from the PN terms has been removed as explained insection 3.2, with the dashed lines showing the residual from the subtraction process for an isolated evolution. The removalprocess works well until the final stages of the merger where the difference in energy starts to grow rapidly. There is a remainingsecular evolution in the Ketju energy which appears to be consistent with the simple semi-analytic scattering models.
We also computed comparison Peters-Quinlan modelsusing the parameter values fitted from the
Ketju runs.For all runs the results are accurate at a similar levelas the run X comparison model using literature valuesdiscussed above, i.e. at the level of a few percent. As thebehavior is so similar, we opt to not display the resultsof these calculations.5.3.
Gravitational Wave SEDs
Next we focus on the GW signal emitted by the SMBHbinaries, and the differences between the signal com-puted from the
Ketju runs and the semi-analytic com-parison models. In computing the GW spectrum fromthe
Ketju runs we use the semi-analytic method ofsubsection 4.2 until the SMBH binary semi-major axisis ∼ R s , after which we switch to the direct wave-form method of subsection 4.3 and compute the isolatedevolution of the binary for the final a (cid:46) R s . Theexact point where the switch happens depends on theeccentricity of the binary, as all the major harmonicsof the initial orbital frequency need to be visible in thespectrum computed from the waveform. Based on sub-section 4.4, the error of the semi-analytic GW methodshould be of the order of a few percent at the switchoverpoint, as the error scales with (cid:15) PN ∝ a − . For computingthe GW signal from the semi-analytic Peters and Peters-Quinlan models, we naturally utilize the semi-analyticmethod throughout.The energy spectrum of gravitational waves as emittedover the lifetime of the binary is shown in Figure 7. The spectra are all fairly similar due to the similar orbitalevolution, and show the characteristic peaked shape ofan initially highly eccentric binary. The differences be-tween the runs are mainly explained by the differenttotal black hole masses and mass ratios, as well as thedifferent eccentricities.At high frequencies, where the binaries have mostlycircularized ( e ∼ R = 6 R s where the Ketju simulationis stopped. The Peters-Quinlan model gives differencesof a similar magnitude for these runs, which is expectedas in this phase these massive binaries are highly rela-tivistic and the differences in the spectrum are mainlydue to PN effects. The effect of the environment in thespectrum is larger for run X, where the Peters modelshows a difference of almost 20% at f = 0 . − , whilethe Peters-Quinlan model agrees at the percent level.This is mainly due to the larger eccentricity of the bi-nary in this frequency range, and in particular the differ-ence in eccentricity in the Peters model, which shifts thepeak of the spectrum. The agreement at slightly higherfrequencies where the binary has circularized is muchbetter, with the semi-analytic models matching almostexactly the direct waveform computation results. At the6 − − − f (yr − )10 d E G W / d f ( M (cid:12) / y r − ) ABCDX 1234 d E G W / d f r a t i o Peters10 − − − f (yr − )0 . . . d E G W / d f r a t i o Peters-Quinlan 0 . . . . Figure 7.
Total emitted GW energy per unit frequency, i.e. the GW SED, for the
Ketju runs (left), and the ratios of thesemi-analytic Peters model (top right) and the semi-analytic Peters-Quinlan model (bottom right) results to the
Ketju result.The zoom-ins highlight the frequency range most relevant for PTA observations ( f (cid:38) − yr − ). − − − − f r (yr − )10 τ G W ( y r ) ABCDX 1 . . . . τ G W r a t i o Peters10 − − − − f r (yr − )0 . . . . τ G W r a t i o Peters-Quinlan 1 . . . . Figure 8.
The gravitational wave timescale τ GW = d t d ln f r of the Ketju runs (left), and the ratios of the semi-analyticPeters model (top right) and the semi-analytic Peters-Quinlan model (bottom right) results to the
Ketju result. The zoom-inshighlight the frequency range most relevant for PTA observations ( f (cid:38) − yr − ). point when the GW calculation is switched to the directwaveform method, the differences to the semi-analyticGW method are of the order of few percent, consistentwith the PN expansion parameter being (cid:15) PN ∼ .
01 atthat point. This leads to the small jump that can beseen in the difference to the semi-analytic comparisonmodels (see the zoom-ins in Figures 7 and 8).At lower frequencies where only the semi-analytic GWcomputations are performed on the
Ketju simulations, the isolated Peters models lead to larger amplitudes withpeaks at lower frequencies. The resulting relative differ-ences are mostly around 50%, but can reach values ofup to 300% in the case of run X. These large differencesare due to the faster hardening of the binary due to thestellar environment, which is not accounted for in thePeters models and in fact the Peters-Quinlan modelswhich include this effect are in much better agreementwith the
Ketju results. The most significant differences
Ws from SMBHs in Galactic Simulations
Ketju and Peters-Quinlan results occur inthe case of run A, where the Peters-Quinlan result isabout 70% lower. In this case the hardening parameterused for the Peters-Quinlan model was already found tobe 5 to 8 times larger than the effective ones determinedfrom the
Ketju run, so the large difference is not sur-prising. The other runs (B–D) with large differences inthe hardening parameters show similar behavior, withdifferences of around 40%. In any case the differencesoccur mainly at relatively low GW frequencies that areunobservable for the foreseeable future.The differences in the GW spectrum are mainly dueto differences in the gravitational wave timescale τ GW ,shown in Figure 8, and also in part due to the differencesin orbital eccentricity. The evolution of eccentricity de-termines how the emitted GWs are spread out over thespectrum, which also determines the location of the peakof the spectrum. Here too the effects of the environ-ment are obvious, with the Peters models having longerand the Peters-Quinlan models shorter GW timescalesthan the Ketju runs. The differences are mostly be-tween 10 and 60 percent at orbital frequencies below0 . − , with run X showing both the highest differenceof ∼ Ketju .5.4.
Gravitational waves and energy conservation
The good agreement at the percent level between thespectra computed using the Peters formulae and the fullPN computation at high frequencies may seem surpris-ing, since the strongest effects of the PN-correctionsshould appear there. If we look at the instantaneousGW flux, shown in Figure 9, the difference between themethods does indeed increase with increasing frequency,reaching tens of percent. However, this is counteredby the slightly slower evolution of the orbit accordingto the Peters formula, which can be seen as a slightlylonger timescale at high frequencies in τ GW in Figure 8.Fundamentally, the similarity of the spectra computedwith the different methods for an isolated binary fol-lows from the conservation of energy, which is exact byconstruction in the Peters evolution and approximatein the PN3.5 evolution, as was also noted in subsec-tion 4.4. Figure 9 also shows the importance of includ-ing the waveform PN1 corrections as well, as otherwisethe flux can be overestimated by over 10%.Another point relating to the energy conservation ofthe system when using PN dynamics to compute thebinary evolution and GW emission in the final stages − f r (yr − )0 . . . . . . . . . d E G W / d t r a t i o ABCDX
Figure 9.
The ratio of the instantaneous GW fluxes com-puted using less accurate approximations to the PN3.5 mo-tion + PN1 waveform used for the final part of the GW SEDcalculations. The solid lines show the result of the PN3.5motion + PN0 (quadrupole) waveform calculation, while thedashed lines show the result of the semi-analytic Peters for-mulae calculation. The differences between the runs are dueto the different mass ratios and slightly different eccentrici-ties. of an SMBH coalescence is the fact that the energy ofthe system is conserved only approximately. This lim-its the applicability of the method, as large violationsof energy conservation clearly signify incorrect results.Figure 10 shows the total relative error in the energyconservation of the system during the final stages of themerger. The orbital energy is computed using the PNaccurate energy of equation (7), while the emitted GWenergy is computed by integrating the instantaneous en-ergy flux computed from the waveform. For most partthe full computation including the PN waveform correc-tions gives an error of under 5%, and only in the veryfinal stages does the error begin to rapidly increase. Thissuggest that the results computed using this method arefairly reliable down to around a ∼ R s or to orbitalfrequencies of f r ∼ (10 M (cid:12) /M ) yr − , which also corre-sponds to the part in the spectrum in Figure 7 where thedirect waveform and the Peters results begin to again di-verge. When the GW emission is computed using onlythe quadrupole formula the error is naturally greater.A similar error in energy conservation for the full PNcomputation has also been observed by Csizmadia et al.(2012). DISCUSSION6.1.
SMBH eccentricity Distribution a ( R s ) − . . . . . . . . . . ( ∆ E + ∆ E G W ) / | E | ABCDX
Figure 10.
The relative error in total energy conservationduring the final isolated evolution of the binaries when theemitted GW energy is computed from the PN1 waveform(solid) or the PN0 quadrupole waveform (dashed).
Since the eccentricity affects strongly both the merg-ing time of the SMBHs (Peters 1964) and the emittedGW spectrum (Enoki & Nagashima 2007), understand-ing the distribution of eccentricities of SMBH binariesis essential for correctly predicting the emitted GWB.As can be seen from Table 1 all the
Ketju runs studiedin this paper show high initial eccentricities in excessof e > . Ketju (Rantala et al.2017, 2018).Previous works studying the evolution of SMBH bi-naries in a galaxy merger using a smaller number ofparticles have also in general found high binary eccen-tricities, e (cid:38) . M • (cid:46) M (cid:12) ).6.2. Relevance for Pulsar Timing Array Predictions
It is important to account for the SMBH binary ec-centricity when computing predictions for PTA observa-tions. In older work this was not typically done, and in-stead the binaries were assumed to be on circular orbits(e.g. Sesana et al. 2008; Kelley et al. 2017a). However,more recent treatments have accounted for this in thecase of a fixed initial eccentricity as well as other effectssuch as stellar scattering and drag from a circumbinaryaccretion disk using semi-analytic models (e.g. Kelleyet al. 2017b; Sesana et al. 2018).These semi-analytic models appear to work quite wellfor the mergers considered in this paper, as far as theprocesses modeled by
Ketju are concerned. The differ-ences in the evolution of orbital parameters (see Figures4 and 5) and the GW spectrum in the PTA frequencyrange (see Figure 7) are brought down to under 5% pro-vided a correct choice of semi-analytic model parame-ters, and even when using incorrect parameter values orignoring the stellar scattering completely the differencesremain typically at a level of few tens of percent.We find that the measured hardening parameters ofthe considered semi-analytic scattering model can besmaller for our
Ketju cored galaxy simulations thancommon literature values by factors up to 8. However,this is partially offset by the PN3.5 accurate equations ofmotion also speeding up the merger process slightly com-pared to models including only the leading quadrupoleradiation reaction. The difference to literature valuesis also smaller for the simulation without a stellar corethat has less extreme SMBH masses (run X), suggest-ing that for most SMBH binaries the commonly usedsemi-analytic models and their parameters are reason-ably accurate.Comparing the eccentricities seen in our simulationsto the results of Kelley et al. (2017b), we find that whileour eccentricities are among the highest considered therethey are not quite high enough to cause the strong sup-pression of the GWB seen for their model with an initialeccentricity of e ∼ . f ∼ . − , depending on the eccentricities of the bi-naries (Kelley et al. 2017b).Whether the error caused by the simplified semi-analytic models leads only to a slightly increased sta-tistical uncertainty, or if some kind of systematic bias Ws from SMBHs in Galactic Simulations
Ketju calculation was achieved when using theproperties of the stellar population derived directly fromthe
Ketju runs. Thus the strong effects of the SMBHbinary on the stellar population are accounted for in thedensity and velocity dispersion, while in the typical usecase of semi-analytic models they are not, as this regimefalls below the gravitational softening length. This issuewill particularly important for cosmological simulationsin which the stellar population cannot be resolved athigh spatial accuracy.At high gravitational wave frequencies affected bydiscretization effects, the effects of higher order PN-corrections can also become significant. This is due tothe fact that the GW spectrum is determined by the in-stantaneous GW flux when there are only a few sources,and as was seen the difference between the PN-accurateinstantaneous flux and the commonly used quadrupoleapproximation is of the order of 30% at f ∼ − for M ∼ M (cid:12) binaries (see Fig 9). Ignoring these effectscan thus lead to a systematic bias in the computed back-ground amplitude, if there are enough massive binariesfor which this effect is at its strongest. However, we notethat the corrections to the GW flux can also be handledanalytically to at least PN3 order (e.g. Arun et al. 2008;Blanchet 2014).6.3. Potential simulation caveats
The derived differences between the accurate
Ketju computations and the simpler semi-analytical modelsare relatively small, on the level of ∼ ∼ real stars.Secondly, the motion of the binary is modeled only ap-proximately using PN3.5 level corrections for the binaryonly, and while the accuracy of this modeling is stateof the art for the post-Newtonian approach, the actualphysics may differ from this.The hybrid nature of Ketju allows us to utilize alarger number of particles than in standard N-bodycodes. However, despite this fact, the actual numberof stellar particles interacting with the SMBH binary inits final stages before the merger is still small comparedto the actual number of stars that they represent. Whenthe SMBH binary merges in our simulations, there aretypically only a few stellar particles left within 10 pcof the binary instead of tens or hundreds of thousands. As a result, strong interactions between the binary andstars happen far more rarely, and when they do happenthey are too strong, leading to step-like changes in thebinary orbital parameters. Fortunately, it appears thatmost of the environmental effects are caused by fairlyweak long range interactions, thus the error in modelingthe strong interactions should not seriously affect thelong term binary evolution.The accuracy of the binary motion also depends onthe accuracy of the PN3.5 equations of motion. Herewe are limited by the fact that higher order equations ofmotion are unavailable for general binaries. Such higherorder corrections would be further complicated by theirdependence on the history of the binary in addition toits instantaneous state, making implementing them im-practical.While formally the higher order corrections should beof order O ( (cid:15) ), in practice the effects may be signifi-cant even for fairly small (cid:15) PN . This was illustrated bythe rather large change in the spectrum when going fromPN3 to PN3.5 equations of motion in Figure 3, eventhough formally the effect is of order (cid:15) . (cid:46) − there.The apparently fairly significant contributions from theunknown higher order corrections are also seen in thefailure of energy conservation when the binary is verynear merger (see Fig. 10). Since the eccentricity of thebinary is fairly small in this phase, results computedfor circular binaries can be used to estimate just howmuch the higher order corrections affect the evolutionof the binary. According to Blanchet et al. (1995), thenon-linear tail-terms that are the leading correction ne-glected here cause the number of orbits an equal-massbinary completes during its last decade in orbital fre-quency to change by almost 10%, which can be consid-ered to be a fairly substantial effect. CONCLUSIONSIn this paper we have studied the final stages of SMBHbinary evolution in simulations of mergers of gas-poormassive early-type galaxies using the hybrid
Ketju code(Rantala et al. 2017, 2018). In general, the evolution ofthe SMBH binaries was found to agree well with predic-tions from simple semi-analytic models commonly usedfor studying gravitational wave emission from mergingSMBHs in a cosmological setting.Differences in the emitted energy in gravitationalwaves between the
Ketju runs and the semi-analyticmodels were found to be of the order 10% in the fre-quency range of ( f (cid:38) − yr − ), which is accessibleby PTA. This relatively good agreement was reachedprovided that the properties of the stellar populationaccounting for the ejection of stars via interactions with0the SMBHs and the binary eccentricity was chosen cor-rectly. The agreement between the semi-analytic modelsand the Ketju runs can be further improved to the levelof a few percent by using parameter values derived fromthe
Ketju runs, showing that the model itself capturesthe relevant dynamics well.Overall, the derived parameters from the
Ketju runsfor the semi-analytic scattering model were found to besomewhat smaller than the literature values commonlyused (Sesana et al. 2006). However, this is largely ex-plained by the fact that our runs A–D describe very mas-sive early-type galaxies with large cores scoured almostempty by very massive SMBH binaries. The simulationwithout a stellar core and much smaller SMBHs (runX) is in much better agreement with the semi-analyticscattering model. Determining how the effects of theSMBH binary on the stellar population affect the scat-tering model parameters in general would allow reducingthe error of the semi-analytic models.The eccentricities of all of our simulated SMBH bina-ries were found to be very high, in excess of e > . a ∼ R s .Thus, based on our simulations the SMBHs enter theinspiral phase on strongly eccentric orbits, with circularorbits being viable only at the very final stages of the or-bit after the completion of a GW driven circularizationprocess (e.g. Sesana et al. 2008; Kelley et al. 2017a,b).In addition, the high eccentricities also result in themerger times of the binaries being quite sensitive tosmall changes in their initial parameters. This is in par- ticular the case for the semi-analytic models as in gen-eral the eccentricity has to be fixed at a relatively largebinary separation, often leading to differences of sometens of percents in the binary evolution when comparedto the resolved Ketju runs.Based on our sample of
Ketju mergers presented inthis paper and also in earlier work (Rantala et al. 2017,2018, 2019) it seems that our binaries do not reach ex-tremely high eccentricities in excess of e (cid:38) .
99, whichcould lead to significant suppression of the gravitationalwave background amplitude detectable with PTA, asdiscussed in Kelley et al. (2017b). However, a largersample of mergers, also including cosmological simu-lations would be required to properly characterize theSMBH eccentricity distribution at the onset of the GWdriven inspiral phase.ACKNOWLEDGMENTSThe numerical simulations were performed on facili-ties hosted by the CSC – IT Center for Science, Fin-land. M.M. acknowledges the financial support by theJenny and Antti Wihuri Foundation. M.M., P.H.J., P.P.and A.R. acknowledge the support by the European Re-search Council via ERC Consolidator Grant KETJU(no. 818930).
Software:
Ketju (Rantala et al. 2017, 2018),GADGET-3 (Springel et al. 2005), NumPy (van derWalt et al. 2011), SciPy (Jones et al. 2001–), Matplotlib(Hunter 2007)REFERENCES
Abbott, B. P., et al. 2016, Physical Review Letters, 116,061102, doi: 10.1103/PhysRevLett.116.061102Amaro-Seoane, P., Audley, H., Babak, S., et al. 2017, arXive-prints, arXiv:1702.00786.https://arxiv.org/abs/1702.00786Arun, K. G., Blanchet, L., Iyer, B. R., & Qusailah, M. S. S.2008, PhRvD, 77, 064035,doi: 10.1103/PhysRevD.77.064035Arzoumanian, Z., Baker, P. T., Brazier, A., et al. 2018,ApJ, 859, 47, doi: 10.3847/1538-4357/aabd3bBansal, K., Taylor, G. B., Peck, A. B., Zavala, R. T., &Romani, R. W. 2017, ApJ, 843, 14,doi: 10.3847/1538-4357/aa74e1Begelman, M. C., Blandford, R. D., & Rees, M. J. 1980,Nature, 287, 307, doi: 10.1038/287307a0Berczik, P., Merritt, D., Spurzem, R., & Bischof, H.-P.2006, ApJL, 642, L21, doi: 10.1086/504426 Berentzen, I., Preto, M., Berczik, P., Merritt, D., &Spurzem, R. 2009, ApJ, 695, 455,doi: 10.1088/0004-637X/695/1/455Bernard, L., Blanchet, L., Faye, G., & Marchand , T. 2018,PhRvD, 97, 044037, doi: 10.1103/PhysRevD.97.044037Binney, J., & Tremaine, S. 2008, Galactic Dynamics:Second Edition (Princeton University Press)Blanchet, L. 2014, Living Reviews in Relativity, 17, 2,doi: 10.12942/lrr-2014-2Blanchet, L., Damour, T., Iyer, B. R., Will, C. M., &Wiseman, A. G. 1995, Physical Review Letters, 74, 3515,doi: 10.1103/PhysRevLett.74.3515Bonetti, M., Haardt, F., Sesana, A., & Barausse, E. 2016,MNRAS, 461, 4419, doi: 10.1093/mnras/stw1590Bonetti, M., Sesana, A., Barausse, E., & Haardt, F. 2018,MNRAS, 477, 2599, doi: 10.1093/mnras/sty874Bulirsch, R., & Stoer, J. 1966, Numerische Mathematik, 8,1, doi: 10.1007/BF02165234
Ws from SMBHs in Galactic Simulations2