Issue 
A&A
Volume 615, July 2018



Article Number  A71  
Number of page(s)  9  
Section  Astrophysical processes  
DOI  https://doi.org/10.1051/00046361/201730489  
Published online  17 July 2018 
Gravitational wave driven mergers and coalescence time of supermassive black holes
^{1}
Department of Space Science, Institute of Space Technology,
PO Box 2750,
Islamabad,
Pakistan
email: khanfazeel.ist@gmail.com
^{2}
Main Astronomical Observatory, National Academy of Sciences of Ukraine, MAO/NASU,
27 Akad. Zabolotnoho St.,
03680
Kyiv,
Ukraine
email: berczik@mao.kiev.ua
^{3}
National Astronomical Observatories and Key Laboratory of Computational Astrophysics, Chinese Academy of Sciences, NAOC/CAS,
20A Datun Rd.,
Chaoyang, Beijing
100012,
PR China
^{4}
Astronomisches RechenInstitut, Zentrum für Astronomie der Universität Heidelberg,
Mönchhofstr. 12–14,
69120
Heidelberg,
Germany
email: just@ari.uniheidelberg.de
Received:
24
January
2017
Accepted:
12
March
2018
Aims. The evolution of supermassive black holes (SMBHs) initially embedded in the centres of merging galaxies realised with a stellar mass function (SMF) is studied from the onset of galaxy mergers until coalescence. Coalescence times of SMBH binaries are of great importance for black hole evolution and gravitational wave detection studies.
Methods. We performed direct Nbody simulations using the highly efficient and massively parallel phiGRAPE+GPU code capable of running on highperformance computer clusters supported by graphic processing units (GPUs). PostNewtonian terms up to order 3.5 are used to drive the SMBH binary evolution in the relativistic regime. We performed a large set of simulations with three different slopes of the central stellar cusp and different random seeds. The impact of a SMF on the hardening rate and the coalescence time is investigated.
Results. We find that SMBH binaries coalesce well within one billion years when our models are scaled to galaxies with a steep cusp at low redshift. Here higher central densities provide a larger supply of stars to efficiently extract energy from the SMBH binary orbit and shrink it to the phase where gravitational wave (GW) emission becomes dominant, leading to the coalescence of the SMBHs. Mergers of models with shallow cusps that are representative of giant elliptical galaxies having central cores result in less efficient extraction of the binary’s orbital energy, due to the lower stellar densities in the centre. However, high values of eccentricity witnessed for SMBH binaries in such galaxy mergers ensure that the GW emission dominated phase sets in earlier at larger values of the semimajor axis. This helps to compensate for the less efficient energy extraction during the phase dominated by stellar encounters resulting in mergers of SMBHs in about 1 Gyr after the formation of the binary. Additionally, we witness mass segregation in the merger remnant resulting in enhanced SMBH binary hardening rates. We show that at least the final phase of the merger in cuspy lowmass galaxies would be observable with the GW detector eLISA.
Key words: galaxies: interactions / galaxies: nuclei / quasars: supermassive black holes / gravitational waves
© ESO 2018
1 Introduction
Gravitational waves (GWs) were predicted by Einstein shortly after presenting his general theory of relativity (GR; Einstein 1916, 1918). Indirect evidence that GWs exist came from studies of the orbital decay of binary pulsars in accordance with GR (Hulse & Taylor 1975). Recently, the first direct measurement of GWs was accomplished by the observation of the merging event of the stellar mass black hole binary GW150914 with LIGO (Abbott et al. 2016). Since then, three more events including a neutron star–neutron star merger have been observed with LIGO and VIRGO (Abbott et al. 2017). Observations of GWs coming from various sources at various cosmic epochs will open an entirely new window to study the Universe, currently beyond the capabilities of electromagnetic probes. Supermassive black hole (SMBH) binaries are considered promising sources of GWs (Begelman et al. 1980). Observations of GWs emitted during the final phase of inspiral would provide the merger history of galaxies as a function of redshift leading to important constraints on SMBH and galaxy formation and evolution scenarios. The International Pulsar Timing Array (IPTA) is searching for GWs from coalescing binary SMBHs in the mass range 10^{7}−10^{9} M_{⊙} up to redshift z = 2 (Hobbs et al. 2010; Desvignes et al. 2016; Reardon et al. 2016; The NANOGrav Collaboration 2015; Verbiest et al. 2016). In the future the (evolved) Laser Interferometer Space Antenna (eLISA, LISA) is expected to detect GWs to much larger redshifts (z ~ 10) (AmaroSeoane et al. 2013, 2017; Gravitational Observatory Advisory Team 2016).
SMBH binaries form as a result of mergers between two sizable galaxies, each hosting a central SMBH, which are ubiquitous in galaxy cores (McConnell & Ma 2013; Kormendy & Ho 2013). Observational searches for SMBH binaries are going on; there are a handful of cases of two wellseparated accreting SMBHs seen as AGNs, as well as circumstantial evidence for bound Keplerian binaries (Komossa et al. 2003; Bogdanović 2015; Graham et al. 2015). In the merger remnant, the evolution of SMBHs happen in three phases, each characterised by a distinct physical process. Dynamical friction (Chandrasekhar 1943; Just et al. 2011) is responsible for the initial sinking of SMBHs in the merger remnant bringing the SMBHs close enough that they form a binary system. Dynamical friction efficiently extracts energy from the binary orbit until the binary gets hard. By this point in time, dynamical friction ceases to be efficient while SMBHs are separated by parsec scale distances. In the second phase, stars on orbits intersecting the binary orbit extract energy from the SMBH binary by the slingshot mechanism during threebody encounters bringing the black holes closer. If it is efficient to bring the SMBHs in the binary close enough (milliparsec separations), GW emission in the third and final phase drains out the remaining energy in the binary orbit leading the SMBHs to coalesce. How efficiently the SMBH binary evolves in the threebody scattering phase such that the separation between the SMBHs shrinks to the GW dominated regime depends strongly on the orbit contents of the host galaxy (merger remnant; Merritt & Poon 2004; Li et al. 2015; Gualandris et al. 2017). Earlier it was shown that for galaxy shapes close to sphericity, SMBH binaries may stay longer than a Hubble time in the threebody scattering regime (Makino & Funato 2004; Berczik et al. 2005). However, strongly flattened or mildly triaxial shapes, a natural product of galaxy mergers, have shown an effective shrinking of the SMBH binary’s semimajor axis to the point where GW emission dominates (Berczik et al. 2006; Khan et al. 2011, 2013; Preto et al. 2011; Gualandris & Merritt 2012; Vasiliev et al. 2015). The hardening rate and eccentricity of the binaries depend strongly on the central stellar density profile and hence the estimated coalescence times of SMBHs in binaries (Khan et al. 2012). In previous studies of SMBH binary evolution in equilibrium galaxy models or in mergers of galaxy bulges, coalescence times were obtained (Khan et al. 2015; HolleyBockelmann & Khan 2015; Sesana & Khan 2015; Rantala et al. 2017) by extrapolating the nearly constant hardening rate of SMBH binaries in the threebody scattering phase to theGW dominated regime and then using orbit averaged expressions (Peters & Mathews 1963) for the hardening by GW emission.
Here we study the SMBH binary evolution in mergers of galaxy spheroids having various stellar density profiles towards the centre, as in Khan et al. (2012) but this time following the binary evolution into the relativistic regime for a complete set of mergers by including postNewtonian (PN) terms up to order 3.5 in the equation of motion of the SMBH binary. The effect of a stellar mass function (SMF) on the SMBH binary evolution has not been well studied. On the one hand, it is known from threebody scattering investigations (Hills & Fullerton 1980; Sesana 2010) that for a uniform stellar population the hardening rate should be independent of the stellar mass in the lowmass regime and reduced for higher intruder masses above 1:10 with respect to the secondary SMBH mass. On the other hand, a mass segregated system should have an enhanced hardening rate because the velocity dispersion of the highmass end is smaller, which leads to an enhanced contribution to the energy extraction. In order to investigate this effect we introduced a SMF for particles in the merging galaxies in order to allow for mass segregation effects.
This paper is arranged as follows: Sect. 2 describes our models and their scalings. It also includes numerical codes and hardware used to perform the galaxy merger and SMBH binary evolution simulations. Sect. 3 presents the results of our study. Convergence tests are described in Sects. 4. and 5 summarises and concludes our study.
2 Simulation setup and numerical techniques
2.1 Galaxy models
We set up our initial galaxies by spherical isotropic distributions of stars such that the density distribution satisfies a Dehnen (1993) profile (1)
where M_{gal} denotes the mass of the galaxy, r_{0} its scale radius, and γ the slope of the inner density profile. We generated galaxy models for three different values of γ, γ = 0.5, 1.0, and 1.5 (series A, B, C, respectively) in order to cover the observed range of density profiles of galactic nuclei. A particle having mass equal to 1% of the mass of the galaxy is placed at the centre to represent the SMBH. This SMBH mass (M_{∙}) is a few times greater than the observed M_{∙} –M_{gal} relation. Therefore, our models can be viewed as a representation of the central parts of galaxies/bulges. Positions and velocities are assigned to the stars by numerically computing the distribution function in the combined potential of black hole and stars such that our models are in dynamical equilibrium. In our “model units” M_{gal}= G = r_{0} = 1 for the primary (i.e. more massive) galaxies. The masses of the smaller secondary galaxies and their SMBHs are scaled down by a factor q. The smaller galaxies follow the same density profile as the primary galaxies, but have smaller masses and different scale radii r_{0,s}. The size ratio of the galaxies and the corresponding mass ratio q are related as . The primary and secondary galaxies both have the same number of particles N. For this study we chose q = 0.25 and N = 200 000. The SMBH binary evolution in the merger remnant realised with 400 000 stellar particles is expected to be Nindependent (Berczik et al. 2006; Preto et al. 2011; Khan et al. 2011). The basic setup is very similar to the models in Khan et al. (2012), but with slightly smaller number of particles and larger SMBH mass in order to speed up the simulations.
For each γ we generated three galaxy models with different random initialisations (seeds) for both major and minor galaxies. In order to study the effect of a SMF and the corresponding mass segregation on the SMBH binary evolution each galaxy model was given a mass function according to a Salpeter IMF in the mass range of 0.08 to 8 M_{⊙} with a mean mass of 0.25 M_{⊙}, also with three different random realisations (see Table 1). In our simulations each particle represents a number of stars with same velocity and mass. Since the evolution times of the SMBHs are of the order of 1 Gyr and we are mainly interested in the possible effect of mass segregation, we did not include mass loss by stellar evolution. For mass segregation the dynamic range of particle masses is the most important parameter. Our choice of a factor of 100 between most and least massive particles corresponds to the mass range of stellar black holes and the lower main sequence of an old population. In this sense the chosen mass function should be seen as a very rough representation of the mass function of a real galaxy.
In each galaxy the mass ratio of the SMBH and the most massive stellar particle is 1:62 and after the merger the maximum mass ratio of the secondary SMBH and the maximum stellar mass is 1:15.5, which is still in the limit of small intruder mass for the threebody scattering events. In the test simulations discussed in Sect. 4 the number of particles is increased by a factor of 5, leading to a maximum mass ratio of 1:77, which is sufficient for quantifying the hardening rate but still problematic for the investigation of the eccentricities due to large fluctuations.
For each γ value (series A, B, C), we have ten runs, the 0th run represents a galaxy merger of galaxies having an equal mass of stellar particles in each galaxy. The remaining nine runs are three galaxy models each having three different random realisations of the SMF (see Table 1).
Galaxy merger runs.
Physical scalings for our galaxy models.
2.2 Initial orbits and scaling to real galaxies
For each run we set two galaxies (primary and secondary galaxy) at apocentre on eccentric orbits with e = 0.75. The initial separation between the centres of the merging galaxies in our simulations is 15 of our model units.
We chose three different galaxies, M87, M31, and the Milky Way (MW), for the physical scaling of our models. We used the observed mass of the SMBH and the velocity dispersion in these galaxies to calculate the sphere of influence r_{h} of the SMBH. Then we compared our model SMBH and its sphere of influence to the observed ones to get the physical scaling of our models (see Table 2). M87 represents giant elliptical galaxies, which typically have very shallow cusps or a core in the centre, and thus represent our γ = 0.5 (series A) shallow cusp models. To scale our series B, which have a γ = 1.0 inner density slope, we chose M31. The physical parameters of the MW centre were used to scale our γ = 1.5 steep cusp models of series C. The details of the scaling parameters are given in Table 2.
2.3 Numerical code
The numerical simulations of the galaxy mergers are performed using an updated version of the direct Nbody code ϕGRAPE, originally designed to run on GRAPE cards. Our updated code (ϕGRAPE+GPU) is capable of running on massively parallel clusters supported by graphic processing units (GPUs). For the pairwise force calculations we use a softening parameter equal to 10^{−5} in model units for the stars and no softening for the SMBHs. For the pairwise forces we apply the rms values of the softening leading to a softening of 7 × 10^{−6} for star–SMBH interactions and no softening for the SMBH–SMBH interaction.
Relativistic effects are taken into account by incorporating PN terms up to order PN3.5 in the SMBH binary’s equation of motion (Blanchet 2006). More details on the simulation code can be found in Khan et al. (2013). We used the Laohu cluster of the National Astronomical Observatories of the Chinese Academy of Science to perform our simulations.
3 SMBHbinary evolution
We discuss first the hardening of the SMBH binaries. Figure 1 shows how the separation between the SMBHs shrinks, initially due to the galaxy merger, later on as a SMBH binary forms and hardens due to dynamical friction, then in the hard binary phase due to threebody scattering, and finally due to GW emission. As the galaxies merge, the separation between the two black holes shrinks below 1 model unit. Galaxies are merged at almost the same time at T ~ 100 after starting the simulations for all models because of the same masses and orbits. Then we witness a fast decay in the SMBH binary separation due to dynamical friction. Dynamical friction is more efficient for steep cusps due to a higher central density resulting in a faster orbital decay for case C than for case A. However, dynamical friction becomes inefficient when the orbital velocity of the binary is significantly larger than the velocity dispersion, i.e. when the binary gets hard. In case C this happens at a smaller separation due to the larger velocity dispersion. We can see that for the models in group A the dynamical friction becomes less efficient at a separation between the two SMBHs of roughly 0.01 in model units, whereas for models C the same happened at a separation that was ten times smaller.
As a consequence of these two competing effects the transition to the threebody scattering phase takes place at roughly the same time T ~ 150. The oscillations of the separation due to the eccentric SMBH binary orbit are not fully resolved in the plots due to the short orbital time compared to the larger output timesteps.
The evolution of the inverse semimajor axis of the SMBH binaries (which is a measure of the binding energy) for all our merger models is shown in Fig. 2. We see that the inverse semimajor axis of the binaries evolve at a constant rate due to threebody scattering. The final phase, where the energy loss by GWs dominates, is characterised by an accelerated hardening. The onset of this phase and the final merging times are consistent with the analytic estimates combining a constant hardening rate s_{3body} and the orbit averaged hardening rate from Peters & Mathews (1963).
We estimated the SMBH binary hardening rates in the stellar dynamical hardening regime by fitting straight lines to a^{−1} (t) (see Table 3). The hardening rates are systematically higher for steeper cusps with higher values of γ as already shown in Khan et al. (2012).
The last rows in Table 3 show the mean hardening rate ⟨s⟩ and the rmsscatter of the simulations with SMF for each series. The hardening rates are approximately 30− 40% higher for SMBH binaries evolving in merger remnants formed as a result of mergers of galaxies having a SMF compared to the singlemass simulations A0, B0, and C0 without SMF.
If the stellar population is well mixed in phase space, there is no impact of a SMF on the hardening rate expected. The hard binary phase, where threebody encounters dominate the energy extraction from the SMBH binary, corresponds to the lowvelocity limit for the intruders in threebody scattering. For lowmass perturbers (mass ratios below 1:10 with respect to the secondary SMBH) the mean energy loss of the SMBH binary is proportional to the intruder mass (see e.g. Hills & Fullerton 1980; Quinlan 1996; Sesana 2010). For higher mass ratios, the energy loss is sublinear, leading to a reduced hardening rate. As a consequence the hardening rate should depend only on the mass density distribution and the kinematic properties and scales as (2)
with the dimensionless hardening parameter H (Sesana & Khan 2015). Density ρ_{f} and velocity dispersion σ_{f} are to be taken at the radius r_{f}, usually the influence radius of the binary (e.g. Sesana 2010). In a first test to find the reason for the enhanced hardening rate, we compared the mass density and velocity dispersion profiles as well as the anisotropy profiles of the singlemass and the SMF simulations, but did not find any differences above the noise level. However, the number density profiles are different, which shows that the mean particle mass increases with decreasing distance from the SMBHs. In order to see whether these increased hardening rates witnessed in SMF runs are due to mass segregation, we investigated the mass profiles of particle species of different mass ranges. We chose three mass species of particles (in numbers): the most massive 12% (M), intermediate 30% (I), and the least massive 58% (L). For Fig. 3 we selected representative runs (A1, B2, and C1) of each γ in our merger simulations. The plots show the fractional contribution of each stellar mass bin to the cumulative mass profile for two different times. We can clearly see that the three species have very different mass profiles showing mass segregation for the most massive species (the estimated mass segregation timescale is ~ 100 time units). The steeper density profile of the highmass stars goes hand in hand with a smaller velocity dispersion and vice versa for the lowmass stars. If we apply Eq. (2) to each mass component separately and add up the contributions to the total hardening rate (adopting a universal H), we find a slightly higher value for s, since the ratio ρ_{f} ∕σ_{f} is lower than the sum ρ_{i}∕σ_{i}. This effect is very small, however, and does not explain the enhanced hardening rate for the SMF case.
It is well known that the quantification of the eccentricity evolution, which is connected to the angular momentum of the binary (e.g. Mikkola & Valtonen 1992; Quinlan 1996; Sesana 2010), is much harder than for the binding energy evolution. The main reason is that angular momentum changes are firstorder perturbations leading to a high sensitivity of the eccentricity on the random properties of the individual encounters. The evolution of the SMBH binary eccentricities e is shown in Fig. 4. In Table 3 the mean eccentricity for each run is listed, as well as averages over the nine SMF runs for each series. We observe strong fluctuations of the individual eccentricities, which are significantly larger for the SMF runs. A detailed discussion is given in Sect. 4. We note that the mean eccentricities of the SMF runs show a large scatter in each series. The eccentricities of the singlemass runs A0 and B0 are at the higheccentricity end, whereas for the steep cusp C0, the eccentricity is at the low end. We do not observe a correlation of the hardening rate with the eccentricity, but for the shallow and intermediate cusp series A and B the coalescence times depend strongly on (1 − e^{2}) as expected from the inset of energy loss by GW emission. For the low eccentricities of case C, the impact of variations in e are small.
We have plotted histograms of the average eccentricities in bins of 0.2 for all models in Fig. 5. For models A, most binaries (7 out of 10) have a very high average eccentricity (in the range 0.8–1.0), whereas for models B, this number becomes 4 for the 0.8–1.0 and the 0.6–0.8 bins. For both series A and B, the number of SMBH binaries in the eccentricity range 0.4–0.6 is 2, whereas there is no binary with an eccentricity below 0.4. The situation is very different for runs C, where the merging galaxies have steep inner profiles (dense nuclei) with γ = 1.5. The number of SMBH binaries peaks (5 out of 10) in the bin 0.2–0.4. There are three binaries in the bin 0–0.2 and two in the bin 0.4–0.6. This systematic trend to more circular binary orbits for steeper cusps is in accordance with the more effective circularisation by dynamical friction in steeper density profiles.
The coalescence times T_{coal} for SMBH binaries are also given in Table 3. The average coalescence time of SMBH binaries in the SMF runs A is 1.25 Gyr. For SMBH binaries in runs B and C it is 0.39 and 0.31 Gyr, respectively. In series A and B the rms variation of ~25% is similar, whereas in series C the scatter is much smaller. The main reason for this variation is the scatter in eccentricity in the higheccentricity regime.
In the late phase of the evolution the SMBH binaries emit lowfrequency GWs, which may be observable with GW detectors like various pulsar timing array experiments (Desvignes et al. 2016; Reardon et al. 2016; The NANOGrav Collaboration 2015; Verbiest et al. 2016) and planned spaceborne GW observatory eLISA or LISA (AmaroSeoane et al. 2013, 2017; Gravitational Observatory Advisory Team 2016). In Fig. 6, we have calculated for the three selected cases A1, B2, and C1 the characteristic strain adopting a redshift of z = 3.0 with a corresponding luminosity distance of D = 26 Gpc. The lowmass SMBHs (in steep cusps) should be visible at the lowfrequency end of eLISA and LISA. The highmass end (with shallow cusps) is close to the frequency range and sensitivity of current PTAs.
Fig. 1
Shrinking of the separation between the SMBHs in galaxy mergers (top panel for models A, middle panel for models B, and bottom panel for models C). Separation and time are in model units. 
Fig. 2
Evolution of the inverse semimajor axis (top panel for models A, middle panel for models B, and bottom panel for models C). The singlemass model 0 is represented with a thicker line in each plot. 
Galaxy merger runs.
Fig. 3
Fractional contribution of each stellar mass bin (M, I, and L for massive, intermediate, and low mass, respectively) to the cumulative mass profiles at an early (T = 150) and late (T = 400) evolution time for representative models A1 (top), B2 (middle), and C1 (bottom). 
Fig. 4
Evolution of SMBH binary eccentricities (top panel for models A, middle panel for models B, and bottom panel for models C). Model 0 is represented with a thicker line. 
4 Convergence tests
We note that there are a few sudden jumps in the 1∕a evolution, especially for steeper cusp mergers. Sometimes there is also a jump in eccentricity, but not always. Our analysis shows that occasionally the SMBH binary and a massive stellar particle at the high end of the SMF form a threebody bound system surviving for a relatively long time. In this phase the slingshot ejection of a fourth body is able to produce these large jumps in 1∕a. To determine the hardening rate, we have corrected for these unrealistic jumps by measuring s in time intervals with no significant jumps. In the eccentricity e evolution there are more jumps in the shallow cusp case A. These jumps do not dominate the overall evolution. Most of the fluctuations occur on longer timescales covering many orbital times and many scattering events. To reduce the effect of this random scatter significantly, much higher particle numbers may be necessary. In this section we test the N dependence of the hardening rate and the fluctuations in the eccentricity.
In order to test to what degree our results depend on the particle number N, we performed aconvergence test for two of our runs where we use up to 2 million particles. We chose a shallow profile of γ = 0.5 because it is expected that the shallow profile mergers take less computational time to complete. We chose two series of runs: A0, A01m, and A02m are singlemass runs with 400k, 1 million, and 2 million particles, and A4, A41m, and A42m are SMF runs with 400k, 1 million, and 2 million particles. For our largest N we have a five times lower maximum mass ratio of the secondary SMBH to the stars of 1:77. The results of this study are presented in Fig. 7. The top panel shows the inverse semimajor axis and the middle panel the corresponding smoothed hardening rates. In the early phase up to T = 150, where dynamical friction is still active and the inner loss cone is not empty, there are differences in the evolution. Here we are interested in the stationary later phase where the threebody encounters dominate the hardening.
The Nindependent evolution of 1∕a witnessed by earlier studies for galaxy mergers is reproduced for the singlemass runs. The stationary phase of constant hardening rate is reached at about T = 200. For the SMF runs a weak trend of decreasing s with increasing N may be present. However, the hardening rate of the 1 million run A41m show fluctuations breaking this trend at T = 200. In addition to this possible trend, we also notice that s is systematically higher for all SMF runs when compared with their counterpart singlemass runs. Since the twobody relaxation time, and thus the mass segregation timescale, depends on N, we expect adependence (delay) of mass segregation and enhanced hardening with increasing N. For the SMF cases we still see a considerable mass segregation at the highmass end as in the 400 k case, but starting a bit later at T = 150 consistent with the enhanced hardening rate. For the 1 and 2 million runs we observe the same density, velocity dispersion, and anisotropy profile as for the 400 k case. There is no significant difference of these global properties in phase space for the SMF and nSMF cases. We calculate the dimensionless hardening parameter H using Eq. (2) for our 2 million particle runs. Parameter H has a value of 14.82 for run A42m and a somewhat lower value of 12.2 for the A02m run, which reflects a similar difference in hardening rates in the two runs. We conclude that mass segregation has a more subtle effect on the distribution functions enhancing the encounter statistics in order to explain the higher hardening rates.
Another important parameter for the merger times is the SMBH binary eccentricity (bottom panel of Fig. 7). For all singlemass runs there is a remarkable consistency in SMBH binary eccentricities of ~ 0.9 in the late stationary phase, even though for the 400 k case A4 the eccentricity starts with a smaller value and increases as theoretically predicted. For the SMF runs we do not witness a strong match, but eccentricities of all binaries are in the most favourable range of 0.6–1.0 for A runs, and there is no systematic trend with increasing particle number. We note that the fluctuations in the eccentricity are not reduced significantly by increasing the particle number by a factor of five, although jumps are now much less prominent and the evolution of e is much smoother. This is a hint that they may arise from inhomogeneities in the distribution of particles leading to an anisotropic flux of interacting stars. In order to reach the same maximum mass ratio of 1:500 as in the singlemass mergers, about 15 million particles would be needed, which is not feasible with the current computer facilities.
Fig. 5
Histograms showing the number of SMBH binaries in eccentricity bins with width 0.2. 
Fig. 6
Comparison of the predicted characteristic strain for models A1, B2, and C1 at redshift z = 3.0 and the sensitivity curves for various GW detectors. 
Fig. 7
Evolution of SMBH binary inverse semimajor axis, hardening rates, and eccentricity for our convergence test runs with N = 400k, 1 and 2 million for the singlemass A0, and the SMF case A4. 
5 Summary and conclusion
We have performed a statistical set of mergers of galaxies with a mass ratio of q = 1∕4 with shallow, intermediate, and steep central density profiles. Each galaxy contains a central SMBH with 1% of the galaxy mass and the particles were realised with a Salpeterlike SMF. We have used independent random realisations for the initial phase space distribution and for the stellar masses. The dynamical evolution of the galaxy pairs was performed with a direct Nbody code including general relativistic effects of the SMBHs using PN corrections. All simulations were performed until the final coalescence of the SMBHs by GW emission.
The total coalescence time is dominated by the length of the threebody encounter phase, where the SMBH hard binary lose energy by slingshot encounters with stellar particles. However, the end of this phase depends strongly on the eccentricity of the SMBH binary by the onset of GW emission. The hardening rate s, the mean eccentricity e, and as a consequence the coalescence time T_{coal} show a significant scatter due to the random realisation of both the phase space distribution and the SMF. The scatter and fluctuations in the eccentricity does not decrease significantly for the tested shallow cusp case when increasing the number of particles by a factor of five for the SMF case. In contrast, the singlemass case shows a remarkable similar eccentricity for all particle numbers N.
In addition to this scatter, we observe significantly higher hardening rates for the steeper profiles due to the larger central densities. The mean eccentricities are lower for steeper profiles which compensates partly for the faster evolution of the SMBH binaries due to the delayed influence of GW emission. Compared to the singlemass systems, the hardening rates of the systems with SMF are higher by ~30−40%, whereas there is no systematic effect on the eccentricity observed.
The enhanced hardening rate for the SMF simulations due to mass segregation is also seen in the dimensionless hardening rate H, because the density and velocity dispersion profiles are similar to the nonSMF cases. The reason must be hidden in a more subtle difference inthe phase space distributions. This requires a much more detailed analysis of the system, which must be postponed to a future investigation.
We admit that the impact of the SMF and mass segregation maybe overestimated since the relaxation time in the simulations is too short due to the small particle number compared to realistic systems. Nevertheless, the tendency of speeding up the SMBH binary evolution is interesting because there are other ways to form a mass segregated galactic nucleus, for example by an inhomogeneous mixture of stellar populations with different ages.
We have applied the simulations to three representative galaxies for the three different density profiles. We found coalescence times of 0.30 ± 0.04 Gyr for the MW (representing a steep slope of γ = 1.5), 0.38 ± 0.10 Gyr for M31 with intermediate slope γ = 1.0, and 1.26 ± 0.36 Gyr for M87 with a shallow slope γ = 0.5. In all cases the coalescence time is short compared to the Hubble time and the expected time between two mergers in the presentday Universe. At high redshifts where galaxies are compact and thus possess denser central regions, SMBH merger times can be orders of magnitude smaller (Khan et al. 2016).
We have also calculated the strength of the GWs emitted in the final phase of the SMBH binary evolution and have shown that lowmass and intermediatemass mergers (M_{∙} ~ 4 × 10^{6}− 2 × 10^{8} M_{⊙}) are visible with eLISA at the lowfrequency end at redshift z = 3 (corresponding to a luminosity distance of D = 26 Gpc), whereas highmass mergers (M_{∙}~ 6 × 10^{9} M_{⊙}) would be close to the frequency and sensitivity limit of current PTAs at that redshift.
Acknowledgements
We are thankful to Alberto Sesana for providing sensitivity curve data for various PTAs. F.M.K. acknowledges support by the Excellenzinitiative II “Mobilitätsmaßnahmen im Rahmen internationaler Forschungskooperationen 2015–16” of Heidelberg University. F.M.K. also acknowledges support by the Higher Education Commission of Pakistan through NRPU grant 4159. The special GPU accelerated supercomputer Laohu funded by NAOC/CAS and through the “Qianren” special foreign experts program of China (Silk Road Project), has been used for the simulations. We thank the anonymous referee for the valuable comments which improved the paper significantly.
References
 Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016, Phys. Rev. Lett., 116, 061102 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017, Phys. Rev. Lett., 119, 161101 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 AmaroSeoane, P., Aoudia, S., Babak, S., et al. 2013, GW Notes, 6, 4 [Google Scholar]
 AmaroSeoane, P., Audley, H., Babak, S., et al. 2017, ESA, submitted, [arXiv:1702.00786] [Google Scholar]
 Begelman, M. C., Blandford, R. D., & Rees, M. J. 1980, Nature, 287, 307 [NASA ADS] [CrossRef] [Google Scholar]
 Berczik, P., Merritt, D., & Spurzem, R. 2005, ApJ, 633, 680 [NASA ADS] [CrossRef] [Google Scholar]
 Berczik, P., Merritt, D., Spurzem, R., & Bischof, H.P. 2006, ApJ, 642, L21 [NASA ADS] [CrossRef] [Google Scholar]
 Blanchet, L. 2006, Liv. Rev. Rel., 9, 4 [CrossRef] [Google Scholar]
 Bogdanović, T. 2015, Astrophys. Space Sci. Proc., 40, 103 [NASA ADS] [CrossRef] [Google Scholar]
 Chandrasekhar, S. 1943, ApJ, 97, 255 [NASA ADS] [CrossRef] [Google Scholar]
 Dehnen, W. 1993, MNRAS, 265, 250 [NASA ADS] [CrossRef] [Google Scholar]
 Desvignes, G., Caballero, R. N., Lentati, L., et al. 2016, MNRAS, 458, 3341 [NASA ADS] [CrossRef] [Google Scholar]
 Einstein, A. 1916, Sitzungsberichte der Königlich Preußischen Akademie der Wissenschaften (Berlin: Seite), 688 [Google Scholar]
 Einstein, A. 1918, Sitzungsberichte der Königlich Preußischen Akademie der Wissenschaften (Berlin: Seite), 154 [Google Scholar]
 Graham, M. J., Djorgovski, S. G., Stern, D., et al. 2015, Nature, 518, 74 [NASA ADS] [CrossRef] [Google Scholar]
 Gravitational Observatory Advisory Team 2016, The ESA L3 Gravitational Wave Mission  Final Report, [http://www.cosmos.esa.int/documents/427239/653121/goatfinalrev1.pdf/5b27a84519484c1e9e3867bf7156dfe4] [Google Scholar]
 Gualandris, A., & Merritt, D. 2012, ApJ, 744, 74 [NASA ADS] [CrossRef] [Google Scholar]
 Gualandris, A., Read, J. I., Dehnen, W., & Bortolas, E. 2017, MNRAS, 464, 2301 [NASA ADS] [CrossRef] [Google Scholar]
 Hills, J. G., & Fullerton, L. W. 1980, AJ, 85, 1281 [NASA ADS] [CrossRef] [Google Scholar]
 Hobbs, G., Archibald, A., Arzoumanian, Z., et al. 2010, Class. Quant. Grav., 27, 084013 [NASA ADS] [CrossRef] [Google Scholar]
 HolleyBockelmann, K., & Khan, F. M. 2015, ApJ, 810, 139 [NASA ADS] [CrossRef] [Google Scholar]
 Hulse, R. A., & Taylor, J. H. 1975, ApJ, 195, L51 [NASA ADS] [CrossRef] [Google Scholar]
 Just, A., Khan, F. M., Berczik, P., Ernst, A., & Spurzem, R. 2011, MNRAS, 411, 653 [NASA ADS] [CrossRef] [Google Scholar]
 Khan, F. M., Just, A., & Merritt, D. 2011, ApJ, 732, 89 [NASA ADS] [CrossRef] [Google Scholar]
 Khan, F. M., Preto, M., Berczik, P., et al. 2012, ApJ, 749, 147 [NASA ADS] [CrossRef] [Google Scholar]
 Khan, F. M., HolleyBockelmann, K., Berczik, P., & Just, A. 2013, ApJ, 773, 100 [NASA ADS] [CrossRef] [Google Scholar]
 Khan, F. M., HolleyBockelmann, K., & Berczik, P. 2015, ApJ, 798, 103 [NASA ADS] [CrossRef] [Google Scholar]
 Khan, F. M., Fiacconi, D., Mayer, L., Berczik, P., & Just, A. 2016, ApJ, 828, 73 [NASA ADS] [CrossRef] [Google Scholar]
 Komossa, S., Burwitz, V., Hasinger, G., et al. 2003, ApJ, 582, L15 [NASA ADS] [CrossRef] [Google Scholar]
 Kormendy, J., & Ho, L. C. 2013, ARA&A, 51, 511 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Li, B., HolleyBockelmann, K., & Khan, F. M. 2015, ApJ, 811, 25 [NASA ADS] [CrossRef] [Google Scholar]
 Makino, J., & Funato, Y. 2004, ApJ, 602, 93 [NASA ADS] [CrossRef] [Google Scholar]
 McConnell, N. J., & Ma, C.P. 2013, ApJ, 764, 184 [NASA ADS] [CrossRef] [Google Scholar]
 Merritt, D., & Poon, M. Y. 2004, ApJ, 606, 788 [NASA ADS] [CrossRef] [Google Scholar]
 Mikkola, S., & Valtonen, M. J. 1992, MNRAS, 259, 115 [NASA ADS] [CrossRef] [Google Scholar]
 Peters, P. C., & Mathews, J. 1963, Phys. Rev., 131, 435 [NASA ADS] [CrossRef] [Google Scholar]
 Preto, M., Berentzen, I., Berczik, P., & Spurzem, R. 2011, ApJ, 732, L26 [NASA ADS] [CrossRef] [Google Scholar]
 Quinlan, G. D. 1996, New Ast., 1, 35 [NASA ADS] [CrossRef] [Google Scholar]
 Rantala, A., Pihajoki, P., Johansson, P. H., et al. 2017, ApJ, 840, 53 [NASA ADS] [CrossRef] [Google Scholar]
 Reardon, D. J., Hobbs, G., Coles, W., et al. 2016, MNRAS, 455, 1751 [NASA ADS] [CrossRef] [Google Scholar]
 Sesana, A. 2010, ApJ, 719, 851 [NASA ADS] [CrossRef] [Google Scholar]
 Sesana, A., & Khan, F. M. 2015, MNRAS, 454, L66 [NASA ADS] [CrossRef] [Google Scholar]
 The NANOGrav Collaboration (Arzoumanian, Z., et al.) 2015, ApJ, 813, 65 [NASA ADS] [CrossRef] [Google Scholar]
 Vasiliev, E., Antonini, F., & Merritt, D. 2015, ApJ, 810, 49 [NASA ADS] [CrossRef] [Google Scholar]
 Verbiest, J. P. W., Lentati, L., Hobbs, G., et al. 2016, MNRAS, 458, 1267 [NASA ADS] [CrossRef] [Google Scholar]
All Tables
All Figures
Fig. 1
Shrinking of the separation between the SMBHs in galaxy mergers (top panel for models A, middle panel for models B, and bottom panel for models C). Separation and time are in model units. 

In the text 
Fig. 2
Evolution of the inverse semimajor axis (top panel for models A, middle panel for models B, and bottom panel for models C). The singlemass model 0 is represented with a thicker line in each plot. 

In the text 
Fig. 3
Fractional contribution of each stellar mass bin (M, I, and L for massive, intermediate, and low mass, respectively) to the cumulative mass profiles at an early (T = 150) and late (T = 400) evolution time for representative models A1 (top), B2 (middle), and C1 (bottom). 

In the text 
Fig. 4
Evolution of SMBH binary eccentricities (top panel for models A, middle panel for models B, and bottom panel for models C). Model 0 is represented with a thicker line. 

In the text 
Fig. 5
Histograms showing the number of SMBH binaries in eccentricity bins with width 0.2. 

In the text 
Fig. 6
Comparison of the predicted characteristic strain for models A1, B2, and C1 at redshift z = 3.0 and the sensitivity curves for various GW detectors. 

In the text 
Fig. 7
Evolution of SMBH binary inverse semimajor axis, hardening rates, and eccentricity for our convergence test runs with N = 400k, 1 and 2 million for the singlemass A0, and the SMF case A4. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.