Gravitational Waves from the Inspiral of Supermassive Black Holes in Galactic-scale Simulations

The simulations analyzed in this paper have been run using the recently developed Ketju code (see Rantala et al. 2017, 2018 for full code details), which is an extension of the tree-SPH simulation code GADGET-3 (Springel et al. 2005). The central idea of the code is the inclusion of a regularized region around every SMBH particle, in which the non-softened gravitational dynamics is computed using the regularized AR-CHAIN (Mikkola & Merritt 2008) integrator while the dynamics of the remaining particles is computed with the GADGET-3 leapfrog using the tree force calculation method.

In practice the code operates by dividing simulation particles into three categories. The SMBH and all stellar particles that lie within a user-defined chain radius $({r}_{\mathrm{chain}})$ are marked as chain subsystem particles. Particles that lie just outside the chain radius but induce a strong tidal perturbation on the chain system are marked as perturber particles. Finally, all the remaining particles that are far from any SMBHs are treated as ordinary GADGET-3 particles with respect to the force calculation. The Ketju code also allows for both multiple simultaneous chain subsystems and several SMBHs in a given single subsystem.

During each global GADGET-3 time step the particles in the chain subsystems are propagated using the AR-CHAIN algorithm (Mikkola & Merritt 2008; Rantala et al. 2017). This algorithm has three main aspects: algorithmic regularization, the use of relative distances to reduce round-off errors by organizing the particles in a chain (e.g., Mikkola & Aarseth 1993), and finally the use of a Gragg–Bulirsch–Stoer (GBS); (Gragg 1965; Bulirsch & Stoer 1966) extrapolation method that yields high numerical accuracy in orbit integrations at a preset user-given error tolerance level $({\eta }_{\mathrm{GBS}})$. In essence, algorithmic regularization works by transforming the equations of motion by introducing a fictitious time variable such that integration by the common leapfrog method yields exact orbits for a Newtonian two-body problem including two-body collisions (e.g., Mikkola & Tanikawa 1999; Preto & Tremaine 1999).

The AR-CHAIN algorithm within the Ketju code features an extension of phase space with the help of an auxiliary velocity variable (Hellström & Mikkola 2010; Pihajoki 2015). This allows for the efficient implementation of the velocity-dependent post-Newtonian corrections in the motion of the SMBH particles (e.g., Will 2006) using an explicit integrator.

Schematically the PN-corrected acceleration can be written as

$$\begin{eqnarray}&&{\boldsymbol{a}}={{\boldsymbol{a}}}_{{\rm{N}}}+{{\boldsymbol{a}}}_{\mathrm{PN}1}+{{\boldsymbol{a}}}_{\mathrm{PN}2}+{{\boldsymbol{a}}}_{\mathrm{PN}3}+{{\boldsymbol{a}}}_{\mathrm{PN}2.5}+{{\boldsymbol{a}}}_{\mathrm{PN}3.5},\end{eqnarray} \tag{ 1 }$$

where the Newtonian acceleration ${{\boldsymbol{a}}}_{{\rm{N}}}$ is computed including the surrounding stellar particles, while the PN terms only include contributions from other SMBHs. The PN-correction terms are labeled so that they are proportional to the corresponding power of the formal PN expansion parameter ${\epsilon }_{\mathrm{PN}}$, i.e.,

$$\begin{eqnarray}&&\left|{{\boldsymbol{a}}}_{i\mathrm{PN}}\right|\propto {\epsilon }_{\mathrm{PN}}^{i}\sim {\left(\displaystyle \frac{v}{c}\right)}^{2i}\sim {\left(\displaystyle \frac{{R}_{s}}{R}\right)}^{i},\end{eqnarray} \tag{ 2 }$$

where v and R are the relative velocity and separation of a pair of SMBHs, R_s = 2GM/c² is the Schwarzschild radius corresponding to the total binary mass M, c is the speed of light, and G the gravitational constant. The PN terms of integer order are conservative and are associated with conserved energy and angular momentum, while the half-integer order terms are dissipative radiation reaction terms caused by the emission of gravitational radiation.

In the Ketju simulations studied in this paper we use PN-correction terms up to order PN3.5 derived for a binary system of arbitrary eccentricity in the modified harmonic gauge as given in Mora & Will (2004). Terms depending on the spin of the black holes as well as the lowest order cross term corrections for a system consisting of more than two particles have also been implemented in the Ketju code, but are not used in this present work. The spin terms are ignored due to the fact that the spin of the SMBHs is largely determined by interactions with gas, which is not included in this study. Similarly, the cross terms as well as other PN corrections for stellar particles are ignored due to the unphysically strong corrections resulting from stellar particles with very large masses $({m}_{\star }\gg {M}_{\odot })$. Thus post-Newtonian corrections are only included in the interactions between SMBHs and only when they are within the same chain region.

We note that the missing higher order PN corrections may potentially lead to significant effects over long periods of time. Corrections up to the PN4-level have been recently derived (e.g., Bernard et al. 2018), but these are highly impractical to implement due to the appearance of the so-called "tail"-terms that depend on the entire history of the binary in addition to its instantaneous state. The lack of these higher order terms also leads to the simulated motion of the SMBHs becoming unreliable when ${\epsilon }_{\mathrm{PN}}\sim 0.1$ (e.g., Csizmadia et al. 2012). This is most apparent when considering the conservation of energy in the system, which we will discuss at the end of Section 5. Due to the increasing unreliability of the PN corrections, we stop simulating the binaries at a separation of R = 6R_s, where the equations of motion are still well behaved. However, as will be discussed in Section 5, the actual limit of reliability for our purposes is closer to R ∼ 10R_s.

Our merger progenitor galaxies consist of spherically symmetric, isotropic multicomponent systems consisting of a stellar bulge, a dark matter halo and a central SMBH. The stellar component follows a Dehnen profile (Dehnen 1993) with γ = 3/2, whereas the dark matter halo is modeled by a Hernquist (γ = 1) profile (Hernquist 1990). The multicomponent initial conditions are 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 for the different mass components are obtained from the corresponding density profile ρ_i and the total gravitational potential Φ_T is then derived using Eddington's formula (Binney & Tremaine 2008).

In this paper we primarily study the orbits and merging of SMBHs in unequal-mass mergers. The primary galaxy is identical to the "γ-1.5-BH-6" initial condition of Rantala et al. (2018, 2019), which was used to model the major merger progenitor of NGC 1600 (Thomas et al. 2016). This model corresponds to a massive gas-free early-type galaxy with a stellar mass ${M}_{\star }=4.15\times {10}^{11}{M}_{\odot }$, with an effective radius ${R}_{{\rm{e}}}\,=7\,\mathrm{kpc}$, and a dark matter halo mass ${M}_{\mathrm{DM}}=7.5\times {10}^{13}{M}_{\odot }$. The dark matter fraction within the effective radius is set to f_DM(R_e) = 0.25 and the central SMBH has a mass ${M}_{\bullet }\,=8.5\times {10}^{9}{M}_{\odot }$. The secondary galaxy is a scaled down version of the primary galaxy, with the masses of all components divided by a factor of 5, i.e., ${M}_{\star }=8.3\times {10}^{10}{M}_{\odot }$, M_DM = 1.5 × 10¹³M_☉, ${M}_{\bullet }=1.7\times {10}^{9}{M}_{\odot }$ and a resulting effective radius R_e = 3.5 kpc (this corresponds to IC-1 in Rantala et al. 2019). For this simulation set the masses of stellar particles are set to ${m}_{\star }={10}^{5}{M}_{\odot }$ and m_DM = 7.5 × 10⁶M_☉ for the dark matter particles, resulting in ${N}_{\star }=4.15\times {10}^{6}$ stellar and N_DM = 1.0 × 10⁷ dark matter particles for the primary galaxy, the number of particles being a factor of five lower for the secondary galaxy.

We first simulate a 5:1 minor merger between the primary galaxy (Simulation A in Table 1). Next we continue the 5:1 merger run with subsequent merger generations until a total of four minor mergers are completed (Simulations B–D in Table 1). The fifth merger generation simulation could not be used for the present study, as the black holes did not merge within a reasonable simulation time due to the formation of a very low density core in this system, with the missing stellar mass effectively giving rise to a final-parsec problem.

Table 1.  Summary of the Studied Merger Runs

Run	${m}_{\bullet 1}$	${m}_{\bullet 2}$	${M}_{\bullet }$	${M}_{\mathrm{DM}}$	${M}_{\star }$	${N}_{\mathrm{DM}}$	${N}_{\star }$	q	a0	e0	t0	${t}_{\mathrm{merger}}$
	$({10}^{9}{M}_{\odot })$	$({10}^{9}{M}_{\odot })$	$({10}^{10}{M}_{\odot })$	$({10}^{13}{M}_{\odot })$	$({10}^{10}{M}_{\odot })$	$(\times {10}^{6})$	$(\times {10}^{6})$		(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.30
B	10.2	1.7	1.19	10.5	58.1	14.0	5.81	6:1	5.45	0.971	1.09	1.22
C	11.9	1.7	1.36	12.0	66.4	16.0	6.64	7:1	4.68	0.961	1.37	1.50
D	13.6	1.7	1.53	13.5	74.7	18.0	7.47	8:1	5.32	0.954	2.10	2.38
X	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_•1, ${m}_{\bullet 2}$), the total binary mass (${M}_{\bullet }$), the total dark matter and stellar masses and particle numbers of the merger remnants (${M}_{\mathrm{DM},\star }$, ${N}_{\mathrm{DM},\star }$), the merger mass ratio q, the SMBH binary orbital parameters a₀ and e₀ at the time t₀, as well as the time of the SMBH binary merger t_merger. The times are measured from the start of the simulation, and the time period considered in this work spans the interval (t₀, t_merger).

Download table as:  ASCIITypeset image

Finally, in addition to the unequal-mass mergers we simulated an equal-mass major merger between two identical galaxies with significantly lower masses of ${M}_{\star }=8.41\times {10}^{10}{M}_{\odot }$, M_DM = 9.71 × 10¹²M_☉, and ${M}_{\bullet }=4\times {10}^{8}{M}_{\odot }$, and an effective radius ${R}_{{\rm{e}}}=4.0\,\mathrm{kpc}$ (this corresponds to the IC in Eisenreich et al. 2017 excluding the hot gas halo). For this run we used ${N}_{\star }\,=8.41\times {10}^{5}$ stellar and ${N}_{\mathrm{DM}}=9.71\times {10}^{6}$ dark matter particles for each galaxy, resulting in particle masses ${m}_{\star }={10}^{5}{M}_{\odot }$ and m_DM = 10⁶ M_☉ for the stellar and dark matter components, respectively, and a final merged SMBH mass ${M}_{\bullet }=8\times {10}^{8}{M}_{\odot }$ (simulation X in Table 1). The motivation for this run was to study the effect of the stellar environment on the SMBH merging process in a setting where the black hole masses are lower by a factor of ∼10 compared to the high SMBH mass A–D simulation set.

All merger orbits are nearly parabolic, with the pericenter distance set to r_p ∼ 0.5R_e of the primary galaxy. After each minor merger (runs A–D), the merger remnant is reoriented so that the satellite galaxies fall in from random directions with respect to the principal axis of the primary galaxy. In Figure 1 we show a sequence of illustrative snapshots of run A, depicting the evolution of the SMBH binary from the dynamical friction phase (left panel), through the stellar scattering stage (middle panel), and into the final gravitational wave driven inspiral phase (right panel).

Zoom In Zoom Out Reset image size

Following Rantala et al. (2018, 2019) we set the GADGET-3 integrator error tolerance to η = 0.002 and the force accuracy to α = 0.005, using the standard cell opening criterion (Springel et al. 2005). The chain radius is set to ${r}_{\mathrm{chain}}=10\,\mathrm{pc}$ and the perturber radius to twice this value, i.e., r_pert = 2 r_chain, for all the simulation runs. The gravitational softening lengths are set to ${\epsilon }_{\star }=3.5\,\mathrm{pc}$ and _DM = 100 pc, respectively. The softening lengths are chosen to fulfill the criterion r_chain > 2.8 ×  (Rantala et al. 2017). In interactions between particles with different softening lengths GADGET-3 uses the larger softening length.

The mergers are initially evolved using Ketju with a chosen GBS tolerance of η_GBS = 10⁻⁶ until the binary black hole semimajor axis is around a ∼ 5000R_s. This corresponds to about ${a}_{0}\sim 5\,\mathrm{pc}$ for runs A–D and ${a}_{0}\sim 0.5\,\mathrm{pc}$ for the lower mass run X, with all SMBH binaries having high eccentricities, in excess of e₀ ≳ 0.9. These separations are typically reached within t₀ ∼ 1–2 Gyr after the start of the simulation (see Table 1).

At this point in time, we restart our Ketju runs with a stricter GBS tolerance value of η_GBS = 10⁻⁹. The increased accuracy is required for calculating the resulting GW spectrum and resolving in detail the post-Newtonian orbital motions of the individual SMBHs all the way down to their final coalescence, which typically occurs Δt_merger ∼ 100–300 Myr later (see Table 1). The simulation of the binary in these runs ends when the estimated merger timescale of the binary is smaller than the time step used for the global simulation. At this point the separation of the binary is typically 50–70 Schwarzschild radii.

We use the SMBH positions and velocities computed in the full Ketju runs that include the influence of the stellar background for our GW calculations until the separation of the SMBHs is ≲100R_s. This distance corresponds to $\sim 0.1\,\mathrm{pc}$ for runs A–D, and to $\sim 8\times {10}^{-3}\,\mathrm{pc}$ for run X. The final part of the merger is then run separately without the surrounding stellar particles at a much higher time resolution to facilitate direct waveform computations as described in Section 4. The transition separation of ∼100R_s was chosen as it provides a good balance between computational cost and the accuracy of these calculations. While in principle stellar particles could also be included in this phase, we find that ignoring the environment in this very final stage does not cause significant changes in the SMBH evolution, as only very few stellar particles are close enough to influence the SMBH binary at our mass resolution. This has also been checked by comparing the evolution of the isolated binaries to the results from the full Ketju run down to separations of ∼70R_s, where the lower time resolution Ketju data are still available.

For an analysis of the evolution of the SMBH binary orbit, it is necessary to have a parameterization of the orbit that changes slowly and can be recovered from the instantaneous position and velocity of the binary. Such a parameterization also allows for computing the gravitational wave emission of the system from only sparsely stored data points, significantly reducing the computational effort.

In Newtonian gravity the orbit of a binary system can be described by the Keplerian orbital elements including the semimajor axis a and the eccentricity e, which are constant or evolve only slowly due to external perturbations. However, due to the inclusion of PN corrections to the motion of the SMBH binary, the standard Keplerian elements are no longer even approximately constant over an orbit, but rather oscillate quite strongly especially near the pericenter of an eccentric orbit, as is illustrated in Figure 2.

Zoom In Zoom Out Reset image size

As a result of this, the Keplerian parameters computed from sparsely stored data suffer from essentially random fluctuations due to the varying phase of the orbit at the data points, which is problematic both for analyzing the evolution of the orbit as well as for GW computations. For small binary separations, the Keplerian orbital parameters are also no longer related to the geometry of the orbit, as is apparent from the increase of the Keplerian eccentricity for a nearly circular binary in the right panel of Figure 2.

The problems associated with the Keplerian parameters can be mostly remedied by using a quasi-Keplerian parameterization that accounts for the effects of the PN corrections. Several such parameterizations accurate to different orders have been developed, and in this work we utilize the PN3 accurate quasi-Keplerian orbital parameterization of Memmesheimer et al. (2004). This parameterization gives an approximate solution to the conservative part of the PN3 accurate equations of motion of a nonspinning binary in the form

$$\begin{eqnarray}&&R=a(1-{e}_{R}\cos u),\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}\begin{array}{rcl}n(t-{t}_{0}) & = & u-{e}_{t}\sin u+\left(\displaystyle \frac{{g}_{4t}}{{c}^{4}}+\displaystyle \frac{{g}_{6t}}{{c}^{6}}\right)(v-u)\\ & & +\left(\displaystyle \frac{{f}_{4t}}{{c}^{4}}+\displaystyle \frac{{f}_{6t}}{{c}^{6}}\right)\sin v\\ & & +\displaystyle \frac{{i}_{6t}}{{c}^{6}}\sin 2v\\ & & +\displaystyle \frac{{h}_{6t}}{{c}^{6}}\sin 3v,\end{array}\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{2\pi }{{\rm{\Phi }}}(\phi -{\phi }_{0}) & = & v+\left(\displaystyle \frac{{f}_{4\phi }}{{c}^{4}}+\displaystyle \frac{{f}_{6\phi }}{{c}^{6}}\right)\sin 2v\\ & & +\left(\displaystyle \frac{{g}_{4\phi }}{{c}^{4}}+\displaystyle \frac{{g}_{6\phi }}{{c}^{6}}\right)\sin 3v+\displaystyle \frac{{i}_{6\phi }}{{c}^{6}}\sin 4v\\ & & +\displaystyle \frac{{h}_{6\phi }}{{c}^{6}}\sin 5v,\end{array}\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&\tan \displaystyle \frac{v}{2}=\sqrt{\displaystyle \frac{1+{e}_{\phi }}{1-{e}_{\phi }}}\tan \displaystyle \frac{u}{2}.\end{eqnarray} \tag{ 6 }$$

The coordinates R and ϕ are the relative separation and angle in the orbital plane, with ϕ₀ fixing the direction of pericenter at time t = t₀. As in the Keplerian case, the orbit is an ellipse with semimajor axis a, but now the frequency of radial oscillations f_r = n/2π is not the same as the frequency of angular motion. Rather, during a single radial period the direction angle ϕ traverses an angle Φ > 2π leading to the precession of the orbit. The angles u and v are generalizations of the eccentric and true anomaly, although their interpretation in terms of the geometry of the orbit is not as straightforward as in the Keplerian case due to the appearance of the additional series expansion factors g, f, i, and h in the generalization of Kepler's Equation (4) and in Equation (5). In addition the single eccentricity of the Keplerian parameterization is replaced by three different eccentricities e_R, e_t, and e_ϕ. In the limit where the post-Newtonian corrections are negligible this parameterization reduces to the standard Keplerian parameterization.

The orbital parameters appearing in the above parameterization can be computed from the binary separation vector ${\boldsymbol{R}}$ and relative velocity ${\boldsymbol{V}}$ in the modified harmonic gauge used by Ketju via PN3 accurate approximations of the conserved energy E and angular momentum J in the form

$$\begin{eqnarray}&&E=\mu \left({E}_{0}+{c}^{-2}{E}_{1}+{c}^{-4}{E}_{2}+{c}^{-6}{E}_{3}\right)\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&J=\mu \left|{\boldsymbol{R}}\times {\boldsymbol{V}}\right|\left({J}_{0}+{c}^{-2}{J}_{1}+{c}^{-4}{J}_{2}+{c}^{-6}{J}_{3}\right),\end{eqnarray} \tag{ 8 }$$

where μ = m₁m₂/M is the reduced mass and M = m₁ + m₂ is the total mass of the binary. The expressions for the E_i, J_i in terms of ${\boldsymbol{R}}$ and ${\boldsymbol{V}}$ are quite lengthy, and can be found in Memmesheimer et al. (2004). These approximations of the conserved energy and angular momentum are not constant along an orbit even when the radiation reaction terms are not included. Instead they oscillate slightly, and this oscillatory behavior carries over to the orbital parameters, which are themselves given as PN3 accurate expressions in terms of the approximated E and J. The inclusion of the radiation reaction terms also introduces additional periodic oscillations (Königsdörffer & Gopakumar 2006), but all of the oscillations are typically much smaller than the resulting secular evolution, and thus of no significant consequence.

The appearance of three different eccentricities causes some ambiguity in their interpretation, as they do not correspond to a single orbital ellipse as in the Keplerian case. The radial eccentricity e_R appears to be most directly related to the dimensions of the orbit, but unfortunately the PN3 formulae for it and e_ϕ fail for nearly circular relativistic orbits, resulting in ${e}_{R}^{2}\lt 0$ and e_ϕ² < 0. The time eccentricity e_t remains real, but it shows a small increase in this regime which does not correspond to an actual increased eccentricity of the orbit. A measure of eccentricity that describes the actual geometry of the orbit even in the strongly relativistic regime can be defined as (e.g., Csizmadia et al. 2012)

$$\begin{eqnarray}&&{e}_{g}=\displaystyle \frac{{R}_{\max }-{R}_{\min }}{{R}_{\max }+{R}_{\min }},\end{eqnarray} \tag{ 9 }$$

where R_min, R_max are the peri- and apocenter separations of the binary. This definition is impractical since it cannot be computed from the instantaneous relative position and velocity of the binary, but it serves as a useful comparison.

The right panel of Figure 2 shows the near-merger behavior of the different eccentricities for an isolated binary with an eccentricity that is comparable to the ones seen in the main Ketju runs at this separation. The time eccentricity e_t can be seen to be quite close to the geometric eccentricity e_g for the most part, while both e_R and e_ϕ go to zero too early, which is caused by them becoming imaginary at that point. In regimes where the PN expansion parameter ${\epsilon }_{\mathrm{PN}}$ is small the differences between the different eccentricities are negligible, so in the following we elect to use $e\equiv {e}_{t}$ as the single parameter describing the eccentricity of the orbit. However, one should keep in mind that for the final stages of the binary inspiral the eccentricity defined this way does not exactly describe the shape of the orbit. The other orbital elements likely suffer from similar deterioration in the very final stages of the inspiral, but for the vast majority of the binary orbital evolution they work reliably.

The stellar environment affects the evolution of the binary orbit in addition to the evolution caused by GW emission. In the final stages of the binary evolution considered in this work these effects are however fairly small, and to analyze them it is necessary to separate them from the main GW driven orbital evolution. The time derivative of a given quantity X, which can be, e.g., the energy or the eccentricity, can be split into a part corresponding to the effect from the PN forces including the GW reaction and a part corresponding to additional effects caused by interactions with the stellar environment:

$$\begin{eqnarray}&&\displaystyle \frac{{dX}}{{dt}}={\left.\displaystyle \frac{{dX}}{{dt}}\right|}_{\mathrm{PN}}+{\left.\displaystyle \frac{{dX}}{{dt}}\right|}_{\mathrm{env}}.\end{eqnarray} \tag{ 10 }$$

From this we find the cumulative effect of the environment as

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{\mathrm{env}}X(t)={\int }_{{t}_{0}}^{t}{\left.\displaystyle \frac{{dX}}{{dt}}\right|}_{\mathrm{env}}{dt}=X(t)-{\int }_{{t}_{0}}^{t}{\left.\displaystyle \frac{{dX}}{{dt}}\right|}_{\mathrm{PN}}{dt}.\end{eqnarray} \tag{ 11 }$$

This quantity is more robust than the instantaneous derivative ${\left.{dX}/{dt}\right|}_{\mathrm{env}}$, as numerically differentiating the X(t) derived from a Ketju run to compute the instantaneous derivative amplifies any numerical noise present in the data, which can easily hide the small effects of the stellar environment.

To compute the effect of the PN forces at a given time t', i.e., ${\left.{dX}/{dt}\right|}_{\mathrm{PN}}(t^{\prime} )$, we perform a short integration of an isolated binary with initial conditions taken from the Ketju computed SMBH positions and velocities. The integration is performed for about 30 orbits of the binary, and the derivative at the middle of the integration interval is estimated by fitting a line to the densely sampled evolution of X over this period. The initial conditions are chosen so that the middle of the integration interval coincides with the desired time t'.

We evaluate the derivatives at about 50,000 points spaced evenly in the semimajor axis a over the considered section of the Ketju results. The cumulative effect ${\int }_{{t}_{0}}^{t}{\left.{dX}/{dt}\right|}_{\mathrm{PN}}{dt}$ is then found by integrating the instantaneous derivatives numerically. With these choices this method of computing the instantaneous derivatives caused by the PN terms appears to work well until the binary separation is under ∼100R_s, after which it fails as the evolution of the parameters starts to become nonlinear over the integration period.

For comparison against the SMBH binary evolution computed with Ketju we consider two semi-analytic models used to describe binaries in simulations that cannot resolve them.

The first model is the widely used Peters (1964) result, which considers the secular evolution of the Keplerian orbital parameters of an isolated binary due to the leading radiation reaction term at PN2.5 level. The binary is assumed to be otherwise Keplerian, and the orbital parameters a and e are assumed to change only slowly. The evolution of the semimajor axis a is then given by

$$\begin{eqnarray}&&\displaystyle \frac{{da}}{{dt}}=-\displaystyle \frac{64}{5}\displaystyle \frac{{G}^{3}{m}_{1}{m}_{2}M}{{c}^{5}{a}^{3}{\left(1-{e}^{2}\right)}^{7/2}}\left(1+\displaystyle \frac{73}{24}{e}^{2}+\displaystyle \frac{37}{96}{e}^{4}\right)\end{eqnarray} \tag{ 12 }$$

and the eccentricity e evolution by

$$\begin{eqnarray}&&\displaystyle \frac{{de}}{{dt}}=-\displaystyle \frac{304}{15}\displaystyle \frac{{G}^{3}{m}_{1}{m}_{2}M}{{c}^{5}{a}^{4}{\left(1-{e}^{2}\right)}^{5/2}}e\left(1+\displaystyle \frac{121}{304}{e}^{2}\right).\end{eqnarray} \tag{ 13 }$$

The main effect of the stellar environment on the evolution of the SMBH binary is through the scattering of individual stars during close encounters. These scattering events remove energy and angular momentum from the binary, decreasing the semimajor axis a and increasing the eccentricity e. The second semi-analytic model considered here is a commonly used (e.g., Sesana 2010; Kelley et al. 2017b) model describing this binary hardening process (Hills & Fullerton 1980; Quinlan 1996), which gives the evolution of the semimajor axis and eccentricity as

$$\begin{eqnarray}&&\displaystyle \frac{d}{{dt}}\left(\displaystyle \frac{1}{a}\right)=\displaystyle \frac{G\rho }{\sigma }H\end{eqnarray} \tag{ 14 }$$

and

$$\begin{eqnarray}&&\displaystyle \frac{{de}}{{dt}}=-{{Ka}}^{-1}\displaystyle \frac{{da}}{{dt}}=K\displaystyle \frac{G\rho }{\sigma }{Ha},\end{eqnarray} \tag{ 15 }$$

where ρ is the stellar density and σ the velocity dispersion. Equation (14) can also be written in terms of the energy of the binary as

$$\begin{eqnarray}&&\displaystyle \frac{{dE}}{{dt}}=-\displaystyle \frac{{GM}\mu }{2}\displaystyle \frac{G\rho }{\sigma }H,\end{eqnarray} \tag{ 16 }$$

where M = m₁ + m₂ is the total and μ = m₁m₂/M the reduced mass of the binary. This form is more suited for our applications as the post-Newtonian corrections modify the relation between the energy E and semimajor axis a compared to the nonrelativistic Keplerian case.

The constants H and K are usually determined by three-body scattering experiments (Quinlan 1996; Sesana et al. 2006), but Sesana & Khan (2015) find that similar values for the constants can also be recovered in N-body simulations when the stellar density and velocity dispersion are evaluated at the influence radius of the SMBH binary. Traditionally, the influence radius ${r}_{\inf }$ is defined by the condition that the enclosed stellar mass is twice the binary mass, i.e., that ${M}_{* }(r\lt {r}_{\inf })=2{M}_{\bullet }$.

We focus our gravitational wave computations on quantities relevant for observations with PTAs, which can detect the low frequency GWs emitted by the massive systems studied in this work. The main observable quantity describing the stochastic GWB is the characteristic strain h_c(f), which can be computed using either semi-analytic or Monte Carlo (MC) methods when the distribution of sources and their GW emission is known.

In the semi-analytic method, the characteristic strain h_c(f) of the GWB at frequency f is given by (Phinney 2001; Enoki & Nagashima 2007)

$$\begin{eqnarray}&&{h}_{c}^{2}({f}_{o})=\displaystyle \frac{4G}{\pi {c}^{2}f}\int {dz}\,{n}_{c}{\left.\displaystyle \frac{{{dE}}_{\mathrm{GW}}}{{df}}\right|}_{f=(1+z){f}_{o}},\end{eqnarray} \tag{ 17 }$$

where n_c is the comoving number density of sources and the source distribution is assumed to be spatially uniform. As the source density n_c is determined by the large-scale evolution of the universe, here we compute only ${{dE}}_{\mathrm{GW}}/{df}$, the GW spectral energy density (SED) of a source integrated over its lifetime, which is thus affected by any processes speeding up or delaying the merger of the SMBH binaries. In the MC methods of computing the characteristic strain (e.g., Sesana et al. 2008), the characteristic strain is instead computed by simply summing the instantaneous contributions from each source. Thus all the environmental effects are included in the source distribution which we do not consider here.

For computing the instantaneous GW flux and the SED of an SMBH binary we use two different methods depending on the regime: in the weakly relativistic regime we save computational effort by applying analytically derived formulae to compute the GW spectrum from the orbital elements, while in the regime where the relativistic effects become significant we compute the GW signal directly from the motion of the BHs. This split between methods allows us to take advantage of the full PN3.5 accurate dynamics used in Ketju while saving computational effort where possible without sacrificing numerical accuracy.

When the post-Newtonian effects are small, it is possible to compute the contribution to the SED to good accuracy using analytic formulae derived assuming a slowly varying Keplerian orbit.

A classic and widely used result derived by Peters & Mathews (1963) is the power emitted in GWs by a binary in the lowest order quadrupole approximation. For an eccentric orbit with eccentricity e and orbital frequency f_p, the GWs are emitted in all harmonics of the orbital frequency, with the nth harmonic having the power

$$\begin{eqnarray}&&{P}_{n}=\displaystyle \frac{32{G}^{5/3}{m}_{1}^{2}{m}_{2}^{2}}{5{c}^{5}{a}^{2}}{\left(2\pi {f}_{p}\right)}^{2}g(n,e).\end{eqnarray} \tag{ 18 }$$

The function g(n, e) is given by

$$\begin{eqnarray}\begin{array}{rcl}g(n,e) & = & \displaystyle \frac{{n}^{4}}{32}\left\{\Space{0ex}{0.25ex}{0ex}\left[\Space{0ex}{3.25ex}{0ex}{J}_{n-2}({ne})-2{{eJ}}_{n-1}({ne})\right.\right.\\ & & +{\left.\displaystyle \frac{2}{n}{J}_{n}({ne})+2{{eJ}}_{n+1}({ne})-{J}_{n+2}({ne})\right]}^{2}\\ & & +(1-{e}^{2}){\left[{J}_{n-2}({ne})-2{J}_{n}({ne})+{J}_{n+2}({ne})\right]}^{2}\\ & & +\left.\displaystyle \frac{4}{3{n}^{2}}{J}_{n}{\left({ne}\right)}^{2}\right\},\end{array}\end{eqnarray} \tag{ 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 may simply use the PN quasi-Keplerian parameters in place of the Keplerian parameters, thus taking f_p=f_r.

From the instantaneous GW power and the known evolution of the orbital elements the SED can be computed as (Enoki & Nagashima 2007)

$$\begin{eqnarray}&&\displaystyle \frac{{{dE}}_{\mathrm{GW}}}{{df}}=\displaystyle \sum _{n=1}^{\infty }{\left.\displaystyle \frac{{P}_{n}({f}_{p},e){\tau }_{\mathrm{GW}}({f}_{p},e)}{{{nf}}_{p}}\right|}_{{f}_{p}=f/n},\end{eqnarray} \tag{ 20 }$$

where

$$\begin{eqnarray}&&{\tau }_{\mathrm{GW}}={f}_{p}\displaystyle \frac{{dt}}{{{df}}_{p}}=\displaystyle \frac{{dt}}{d(\mathrm{ln}{f}_{p})}\end{eqnarray} \tag{ 21 }$$

is the gravitational wave timescale. Implementing this formula for computations from numerical simulation outputs is mostly straightforward: for each frequency f we truncate the sum at some sufficiently high number of harmonics and interpolate the eccentricity e and gravitational wave timescale τ_GW to the required orbital frequencies f_p. Here e and τ_GW are treated as functions of the orbital frequency f_p instead of time t, which is possible as generally f_p(t) is monotonic in the relevant phase of the binary evolution. The number of harmonics used for the computation is taken to be 750 in this work, which gives converged results for the eccentricities encountered.

The computation of τ_GW from the numerical data requires some care, as naively numerically differentiating the time-frequency data will amplify any noise caused by, for example, interactions with the environment to unusable levels. To avoid this, we instead compute the timescale by integrating the time spent in logarithmic bins of orbital frequency. If two adjacent data points belong to the same bin the whole interval is assigned to that bin, while for points in different bins the time interval is divided over the spanned frequency interval. We find that using about 50 bins per decade in frequency gives results that match well with the differentiation of smooth data.

The error in applying these leading order or PN0 results is expected to be of order PN1, and in Section 4.4 we find that the error between this method and the more accurate method described in the next subsection is of the order of a few percent when the PN expansion parameter ${\epsilon }_{\mathrm{PN}}\approx 0.01$, in line with expectations.

Similar methods are commonly extended also beyond the regime where PN effects are minor all the way down to the final orbits before the BHs merge (e.g., Berentzen et al. 2009; Kelley et al. 2017b). This may work well if the orbital evolution is consistently performed in the same approximation, but coupled with PN3.5 accurate orbits this leads to clearly erroneous results. There are also analytic results for PN-accurate spectra available in the literature (e.g., Tessmer & Schäfer 2010; Klein et al. 2018), however these are quite cumbersome to apply, and we find the methods of the next subsection to be more convenient.

When the effects of the PN corrections on the now non-Keplerian binary orbit become significant, the methods of the previous subsection are no longer reliable. In this regime we instead directly compute the emitted waveform, and use a Fourier transform to compute the spectrum.

In principle this method could be used over the entire binary evolution, but for most part the semi-analytic method is sufficiently accurate. Further, if the binary eccentricity is high, the time step required to resolve the highest harmonics of the orbital frequency becomes very small resulting in considerable computational cost. For these reasons we limit this method to the short period before the BH binary merger where it is necessary for accurate results. In this regime the effect of the environment is also negligible, allowing us to easily compute a densely outputted evolution of the SMBH binary in isolation using the Ketju integrator. Inclusion of the environment at this stage would also be possible, but this would require suitable schemes to resample the stellar population to lower mass particles in order to avoid artifacts in the spectrum caused by spurious unphysically strong interactions from massive stellar particles.

In computing the GW waveform emitted by the SMBH binary we follow Poisson & Will (2014). The leading contribution to the strain tensor h_ij is the quadrupole term

$$\begin{eqnarray}&&{h}_{{ij}}=4\displaystyle \frac{{{Gm}}_{1}{m}_{2}}{{c}^{4}{dM}}\left({V}_{i}{V}_{j}-{GM}\displaystyle \frac{{R}_{i}{R}_{j}}{{R}^{3}}\right),\end{eqnarray} \tag{ 22 }$$

where the relative velocity V_i and separation R_i are computed directly from the numerical integration of SMBH trajectories. The observer distance d only scales the result, and can in practice be ignored in these computations. The radiation reaction force corresponding to this part appears as the ${{\boldsymbol{a}}}_{\mathrm{PN}2.5}$ term in the equations of motion, and to maintain consistency with the PN3.5 accurate equations of motion we also include PN corrections to the waveform up to PN1 order beyond the leading PN0 contribution of Equation (22). The formulae for these corrections for a binary system are lengthy, and are not reproduced here but can be found in Poisson & Will (2014). Note that by convention the waveform PN orders are offset from the corresponding effects in the equations of motion by 2.5 orders.

The strain tensor h obtained as a function of time is then projected into the physically observed polarizations h₊ and h_×, which are related to emitted GW power as

$$\begin{eqnarray}&&\displaystyle \frac{{{dE}}_{\mathrm{GW}}}{{dt}}=\displaystyle \frac{{d}^{2}{c}^{3}}{16\pi G}\int \left({\left({\partial }_{\tau }{h}_{+}\right)}^{2}+{\left({\partial }_{\tau }{h}_{\times }\right)}^{2}\right)d{\rm{\Omega }},\end{eqnarray} \tag{ 23 }$$

where τ = t − d/c is the retarded time, which can here be replaced with the coordinate time t due to the arbitrary observer distance. In terms of the Fourier coefficients ${\tilde{h}}_{+,\times }(f)$ the time-dependent signals can be written as

$$\begin{eqnarray}&&{h}_{+,\times }(t)=\displaystyle \sum _{k}\sqrt{2}\left|{\tilde{h}}_{+,\times }({f}_{k})\right|\sin (2\pi {f}_{k}t+{\phi }_{k}),\end{eqnarray} \tag{ 24 }$$

so the total energy emitted over a time period of length T can be written as

$$\begin{eqnarray}\begin{array}{rcl}{\rm{\Delta }}{E}_{\mathrm{GW}} & = & {\int }_{{t}_{0}}^{{t}_{0}+T}{dt}\displaystyle \frac{{{dE}}_{\mathrm{GW}}}{{dt}}\\ & = & \displaystyle \frac{{d}^{2}{c}^{3}T}{16\pi G}\displaystyle \sum _{\displaystyle \genfrac{}{}{0em}{}{k}{p=+,\times }}{\left(2\pi {f}_{k}\right)}^{2}\int {\left|{\tilde{h}}_{p}({f}_{k})\right|}^{2}d{\rm{\Omega }}.\end{array}\end{eqnarray} \tag{ 25 }$$

As this can also be written as

$$\begin{eqnarray}&&{\rm{\Delta }}{E}_{\mathrm{GW}}=\displaystyle \sum _{k}\displaystyle \frac{{\rm{\Delta }}{E}_{\mathrm{GW}}({f}_{k})}{{\rm{\Delta }}f}{\rm{\Delta }}f,\end{eqnarray} \tag{ 26 }$$

the SED at the frequency f_k is given by the corresponding term in the sum

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{{dE}}_{\mathrm{GW}}({f}_{k})}{{df}} & \approx & \displaystyle \frac{{\rm{\Delta }}{E}_{\mathrm{GW}}({f}_{k})}{{\rm{\Delta }}f}\\ & = & \displaystyle \frac{{c}^{3}{\left(2\pi {f}_{k}{Td}\right)}^{2}}{16G}\displaystyle \sum _{p=+,\times }\int {\left|{\tilde{h}}_{p}({f}_{k})\right|}^{2}d{\rm{\Omega }},\end{array}\end{eqnarray} \tag{ 27 }$$

where Δf = f_k+1−f_k = 1/T is the frequency resolution of a discrete Fourier transform of a signal with evenly spaced samples.

To numerically compute these quantities we use the fast Fourier transform (FFT) algorithm, which allows for computing the coefficients ${\tilde{h}}_{+,\times }$ efficiently from evenly spaced samples of the waveform. The integral over the solid angle Ω is conveniently performed numerically, although at least when ignoring the direction-dependent PN corrections to the waveform it can also be treated analytically. The changing frequency of the signal leads to rapid oscillations in the FFT results, which we suppress by averaging them over a small window. Additionally, the abrupt changes at the ends of the signal lead to some artifacts. These could be reduced by windowing, but this would also reduce the accuracy of the computed energy of the GW signal. From the Fourier transformed signals we may also evaluate the instantaneous power from Equation (23) using the relation between derivatives and Fourier transforms. Evaluating the time derivatives via Fourier transforms is also more accurate for periodic signals than using finite difference techniques.

In order to confirm that the semi-analytic method gives good results in the mildly relativistic regime and to find the limits of its applicability we compare it to the direct waveform calculation over short sections of the evolution of an isolated binary. This allows also for the evaluation of the effect of including the additional PN corrections to the waveform instead of simply using the quadrupole formula (22).

For the mildly relativistic regime we compute the orbit of an isolated SMBH binary with a mass ratio q = 5:1 and initial semimajor axis a = 250R_s and eccentricity e=0.2 over $2\times {10}^{4}$ orbital periods. This fairly large number of orbits allows the orbital frequency to evolve enough for the GW spectrum to be resolved instead of effectively consisting of δ-spikes which are problematic numerically. The post-Newtonian expansion parameter is approximately ${\epsilon }_{\mathrm{PN}}\approx 1/200$ in this case, so the effect of the additional PN effects not taken into account in the semi-analytic model can be expected to be of the order of half a percent.

The left panel of Figure 3 shows the resulting spectrum around the strongest n = 2 harmonic. Apart from the small shift in frequency, which is of the order ${\epsilon }_{\mathrm{PN}}$, the semi-analytic result matches well with the direct waveform computation using the full PN1 accurate waveform. This agreement with the PN1 waveform may seem surprising, as the semi-analytic method only includes the PN0 quadrupole radiation, but this can be simply explained based on energy conservation. In this regime the binary is still close to Newtonian and the main effect of the PN1 waveform corrections is to maintain energy balance with the PN3.5 radiation reaction terms. As a result, the semi-analytic method, which is also constructed around energy conservation of a Newtonian binary, gives matching results. The frequency shift is caused by the precession of the binary, which leads to the GW frequencies differing from integer multiples of the radial frequency used for the calculation. While this shift could in principle be corrected for, we find that it is insignificant when the spectrum is computed over the entire binary evolution.

Zoom In Zoom Out Reset image size

To evaluate the methods in the highly relativistic regime we compute the evolution of the binary from an initial semimajor axis a = 40R_s down to a = 6R_s, where the PN equations of motion become unreliable. The initial eccentricity is set to e ≈ 0.1, which is comparable to the eccentricities found in our Ketju runs at this binary separation. The resulting spectra in the right panel of Figure 3 clearly show the increasing significance of consistently computing the waveform from the actual binary motion as well as including the PN corrections to the waveform, with the differences between the methods growing to tens of percents toward the end of the simulation.

The error of the semi-analytic method compared to the full PN waveform is about 5% at the point where all the major harmonics are visible in the spectrum, but somewhat lower when compared to the PN0 waveform. In this regime it is also of interest to look at the significance of the higher order radiation reaction terms at PN3.5 level. The right panel of Figure 3 includes also the quadrupole spectrum from an otherwise identical binary evolved using PN corrections only up to PN3 level. The amplitude of the spectrum is clearly different from the full PN3.5 result, showing the importance of including the higher order radiation reaction term. Interestingly, the semi-analytic method yields essentially exact results compared to the PN0 waveform in this case. This indicates that the differences between the semi-analytic method and the direct waveform methods are at least partly due to the PN3.5 radiation reaction term causing significant changes in the motion of the binary.

We begin by studying the evolution of the orbital elements of the SMBH binaries in the Ketju runs and compare them to the parameters predicted by semi-analytic models describing both isolated binaries and binaries interacting with a population of stars. The semi-analytic comparison models are started from the initial conditions at the start of the analyzed section, i.e., using the values of a₀ and e₀ listed in Table 1. This is similar to what would be done in a standard softened simulation when the SMBH binary separation has decreased to the order of a gravitational softening length.

Isolated binaries are modeled using the model described in Section 3.3, hereafter the Peters model, which describes a Keplerian binary perturbed by the leading gravitational radiation reaction term.

The interaction with stars is included by combining the Peters model with the scattering model of Section 3.4, to produce what we call the Peters-Quinlan model. This model requires as inputs the stellar density and velocity dispersion, which are computed from the Ketju run results in a sphere within the influence radius of the binary. For these calculations we use the values of H and K parameters determined by the fitting formulae of Sesana et al. (2006), which allows us to see how much the behavior of the Ketju runs differs for typically used literature values. The values of H are computed using the fit for a circular binary, as Sesana et al. (2006) do not tabulate the coefficients for the H fit for eccentric binaries. Their results however indicate that at higher eccentricities the value of H should tend to increase slightly.

The time evolution of the semimajor axis of the orbits in the Ketju runs and the semi-analytic comparison models is shown in Figure 4. The evolution of the binary is mainly driven by GW emission at this stage, but the environment still has a nonnegligible effect. This can be seen by studying the isolated Peters models, in which the SMBHs merge significantly later than in the Ketju runs. On the other hand, the Peters-Quinlan models merge significantly earlier than the Ketju runs, which is caused by the fact that the literature values of the H and K parameters from Sesana et al. (2006) significantly overestimate the effects of stellar scattering for our Ketju runs with depleted stellar cores (runs A–D). This is most obvious in the case of run A, while in the case of run X the match between the Ketju run and the Peters-Quinlan model is reasonably good.

Zoom In Zoom Out Reset image size

The eccentricity evolution of the runs is shown in Figure 5 as a function of the orbital frequency ${f}_{r}\approx \sqrt{{GM}/{a}^{3}}/(2\pi )$, as this allows a comparison between the models with different merger times and also relates the eccentricity to the frequency of emitted GWs. Runs A–D all show fairly similar behavior with high eccentricities in the range e ≈ 0.96–0.98 at the initial semimajor axis $a\sim 5\,\mathrm{pc}$ corresponding to an orbital frequency ${f}_{r}\sim {10}^{-4}\,{\mathrm{yr}}^{-1}$ where we begin considering these runs. Run X has a slightly lower eccentricity at its initial $a\approx 0.5\,\mathrm{pc}$, but it too is over 0.9.

Zoom In Zoom Out Reset image size

The differences between the two semi-analytic models reveal an additional reason for the different merger times, in addition to the direct decay of the semimajor axis due to stellar scattering: the eccentricities of the Peters model are considerably lower and those of the Peters-Quinlan model higher than in the Ketju runs. As can be seen from Equation (12), the dependence of the evolution on the eccentricity is highly nonlinear, and especially at high eccentricities even small differences lead to large effects in the merger time. Notably the difference in the evolution of the eccentricity is smallest in the case of the Peters-Quinlan model run X and largest in the case of run A, which show also correspondingly the smallest and largest difference in the SMBH merger times, respectively.

To quantify the strength of the environmental effects, we use the method of Section 3.2 to subtract the dominant effects of the PN corrections and recover the residual effect caused directly by interactions with the stellar environment. The results for the energy and eccentricity 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 the correct sign, but the effect is fairly weak.

Zoom In Zoom Out Reset image size

We perform linear fits to the data in Figure 6 to determine the values of the H and K parameters of the scattering model that would correctly describe the binary evolution when the stellar densities and velocity dispersions are derived from the Ketju runs. Note that while Equation (16) implies that the change in energy is linear in time, Equation (15) implies that the change is eccentricity is linear in the integrated semimajor axis:

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{\mathrm{env}}e(t)={KH}\displaystyle \frac{G\rho }{\sigma }{\int }_{{t}_{0}}^{t}{adt}.\end{eqnarray} \tag{ 28 }$$

The results are given in Table 2, which also shows the H_S and K_S values used for the Peters-Quinlan comparison models, derived from Sesana et al. (2006). In the case of run X the fitted value of H is close to the value determined from the Sesana et al. (2006) fitting formulae, being about 25% lower. This is comparable to the results 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 Sesana et al. (2006).

For runs A–D the agreement between the fitted H parameter and the one used for the Peters-Quinlan model is much poorer, with up to a factor of eight difference in the case of run A. This is not entirely unexpected as these runs describe massive early-type galaxies with very large low density stellar cores. Thus, the derived small values of H are caused by the fairly empty loss cones in these runs, as most of the stars on radial trajectories have already been ejected at this stage by the strong core 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. However, the differences are not quite as large as in the case of the H parameter. On the other hand, the fitted value of K for run X is almost twice as large as the literature value used for comparison. This is likely partly due to effects of eccentricity on the parameters, as the comparison value was computed for a circular binary. For run X the eccentricity evolution between the Peters-Quinlan model and the Ketju run agrees very well, despite slight differences in the scattering parameters (see Table 2). This is made possible since there is a certain degree of degeneracy between the H and K parameters in the Peters-Quinlan model, allowing very similar evolutions to be produced for a range of parameter values.

Table 2.  Fitted and Values from the Literature of the Hardening Parameters

Run	${H}_{{\rm{K}}{\rm{ETJU}}}$	${K}_{{\rm{K}}{\rm{ETJU}}}$	HS	KS	Gρ/σ
					(${10}^{-11}\,{\mathrm{pc}}^{-1}\,{\mathrm{yr}}^{-1}$)
A	2.18	0.032	16.2	0.056	5.49
B	4.38	0.049	16.3	0.059	2.92
C	5.98	0.052	16.4	0.068	1.68
D	8.44	0.027	16.6	0.068	1.03
X	11.5	0.098	14.4	0.054	24.4

Note. The values of the hardening parameters ${H}_{{\rm{K}}{\rm{ETJU}}}$ and ${K}_{{\rm{K}}{\rm{ETJU}}}$ of Section 3.4 have been computed by linear fits to the Ketju data with the effects of the PN3.5 level corrections subtracted (see Section 3.2). The values used for the Peters-Quinlan comparison models (H_S, K_S) computed from the fitting formulae of Sesana et al. (2006) are also listed, as well as the value of Gρ/σ measured from the simulations at the influence radii of the SMBH binaries.

Download table as:  ASCIITypeset image

We also computed comparison Peters-Quinlan models using the parameter values fitted from the Ketju runs. For all runs the results are accurate to a similar level as the run X comparison model using literature values discussed above, i.e., at the level of a few percent. As the behavior is so similar, we opt to not display the results of these calculations.

Next we focus on the GW signal emitted by the SMBH binaries, and the differences between the signals computed from the Ketju runs and the semi-analytic comparison models. In computing the GW spectrum from the Ketju runs we use the semi-analytic method of Section 4.2 until the SMBH binary semimajor axis is ∼100R_s, after which we switch to the direct waveform method of Section 4.3 and compute the isolated evolution of the binary for the final a ≲ 100 R_s. The exact point where the switch happens depends on the eccentricity of the binary, as all the major harmonics of the initial orbital frequency need to be visible in the spectrum computed from the waveform. Based on Section 4.4, the error of the semi-analytic GW method should be of the order of a few percent at the switchover point, as the error scales with ${\epsilon }_{\mathrm{PN}}\propto {a}^{-1}$. For computing the GW signal from the semi-analytic Peters and Peters-Quinlan models, we naturally utilize the semi-analytic method throughout.

The energy spectrum of GWs as emitted over the lifetime of the binary is shown in Figure 7. The spectra are all fairly similar due to the similar orbital evolution, and show the characteristic peaked shape of an initially highly eccentric binary. The differences between the runs are mainly explained by the different total black hole masses and mass ratios, as well as the different eccentricities.

Zoom In Zoom Out Reset image size

At high frequencies, where the binaries have mostly circularized (e ∼ 0), the spectra are fairly similar. For runs A–D the full PN FFT computation of the spectrum gives results that differ by only a few percent from the results computed using the Peters evolution combined with the semi-analytic GW method until just before the highest modeled frequencies when the binary is at a separation of R = 6R_s where the Ketju simulation is stopped. The Peters-Quinlan model gives differences of a similar magnitude for these runs, which is expected as in this phase these massive binaries are highly relativistic and the differences in the spectrum are mainly due to PN effects. The effect of the environment in the spectrum is larger for run X, where the Peters model shows a difference of almost 20% at $f=0.2\,{\mathrm{yr}}^{-1}$, while the Peters-Quinlan model agrees at the percent level. This is mainly due to the larger eccentricity of the binary in this frequency range, and in particular the difference in eccentricity in the Peters model, which shifts the peak of the spectrum. The agreement at slightly higher frequencies where the binary has circularized is much better, with the semi-analytic models matching almost exactly the direct waveform computation results. At the point when the GW calculation is switched to the direct waveform method, the differences to the semi-analytic GW method are of order a few percent, consistent with the PN expansion parameter being ${\epsilon }_{\mathrm{PN}}\sim 0.01$ at that point. This leads to the small jump that can be seen in the difference to the semi-analytic comparison models (see the zoom-ins in Figures 7 and 8).

Zoom In Zoom Out Reset image size

At lower frequencies where only the semi-analytic GW computations are performed on the Ketju simulations, the isolated Peters models lead to larger amplitudes with peaks at lower frequencies. The resulting relative differences are mostly around 50%, but can reach values of up to 300% in the case of run X. These large differences are due to the faster hardening of the binary due to the stellar environment, which is not accounted for in the Peters models, and in fact the Peters-Quinlan models which include this effect are in much better agreement with the Ketju results. The most significant differences between the Ketju and Peters-Quinlan results occur in the case of run A, where the Peters-Quinlan result is about 70% lower. In this case the hardening parameter used for the Peters-Quinlan model was already found to be 5–8 times larger than the effective ones determined from the Ketju run, so the large difference is not surprising. The other runs (B–D) with large differences in the hardening parameters show similar behavior, with differences of around 40%. In any case the differences occur mainly at relatively low GW frequencies that are unobservable for the foreseeable future.

The differences in the GW spectrum are mainly due to differences in the gravitational wave timescale τ_GW, shown in Figure 8, and also in part due to the differences in orbital eccentricity. The evolution of eccentricity determines how the emitted GWs are spread out over the spectrum, which also determines the location of the peak of the spectrum. Here too the effects of the environment are obvious, with the Peters models having longer and the Peters-Quinlan models shorter GW timescales than the Ketju runs. The differences are mostly between 10% and 60% at orbital frequencies below $0.1\,{\mathrm{yr}}^{-1}$, with run X showing both the highest difference of ∼150% to the Peters model and almost no difference to the Peters-Quinlan model. At high frequencies both of the semi-analytic Peters formulae–based comparison models show slightly longer timescales, which is due to the inclusion of higher order PN effects in Ketju.

The good agreement at the percent level between the spectra computed using the Peters formulae and the full PN computation at high frequencies may seem surprising, since the strongest effects of the PN corrections should appear there. If we look at the instantaneous GW flux, shown in Figure 9, the difference between the methods does indeed increase with increasing frequency, reaching tens of percent. However, this is countered by the slightly slower evolution of the orbit according to the Peters formula, which can be seen as a slightly longer timescale at high frequencies in τ_GW in Figure 8. Fundamentally, the similarity of the spectra computed with the different methods for an isolated binary follows from the conservation of energy, which is exact by construction in the Peters evolution and approximate in the PN3.5 evolution, as was also noted in Section 4.4. Figure 9 also shows the importance of including the waveform PN1 corrections as well, as otherwise the flux can be overestimated by over 10%.

Zoom In Zoom Out Reset image size

Another point relating to the energy conservation of the system when using PN dynamics to compute the binary evolution and GW emission in the final stages of an SMBH coalescence is the fact that the energy of the system is conserved only approximately. This limits the applicability of the method, as large violations of energy conservation clearly signify incorrect results. Figure 10 shows the total relative error in the energy conservation of the system during the final stages of the merger. The orbital energy is computed using the PN-accurate energy of Equation (7), while the emitted GW energy is computed by integrating the instantaneous energy flux computed from the waveform. For most part the full computation including the PN waveform corrections gives an error of under 5%, and only in the very final stages does the error begin to rapidly increase. This suggest that the results computed using this method are fairly reliable down to around a ∼ 10R_s or to orbital frequencies of ${f}_{r}\sim ({10}^{10}{M}_{\odot }/M)\,{\mathrm{yr}}^{-1}$, which also corresponds to the part in the spectrum in Figure 7 where the direct waveform and the Peters results begin to again diverge. When the GW emission is computed using only the quadrupole formula the error is naturally greater. A similar error in energy conservation for the full PN computation has also been observed by Csizmadia et al. (2012).

Zoom In Zoom Out Reset image size

Since the eccentricity affects strongly both the merging time of the SMBHs (Peters 1964) and the emitted GW spectrum (Enoki & Nagashima 2007), understanding the distribution of eccentricities of SMBH binaries is essential for correctly predicting the emitted GWB. As can be seen from Table 1 all the Ketju runs studied in this paper show high initial eccentricities in excess of e > 0.9 at the point where GW emission and other relativistic effects start to become important. The encounter orbits of our galaxies are nearly parabolic as motivated by cosmological simulations (Khochfar & Burkert 2006) and high initial eccentricities have also been found in earlier studies employing Ketju (Rantala et al. 2017, 2018).

Previous works studying the evolution of SMBH binaries in a galaxy merger using a smaller number of particles have also in general found high binary eccentricities, e≳0.8 (e.g., Berentzen et al. 2009; Khan et al. 2011, 2018a; Preto et al. 2011). The initial binary eccentricity is affected by the initial condition setup, with strongly radial merger orbits producing in general higher eccentricity binaries. Another factor that could affect the eccentricity distribution is numerical resolution, as higher resolution simulations tend to produce higher eccentricities with less scatter (see Figure 12 in Rantala et al. 2017). Finally, the interaction of circumbinary gas will likely affect the eccentricity evolution (e.g., Roedig et al. 2011), with this effect being particularly important for somewhat lower mass SMBHs $({M}_{\bullet }\lesssim {10}^{8}{M}_{\odot }$).

It is important to account for the SMBH binary eccentricity when computing predictions for PTA observations. In older work this was not typically done, and instead 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 the case of a fixed initial eccentricity as well as other effects such as stellar scattering and drag from a circumbinary accretion disk using semi-analytic models (e.g., Kelley et al. 2017b; Sesana et al. 2018).

These semi-analytic models appear to work quite well for the mergers considered in this paper, as far as the processes modeled by Ketju are concerned. The differences in the evolution of orbital parameters (see Figures 4 and 5) and the GW spectrum in the PTA frequency range (see Figure 7) are brought down to under 5% provided a correct choice of semi-analytic model parameters, and even when using incorrect parameter values or ignoring the stellar scattering completely the differences remain typically at a level of few tens of percent.

We find that the measured hardening parameters of the considered semi-analytic scattering model can be smaller for our Ketju cored galaxy simulations than common literature values by factors up to 8. However, this is partially offset by the PN3.5 accurate equations of motion also speeding up the merger process slightly compared to models including only the leading quadrupole radiation reaction. The difference to literature values is also smaller for the simulation without a stellar core that has less extreme SMBH masses (run X), suggesting that for most SMBH binaries the commonly used semi-analytic models and their parameters are reasonably accurate.

Comparing the eccentricities seen in our simulations to the results of Kelley et al. (2017b), we find that while our eccentricities are among the highest considered there they are not quite high enough to cause the strong suppression of the GWB seen for their model with an initial eccentricity of e ∼ 0.99.

For PTA observations the uncertainties for the contribution from very massive binaries are probably dominated by discretization effects, as there are likely to be only a few very massive systems in the final stages of a merger, with high enough GW frequencies to be observable. These effects affect also lower mass binaries, and for a realistic population of SMBHs, discretization effects tend to become significant at frequencies above $f\sim 0.3\,{\mathrm{yr}}^{-1}$, depending on the eccentricities of the binaries (Kelley et al. 2017b).

Whether the error caused by the simplified semi-analytic models leads only to a slightly increased statistical uncertainty, or if some kind of systematic bias is introduced cannot be answered with our limited sample of runs. However, one should keep in mind that the good agreement between the semi-analytic models and the full Ketju calculation was achieved when using the properties of the stellar population derived directly from the Ketju runs. Thus the strong effects of the SMBH binary on the stellar population are accounted for in the density and velocity dispersion, while in the typical use case of semi-analytic models they are not, as this regime falls below the gravitational softening length. This issue will be particularly important for cosmological simulations in which the stellar population cannot be resolved at high spatial accuracy.

At high gravitational wave frequencies affected by discretization effects, the effects of higher order PN corrections can also become significant. This is due to the fact that the GW spectrum is determined by the instantaneous GW flux when there are only a few sources, and as was seen the difference between the PN-accurate instantaneous flux and the commonly used quadrupole approximation is of the order of 30% at $f\sim 1\,{\mathrm{yr}}^{-1}$ for M ∼ 10¹⁰M_☉ binaries (see Figure 9). Ignoring these effects can thus lead to a systematic bias in the computed background amplitude, if there are enough massive binaries for which this effect is at its strongest. However, we note that the corrections to the GW flux can also be handled analytically to at least PN3 order (e.g., Arun et al. 2008; Blanchet 2014).

The derived differences between the accurate Ketju computations and the simpler semi-analytical models are relatively small, on the level of ∼10%. Thus, it is important to study various potential caveats to the numerical simulations in order to assess whether differences on this level can be reliably determined. First, the stellar environment is not actually resolved, with each stellar particle representing about ∼10⁵ real stars. Second, the motion of the binary is modeled only approximately using PN3.5 level corrections for the binary only, and while the accuracy of this modeling is state of the art for the post-Newtonian approach, the actual physics may differ from this.

The hybrid nature of Ketju allows us to utilize a larger number of particles than in standard N-body codes. However, despite this fact, the actual number of stellar particles interacting with the SMBH binary in its final stages before the merger is still small compared to the actual number of stars that they represent. When the SMBH binary merges in our simulations, there are typically only a few stellar particles left within $10\,\mathrm{pc}$ of the binary instead of tens or hundreds of thousands. As a result, strong interactions between the binary and stars happen far more rarely, and when they do happen they are too strong, leading to step-like changes in the binary orbital parameters. Fortunately, it appears that most of the environmental effects are caused by fairly weak long-range interactions, thus the error in modeling the strong interactions should not seriously affect the long-term binary evolution.

The accuracy of the binary motion also depends on the accuracy of the PN3.5 equations of motion. Here we are limited by the fact that higher order equations of motion are unavailable for general binaries. Such higher order corrections would be further complicated by their dependence on the history of the binary in addition to its instantaneous state, making implementing them impractical.

While formally the higher order corrections should be of order ${ \mathcal O }({\epsilon }_{\mathrm{PN}}^{4})$, in practice the effects may be significant even for fairly small ${\epsilon }_{\mathrm{PN}}$. This was illustrated by the rather large change in the spectrum when going from PN3 to PN3.5 equations of motion in Figure 3, even though formally the effect is of order ${\epsilon }_{\mathrm{PN}}^{3.5}\lesssim {10}^{-3}$ there. The apparently fairly significant contributions from the unknown higher order corrections are also seen in the failure of energy conservation when the binary is very near merger (see Figure 10). Since the eccentricity of the binary is fairly small in this phase, results computed for circular binaries can be used to estimate just how much the higher order corrections affect the evolution of the binary. According to Blanchet et al. (1995), the nonlinear tail-terms that are the leading correction neglected here cause the number of orbits an equal-mass binary completes during its last decade in orbital frequency to change by almost 10%, which can be considered to be a fairly substantial effect.