1. Introduction
During their lifetimes, most massive galaxies undergo various major merger events, either with other individual galaxies or with groups of galaxies. This has been the basis of the hierarchical cosmic structure formation in the
$\Lambda$
CDM cosmology (Peebles and Yu Reference Peebles and Yu1970; White and Rees Reference White and Rees1978; White and Frenk Reference White and Frenk1991). Galaxies gain mass through hierarchical mergers, and most elliptical and spherical galaxies host a central supermassive black hole (hereafter, SMBH). The empirical scaling relations suggest that the host galaxies and their central SMBHs grow together and suggest a co-evolutionary connection between them like the
$M_{BH}$
-
$\sigma$
relation which links the SMBH mass
$M_{BH}$
to the stellar velocity dispersion
$\sigma$
of the galactic bulge (Ferrarese and Merritt Reference Ferrarese and Merritt2000; Gebhardt et al. Reference Gebhardt, Bender and Bower2000). When the galaxy grows,
$\sigma$
increases and hence subsequently the
$M_{BH}$
. The SMBHs undergo stellar dynamical evolution through mergers of galaxies that leads to the formation of a SMBH binary and its eventual coalescence with the probable emission of gravitational waves detectable by current and future observatories (Amaro-Seoane et al. Reference Amaro-Seoane, Audley and Babak2017).
The dynamical evolution of the SMBH is commonly divided into three distinct stages, governed by different physical mechanisms (Begelman et al. Reference Begelman, Blandford and Rees1980). First, in the post merger phase of the galaxy mergers, the SMBHs sink to the center of the galaxy through dynamical friction with other stellar bodies and particles (Chandrasekhar Reference Chandrasekhar1943). The gravitational drag decelerates the SMBHs, causing them to lose orbital energy and angular momentum. This process operates efficiently on kiloparsec to parsec scales and continues until the two SMBHs form a bound binary. Once the SMBHs form a bound binary, further orbital decay is dominated by interactions with the surrounding stellar environment and this is known as the stellar hardening phase. The pair thus forms a hard binary and continue ejecting stellar particles through three-body encounters (as gravitational slingshot), carrying away energy and angular momentum (Quinlan Reference Quinlan1996; Sesana Reference Sesana2010). This process gradually reduces the binary separation, a stage often referred to as “hardening” and causes the binary to shrink. Finally, at sub-parsec scales, the emission of gravitational waves dominates, rapidly driving the binary to coalescence (Peters and Mathews Reference Peters and Mathews1963; Peters Reference Peters1964).
The timescale from the merger of galaxies to the final coalescence of their central SMBHs can vary significantly, depending on the remnant’s matter distribution and the binary’s orbital parameters, ranging from tens of Myr to nearly a Gyr or longer (Khan et al. Reference Khan, Preto and Berczik2012, Reference Khan, Fiacconi and Mayer2016; Colpi Reference Colpi2014). During this evolution, it is possible that a third galaxy joins this galaxy, which already contains a SMBH binary and create a triple SMBH system, or even three (or more) galaxies directly merge simultaneously. However, the long term dynamics of triple SMBH systems still remains poorly understood. In previous work by Koehn et al. (Reference Koehn, Just and Berczik2023), triple SMBH systems, where two subsequent merger takes places within a Gyr, were extracted from the merger tress of SMBHs and it was found that the triple interactions plays a crucial role in the outcome of the dynamical evolution of the SMBHs. In these simulations spherical galaxy systems were used with the initial conditions taken from the ROMULUS25 cosmological simulation (Tremmel et al. Reference Tremmel, Karcher and Governato2017), and it turned out that the two most massive SMBHs merge first. In this work, we extend these models on triple SMBHs by investigating the role of triaxiality in shaping the dynamical evolution of three SMBHs systems to find the common dynamical evolution pattern of the triple system and predict the typical coalescence time scale using high resolution gravitodynamical N-body simulations.
2. Methods and initial conditions
We use the spherical galaxy models of Koehn et al. (Reference Koehn, Just and Berczik2023) in order to create different models of 3-galaxy triaxial initial conditions with dark matter, gas and stars taking the 3 galaxy radial Dehnen profile parameters along with the SMBH initial conditions from Table 2 and Table 3 of their paper. The masses of the SMBHs are
$88.4\times10^7M_\odot , 13.3\times10^7M_\odot\ \mathrm{and}\ 3.6\times10^7M_\odot$
respectively for BH1, BH2 and BH3. All the three models A, B and C in this work (corresponding to model A1, A2 and A3 in Koehn et al. (Reference Koehn, Just and Berczik2023)) differ in the initial galactic density profiles of the stellar components. We then add the triaxial shape taken from Butsky et al. (Reference Butsky, Macciò and Dutton2016) to these galaxies using AGAMA galaxy modeling software (Vasiliev Reference Vasiliev2019). Here the potential of the SMBHs was not taken into account in order to create in a few dynamical time scales of the nucleus in the simulation, a central cusp of the central potential of the SMBH.
We use different high resolution N-body codes for our simulations. For the initial 3-body interaction phase we start with the BONSAI2 Tree code (Bédorf et al. Reference Bédorf, Gaburov and Portegies Zwart2012) with
$\sim 10$
million particles for each of the galaxy for sufficient interaction of the particles during the early evolution of the global galaxy merger phase and the decay of the SMBHs to the center of the merger remnant due to dynamical friction. Once the binary starts to harden and the separation of one black hole pair fall below
$\sim$
100 pc, the risk of close pericenter passages not being properly resolved due to the lack of uniform timestep starts to grow and we get a time resolution problem. Therefore, at this point we shift to the
$\varphi$
-GPU Direct N-body code (Berczik et al. Reference Berczik, Nitadori and Zhong2011) with an individual black timestep scheme for the hardening phase of SMBHs and correctly resolve the 3-body dynamics close to SMBHs until we get a bound system with a hardening binary. In this step we also reduce the total number of particles from
$\sim30\times10^6$
to
$\sim4.8\times10^6$
particles close to the inner 10 kpc of the SMBHs. We also include Post-Newtonian (PN) approximation up to 2.5 PN order when the black holes reach a separation of 1000
$R_s$
(schwarzschild radius) and move towards the relativistic regime.
3. Results
The orbital evolution of the triple SMBH system showed that they form a binary orbits around the most massive BHs. From time to time there are three body encounters of the SMBHs which switch the members of the binary for all our models (see Fig. 1).
Orbital evolution of the three SMBHs where the three BH1,BH2 & BH3 are color-coded as red, blue and green respectively showing the continuous exchange of the SMBHs due to three-body encounters.
For all the three models, initially the BH2 and BH3 approach each other very closely, which then later during the final galactic collision weakens the BH2 slingshot. Hence, the BH2 immediately returns to the galactic center and we get a BH2-BH3 system (blue line) sharing a common halo and orbiting each other, preventing the formation of the BH1-BH3 system (green line). As the separation falls below sub-pc scale we switch to the direct code (marked in white) between 400–500 Myr. Between 500–700 Myr there is a collision of the BH2-BH3 with the first galaxy and we get an exchange and finally a BH1-BH2 binary (red line) forms, as BH3 is sent off to orbit around the galactic center at 100–500 pc. Now, while the BH1-BH2 binary continues to harden, BH3 slowly descends, its orbit circularizes and a stable hierarchical triple system is formed in all the models, where BH3 becomes bound to the inner BH1-BH2 hard binary (the two most massive SMBHs) as shown in Fig. 2.
Black hole orbital parameters across the different simulations. Color denotes to different SMBH pairs. The red lines refer to the orbital parameters around the inner BH1-BH2 binary (two most massive BHs). Background colors mark the codes used during the particular simulation period, light gray for Bonsai2, white for
$\varphi$
-GPU. The left panels show the separation of different SMBH pairs and the right panels show the evolution of the inverse semi-major axis. The black dotted, dashed, and solid lines, mark the influence radius, hardening radius, and the factor 1000 of the Schwarzschild radius (
$R_s$
) of BH1 respectively.
We take the orbital averaged hardening rate by GW emission from the empirical Peters and Mathews (Reference Peters and Mathews1963) formula assuming a constant eccentricity and add the stellar hardening estimated from the simulation and calculate the merging time of the binary. We get a stellar hardening rate between 1-2 kpc
$^{-1}$
Myr
$^{-1}$
for the hard binary and merging time between
$\sim$
3-13 Gyr which is well within the Hubble time. The actual merging rate of the models differs due to the different hardening rate of the binary system, which actually depends on the details of the orbit. Hence, we conclude that triaxiality does not show any strong impact on the interaction and the final outcome of the triple system as compared to SMBHs present in spherical galaxies.
Acknowledgements
The authors gratefully acknowledge the Gauss Center for Supercomputing eV (www.gauss-centre.eu) for funding this project by providing computing time through the John von Neumann Institute for Computing (NIC) on the GCS Supercomputer JUWELS at Jülich Supercomputing Center (JSC). The authors acknowledge support by the High Performance and Cloud Computing Group at the Zentrum für Datenverarbeitung of the University of Tübingen, the state of Baden-Württemberg through bwHPC and the German Research Foundation (DFG) through grant no INST 37/1159-1 FUGG, and the data storage service SDS@hd supported by the Ministry of Science, Research and the Arts Baden-Württemberg (MWK) and the German Research Foundation (DFG) through grant INST 35/1503-1 FUGG. NS acknowledges the support of doctoral scholarship from the German Academic Exchange Service (DAAD) through grant no. 91864951. PB and MS thanks the support from the special program of the Polish Academy of Sciences and the U.S. National Academy of Sciences under the Long-term program to support Ukrainian research teams grant No. PAN.BFB.S.BWZ.329.022.2023.
Supplementary material
To view supplementary material for this article, please visit https://doi.org/10.1017/S1743921325104596
φ
Rs
Open access
