Analysis of Kozai Cycles in Equal-Mass Hierarchical Triple Supermassive Black Hole Mergers in the Presence of a Stellar Cluster

Supermassive black holes (SMBHs) play an important role in galaxy evolution. Binary and triple SMBHs can form after galaxy mergers. A third SMBH may accelerate the SMBH merging process, possibly through the Kozai mechanism. We use N -body simulations to analyze oscillations in the orbital elements of hierarchical triple SMBHs with surrounding star clusters in galaxy centers. We find that SMBH triples spend only a small fraction of time in the hierarchical merger phase (i.e., a binary SMBH with a distant third SMBH perturber). Most of the time, the enclosed stellar mass within the orbits of the innermost or the outermost SMBH is comparable to the SMBH masses, indicating that the influence of the surrounding stellar population cannot be ignored. We search for Eccentric Kozai-Lidov (EKL) oscillations for which (i) the eccentricity of the inner binary and inclination are both oscillate and are anti-phase or in-phase and (ii) the oscillation period is consistent with EKL timescale. We find that EKL oscillations are short-lived and rare: the triple SMBH spends around 3% of its time in this phase over the ensemble of simulations, reaching around 8% in the best-case scenario. This suggests that the role of the EKL mechanism in accelerating the SMBH merger process may have been overestimated in previous studies. We follow-up with three-body simulations, using initial conditions extracted from the simulation, and the result can to some extent repeat the observed EKL-like oscillations. This comparison provides clues about why those EKL oscillations with perturbing stars are short-lived.


INTRODUCTION
1.1 SMBHs and their host galaxies The Milky Way contains an SMBH at its center and there is convincing observational evidence for the existence of SMBHs residing in the nuclei of all massive galaxies (Kormendy & Richstone 1995;Richstone et al. 1998).Various correlations between the central SMBHs and their host galaxy properties have been explored, such as bulge mass and bulge stellar velocity dispersion or dark matter halo mass, indicating that the formation and evolution of SMBHs and host galaxies are linked (see, e.g., Magorrian et al. 1998;Ferrarese & Merritt 2000;Kauffmann & Haehnelt 2000;Ferrarese 2002; see also Kormendy & Ho 2013 for a review).
Within the framework of a cold dark matter cosmology, dark matter structure and galaxies build up via hierarchical merging (e.g., White & Rees 1978;White & Frenk 1991).As a consequence of the frequent merging of dark matter halos (e.g., Fakhouri & Ma 2008;Genel et al. 2009) and their corresponding galaxies (e.g., Guo & White 2008;Boylan-Kolchin, Ma & Quataert 2008), the formation of supermassive black hole binaries (SMBHBs) might be a natural consequence of mergers of two galaxies with pre-existing black holes (Begelman, Blandford & Rees 1980).In particular massive early-type galaxies, hosting the most massive black holes in the Universe might even have experienced a significant number of mergers with other more or less massive galaxies during the late phases of their assembly (e.g., Bell et al. 2006a;De Lucia et al. 2006;Brown et al. 2007;Bell et al. 2006b;Franx et al. 2008;van Dokkum et al. 2010;Hopkins et al. 2010a;Moster, Naab & White 2013;Behroozi, Wechsler & Conroy 2013).These galaxy merger events provide a natural explanation for the structural evolution of the massive early-type galaxy population (e.g., Boylan-Kolchin, Ma & Quataert 2006;Naab, Khochfar & Burkert 2006;Bezanson et al. 2009;Naab, Johansson & Ostriker 2009;Damjanov et al. 2009;van der Wel et al. 2009;Bundy et al. 2009;Robaina et al. 2010;Hopkins et al. 2010b;Oser et al. 2012;Hilz, Naab & Ostriker 2013;Naab et al. 2013); see Naab & Ostriker 2017 for a review.Mergers between massive galaxies occurred regularly in cosmic history.It is plausible to assume that their black holes merged as well.However, the efficiency of the SMBH merger process is actively debated (Milosavljević & Merritt 2003;Khan et al. 2013;Vasiliev, Antonini & Merritt 2014;Gualandris et al. 2017;Berczik et al. 2022).Galaxy mergers can be broadly divided into minor mergers and major mergers, depending on the mass ratio of the galaxies.In a minor merger (mass ratio smaller than 1 : 10), the mass of the satellite galaxy can be so small that dynamical friction becomes inefficient (Johansson, Burkert & Naab 2009;Callegari et al. 2011;Khan et al. 2012a;Van Wassenhove et al. 2014) and the merger timescale might be longer than the age of the Universe (Milosavljević & Merritt 2001;Poon & Merritt 2001;Yu 2002).For major merger events (mass-ratio greater than 1 : 10), gas-rich galaxies merge much more efficiently than gas-poor galaxies, their merger timescale can be relatively short, and their central black holes can sink rapidly to the center of the remnant (Quinlan & Hernquist 1997;Milosavljević & Merritt 2001;Yu 2002;Boylan-Kolchin, Ma & Quataert 2008).For massive and gas-poor early-type galaxy mergers, the merger timescale can be comparable to the age of the Universe.
Their long-lasting merger process can be divided into three evolutionary phases (Begelman, Blandford & Rees 1980).In the first phase, the infalling galaxy centers (and their SMBHs) are slowed down through the dynamical friction exerted by surrounding stars and dark matter (e.g.Barnes 1988), and the SMBHs sink rapidly to the center as a bound binary (Quinlan 1996).Thereafter, slingshot ejection of surrounding stars on centrophilic orbits becomes the dominant process that takes energy and angular momentum away from the SMBHB and thus leads to the shrinking of the separation between the SMBHBs (Peres & Rosen 1962;Bekenstein 1973;Fitchett & Detweiler 1984;Milosavljević & Merritt 2001;Holley-Bockelmann et al. 2002;Merritt & Milosavljević 2005;Berczik et al. 2006;Lu, Yu & Lin 2007;Chen, Liu & Magorrian 2008).Finally, after the SMBHB reaches a separation close enough (e.g., ∼ 0.1 pc for 10 8 M⊙ SMBHB) to trigger significant gravitational wave (GW) emission, the binary quickly loses binding energy and may finally coalesce (Peters 1964;Portegies Zwart & McMillan 2000;Colpi & Dotti 2011).Low-mass black holes grow mostly by gas accretion, but towards higher masses, their host galaxies become more gas-poor, so the direct merging of SMBHs becomes more important.It might even dominate the mass growth with more than 50% for SMBHs beyond 10 9 M⊙ (Malbon et al. 2007;Volonteri & Ciotti 2013;Kulier et al. 2015); see, however, Menou & Haiman (2004) showing that the overall contribution of merging is small.Many black holes are in binary systems, in particular at high redshift (Volonteri, Haardt & Madau 2003).The merging of SMBHs might even support the formation and explain the tightness of relations between black hole mass and galaxy mass properties (Malbon et al. 2007;Hirschmann et al. 2010;Jahnke & Macciò 2011).
Numerical studies of spherically symmetric galaxies suggested slingshot ejection of stars on centrophilic orbits is not sufficient to bring an SMBHB close enough to trigger GW emission within a cosmic time.This so-called 'final parsec problem' (Milosavljević & Merritt 2003) is caused by the depletion of the 'loss cone' (Frank & Rees 1976;Amaro-Seoane & Spurzem 2001;Milosavljević & Merritt 2001), a region in stellar phase space available for slingshot ejection.Further studies suggest that triaxial or axisymmetric galaxies harbor more stars on centrophilic orbits, reducing the SMBH merging time (e.g., Yu 2002;Merritt & Poon 2004;Berczik et al. 2006;Khan, Just & Merritt 2011;Khan et al. 2013;see, however, Vasiliev, Antonini & Merritt 2014, 2015).However, this long period (∼ 0.1−1 Gyr for the coalescence of SMBHBs with mass ∼ 10 8 M⊙; Khan et al. 2012b) would give a chance to the formation of a secondary merger, when the merger timescale is much larger than the interval between two mergers.Since most massive galaxies should have experienced multiple mergers during their lifetime, triple or higher-order SMBH systems may form.When a third galaxy falls towards the ongoing merger of SMBHB, a hierarchical triple system can form (Blaes, Lee & Socrates 2002;Hoffman & Loeb 2007).If we consider the central inner SMBHB as a subsystem with respect to the third SMBH on a wider orbit, we can describe the system as consisting of an inner binary and an outer binary system, in a hierarchical triple system.

Nuclear star clusters
Most if not all galaxies host nuclear star clusters (NSCs) in their center, typically with masses in the range 10 8 ∼ 10 10 M⊙ (cf.for this and the following the review of Neumayer, Seth & Böker 2020).These extremely dense stellar systems are highly luminous sources that stand out above their surroundings.NSCs and massive black holes are found to exist together (Filippenko & Ho 2003;Graham & Spitler 2009), the most well-known example for this is our Milky Way, which is the host of an inactive SMBH with ∼ 4 • 10 6 M⊙ that resides in a NSC with ∼ 10 7 M⊙ (Genzel, Eisenhauer & Gillessen 2010).NSC reside within the central few tens of parsecs of a galaxy; the NSC of our own Galaxy is very well studied, its stellar distribution (Ghez et al. 2003;Genzel et al. 2003;Gillessen et al. 2017), as well as its dynamical properties (Lu et al. 2009;Bartko et al. 2009;Feldmeier-Krause et al. 2015).NSCs have been discovered in many other galaxies throughout the universe (Böker et al. 2002(Böker et al. , 2004;;Kormendy & Kennicutt 2004;Côté et al. 2006), see Neumayer, Seth & Böker 2020 for a review.
The formation of nuclear star clusters is still an area of active research, and there is no universally accepted model for their origin.However, several mechanisms have been proposed to explain their formation.Some NSCs are believed to form in situ, which means that they originate from the material inside their host galaxies (Agarwal & Milosavljević 2011;Fahrion et al. 2021).Another model suggests that the migration of globular clusters toward their galactic centers leads to the formation of an NSC and may also support the growth of an SMBH by intermediate-mass black holes falling in with the clusters (Feldmeier et al. 2014;Navarro et al. 2023).Stone, Küpper & Ostriker 2017 developed an analytical model of the growth of massive black holes in a NSC.Correlations between NSCs and their central SMBHs suggest that their evolution is coupled (Neumayer & Walcher 2012;Naiman et al. 2015).

SMBH triples
The influence brought by a third galaxy on accelerating the merger process of a progenitor SMBHB system is multifaceted.As the SMBH from the infalling galaxy slowly sinks towards the center, it perturbs the progenitor SMBHB repeatedly, exciting the eccentricity of the binary.This may accelerate the merging process by increasing the release of gravitational radiation (Makino & Ebisuzaki 1994).A side effect of the eccentricity growth of the SMBHB is that the cross-section for the slingshot mechanism is slightly altered as well, and therefore more stars are ejected due to the change of their loss cone.Besides the influence on the SMBHB, the gravitational potential of galaxies is changed by the infall of the third galaxy.This produces additional stars on boxy orbits, which refill the loss cone of the progenitor SMBHBs (Naab, Khochfar & Burkert 2006).Also, the presence of additional stars brought in by the third galaxy provides more stars with centrophilic orbits, which are ejected by the progenitor SMBHB, which accelerates the merger process (Roos 1988).Finally, the presence of a tertiary SMBH could also impact the merger process of the binary (Blaes, Lee & Socrates 2002, Xu et al., in prep).As a result, the merger timescale of the progenitor SMBHB can be dramatically reduced.Alternatively, instead of forming interacting triple SMBH systems in a galaxy center, a free-floating SMBH may also emerge when one of the SMBHs in the triple system is ejected out by slingshot interaction (typically the least massive body is ejected; see Valtonen, Mikkola & Heinamaki 1994).A merger with GW emission of two components of the triplet could happen before the final ejection of one component, depending on the properties of the host galaxy centers and the properties of the SMBHs (e.g., mass ratio, separations, and velocities; see Iwasawa, Funato & Makino 2006;Hoffman & Loeb 2007;Iwasawa, Funato & Makino 2008;Amaro-Seoane et al. 2010;Bonetti et al. 2018).
The process described above is complicated, and one natural aspect of simplifying this process is to ignore the other parts of galaxies and focus on the hierarchical SMBH triple only.Even the systems with hierarchical SMBH triple only are complicated, since general relativity (GR) must be considered beside the gravitational triple system.So here we put aside the GR effect first and consider a simple case with hierarchical triple stellar systems instead.In such systems, while the outer binary and inner binary torque each other and exchange angular momentum, their eccentricities experience periodic oscillations over secular timescales under certain conditions.These oscillations usually occur on a secular timescale that is much longer than their orbital periods (Ford, Kozinsky & Rasio 2000).For non-coplanar orbits, the inclination angles also have corresponding oscillations.According to canonical perturbation theory, no oscillation in the semi-major axis is predicted, since the energy exchange averages out during the secular evolution (Heggie 1975).However, Mazeh & Shaham (1979) find that for triple stellar systems, the net effect of tidal dissipation together with eccentricity perturbations may lead to significant decay of the eccentricity and the semi-major axis of the inner binary.
Since SMBHs are massive compact objects, generalrelativistic precession of the inner orbit should be taken into account.When the precession period from general relativity is comparable to that of the Newtonian perturbations, it could either enhance the periodic oscillation generated by the orbital perturbations (Soderhjelm 1984) or it could break down the oscillation (Holman, Touma & Tremaine 1997).Thus, to further simplify the problem, if we ignore general relativity effects, tidal disruption, and the spin of SMBHs, considering only the Newtonian gravitational interactions between the hierarchical triple SMBHs, then the merger process can be simplified as an isolated restricted three-body problem.

The Kozai-Lidov mechanism
When the inner and outer orbital planes are mutually inclined, this system can be considered as a triple system following the Eccentric Kozai-Lidov (EKL) mechanism, which refers to the Standard Kozai-Lidov (SKL) mechanism releasing the test-particle quadrupole (TPQ) approximation and the circular outer orbit constraints (Ford, Kozinsky & Rasio 2000;Katz, Dong & Malhotra 2011;Lithwick & Naoz 2011;Naoz et al. 2011Naoz et al. , 2013)).Von Zeipel, Kozai and Lidov discovered independently that for a hierarchical triple system, with one component of the inner binary harboring negligible mass and the outer binary on a circular orbit, if the orbital plane of the inner binary and outer binary are mutually inclined (39.2 • < i < 140.8 • ), then the orbital eccentricity and inclination of the test particle orbit will oscillate with an opposite phase, the so-called Standard Kozai-Lidov mechanism (von Zeipel 1910;Kozai 1962;Lidov 1962;Ito & Ohtsuka 2019).Furthermore, Harrington (1968) studied the general mass ratio for hierarchical triples and obtained a similar quadrupole approximation.Soderhjelm (1984) derived octupole level equations in the limit of low eccentricities and inclinations and showed that the octupole approximation plays an important role than the quadrupole approximation in this regime.
More detailed studies on the EKL mechanism have been carried out recently and flipping of the inner orbit inclination orientation from prograde to retrograde has been found (see Naoz 2016 for a brief review and Shevchenko 2017 for a longer one).Ivanov, Polnarev & Saha (2005) discussed the changes in angular momentum within an orbit while studying the EKL cycles in tidal disruption of stars near unequal-mass SMBH binaries.Various definitions of the hierarchy have been further explored for which the averaging approximation is valid (Lithwick & Naoz 2011;Antognini et al. 2014;Katz, Dong & Malhotra 2011;Bode & Wegg 2014).Bhaskar et al. (2021) discussed the limit within which the secular approximation is reliable for a mildly hierarchical triple when the outer perturber is in a circular orbit.Hamers et al. (2013); Petrovich (2015); Grishin, Perets & Fragione (2018); Mylläri et al. (2018) and Zhang, Naoz & Will (2023) also explored the stability of the hierarchical triple system and applicability of the double-averaged procedure under different schemes.Antognini et al. (2014) reported failure in simulating the long-term evolution of EKL oscillations using the double-averaging approximation.Similarly, Luo, Katz & Dong (2016) find rapid eccentricity oscillations in addition to secular oscillations, while using direct N -body simulations instead of secular approximations, and thereafter developed corrected equations for the double-averaging approximation, and after the correction it agrees with the N -body simulation result.
For the hierarchical triple SMBH systems with comparable masses and random configurations of their orbits, neither the TPQ approximation nor the circular outer binary orbit constraint is satisfied, but still similar oscillations can be observed (Blaes, Lee & Socrates 2002).
Three-body systems with certain initial conditions can be chaotic systems (e.g.Li et al. 2014;Bhaskar et al. 2021;Hansen & Naoz 2020).Even the tiny differences between secular approximation and direct N -body method on simulating isolated hierarchical triples can lead to different evolutionary results, as well as their oscillation patterns.This implies that it is difficult to study the dynamical evolution of isolated SMBH triples in galaxy centers with surrounding nuclear star clusters or even more extended galactic bulges.However, as what Lidov & Ziglin (1976) called a happy coincidence, this secular problem is integrable at the quadrupolar approximation (Farago & Laskar 2010).Fortunately, statistical studies are proven to be effective (Heinämäki et al. 1999;Stone & Leigh 2019;Manwadkar, Trani & Leigh 2020;Boekholt, Portegies Zwart & Valtonen 2020).In this regime, would those EKL oscillations still be observed as expected in previous studies, or would they become short-lived even no longer exist?Amaro-Seoane et al. (2010) simulated SMBH triplets in galactic centers with a focus on the cumulative eccentricity distribution during the merger of the inner binary and their gravitational wave radiation signal prediction.Although EKL is believed to play a pivotal role during the evolution of hierarchical triple SMBHs in galaxy centers, it is not clear to what extent this process affects the merger probabilities of the inner SMBHBs.Detailed statistics on the properties of these oscillations, for example, how their timescale is comparable to the EKL timescale, how often these oscillations occur, and how long they last, are required.
As in the SKL mechanism, oscillations driven by the EKL mechanism can also repeatedly excite the eccentricity of the inner binary, which may lead to a fast coalescence of the merging SMBHB since gravitational wave emission efficiency is strongly related to the eccentricity of the inner binary (Peters 1964).Furthermore, Zwick et al. (2020) used a revised formula to improve Peter's estimation, showing there is a 1-10 times deviation while using Peter's formula.The effect of the EKL mechanism on accelerating the merger of the inner binaries in isolated hierarchical triple SMBH systems has been studied by Blaes, Lee & Socrates (2002) using the orbital averaging method for three-body integrations, showing that the presence of the third SMBH can shorten the merger timescale of inner SMBH binary by an order of magnitude in more than a half of all cases.
Few active galactic nuclei (AGN) pairs, candidate SMBHBs, with separation below ∼ 1 kpc have been identified so far; see, for example, Komossa, Baker & Liu (2016) for a brief review.However, recent observations of triple massive black hole systems show that triple SMBH systems appear to be more common than previously believed (Deane et al. 2014;Kollatschny et al. 2020;Foord et al. 2021).The final coalescence of an SMBHB is the most energetic event as well as the most powerful GW emitter in the Universe.The detection of GWs from a coalescing SMBHB would be important since it would provide robust evidence of the existence of black holes and examine Einstein's relativity equations in the strong-field limit.It is expected to be measured in the future, for example, by LISA (Hewitson & Consortium 2014;Amaro-Seoane et al. 2017, 2022, 2023), Pulsar Timing Arrays (Jenet et al. 2005;Verbiest et al. 2016), and SKA (Janssen et al. 2015;Square Kilometre Array Cosmology Science Working Group et al. 2020).For SMBHs with a mass over 10 8 M⊙, their GW wavelength will be detectable only by PTAs.
This study aims to analyze the data from the N -body simulation on hierarchical triple SMBHs residing in the center of gas-poor elliptical galaxies, with a focus on the oscillation verification of known mechanisms, such as near-Kozai oscillation and three-body oscillation generated by energy or angular momentum exchanges between the orbits, as well as other oscillation patterns that are not observed in pure three-body interactions but found in those SMBH triples residing in galactic centers.Although the Kozai mechanism is believed to be an important process in many processes related to three-body interaction at SMBH vicinity in previous studies, it is currently difficult to observe directly.This  0.0 0.0 0.0 0.0 0.0 0.0 Table 1.Initial positions (in pc) and velocities (in km/s) for the central triple SMBHs in 6 simulations (models A-F).Subscript 1 denotes the position or velocity for the first SMBH, subscript 2 for the second, and subscript 3 for the third.The initial mass of each SMBH is 10 8 M ⊙ , and the total mass is 10 10 M ⊙ for all 64 000 particles.
work is aimed at analyzing simulated data as a prediction before we can get a direct observation, and get a better understanding of related physical processes.Due to the complications discussed above and in order to study the impact of the surrounding star or star clusters on the SMBHs, we focus mainly on the simplified situation of equal-mass triple SMBH systems embedded in a stellar environment.Thus, we do not include the GR effect and tidal effects in the simulation but we compute the timescale instead.We use period analysis to extract the oscillation information from the dynamical evolution of the detailed N -body simulations of these hierarchical triple SMBHs in galactic centers, especially oscillations related to the semi-major axis and eccentricity of the inner binary.Then we analyze the relatively long-term periodic behavior with respect to their orbital motions and summarise the patterns of the oscillation and their frequencies.This paper is organized as follows.In Section 2 we introduce our simulation and initial parameter selection.In Section 3, different points of aspects on analyzing the evolution process are explored, as well as a detailed discussion on both the properties of different types of oscillations and their influence on the merger process.Isolated comparison three-body simulations with the same initial condition are presented in Section 4. Finally, we summarise and discuss our conclusions in Section 5.

METHOD AND INITIAL CONDITIONS
We carry out direct-summation N -body simulations to determine the dynamical evolution of the SMBH triples.The N -body simulations are initialized with the total particle number of N = 64, 000, including three SMBHs.Each SMBH has a mass of 1% of the total mass.The N⋆ = 63 997 surrounding equal-mass bodies that make up the remaining 97% of the total mass.The (Plummer 1911) model with an n = 5 polytrope is adapted to generate the initial positions and velocities of the stellar system.We adopt N -body units, which are scale-free (Heggie & Mathieu 1986).For our proposed scheme, setting the mass of each SMBH as 10 8 M⊙ and the length scale fitting the center region of massive elliptical galaxies (500 pc for 1 N -body unit), then the total mass of the galactic nucleus is 10 10 M⊙, and the N -body time unit is 1.6938 Myr.All our particles are non-evolving point masses.
The three SMBHs are initialized as binary SMBHs plus a third SMBH (see Figure 1), in which the separation of two SMBHs is about one order of magnitude smaller than the distance to the third SMBH (see Table 1) (Amaro-Seoane et al. 2010).The third SMBH is initially loosely bound or even unbound (e.g. in Model D) to the inner SMBH binary but note that the out-most SMBH will later pulled back by the stellar cluster due to dynamical friction instead of escape.The three black holes are initially placed within the Plummer radius, with velocities below their escape velocities.We distinguish inner and outer binaries by comparing the relative distances of the three SMBHs.The initial positions and initial velocities of the SMBH triple in each simulation are well spread in parameter space in order to represent different initial conditions.
Plummer's model is not a suitable initial density distribution for NSCs (cf.e.g. Frank & Rees 1976; Bahcall & Wolf 1976; Genzel, Eisenhauer & Gillessen 2010); it is chosen here for simplicity and since within a frac-tion of a relaxation time a central density cusp around the innermost SMBHs does form (Preto, Merritt & Spurzem 2004, for a single SMBH); this method has been used also in our previous work (e.g.Panamarev et al. 2019).In case of our system of triple SMBH the cusp forms around the inner pair of SMBH.The stellar density distribution around binary and triple SMBH after a galactic merger is complex, may exhibit rotation and bars and is currently beyond the scope of our work.
Since the merger time scale of dark matter halos and bulges is much shorter than the typical time interval between two consecutive mergers, the merging of these components except SMBHB should be accomplished before the third SMBH approaches the central SMBHB.This also provides a pathway for the inner binary to dynamically evolve, resulting in a low-eccentricity orbit due to the interaction with surrounding matter.Thus, the initial eccentricity of the inner binary should be low, and the initial separation should be close to the hardening radius (Quinlan 1996).Studies on the coalescence of SMBHBs have shown that the separation of the binary would decrease slowly after they reach this third phase, relying on sling-shot ejections of surrounding stars, typically, at several parsec distances for 10 6 M⊙ SMBHBs, or several tens parsec for 10 8 M⊙ SMBHBs, respectively.
In this scheme, the typical initial semi-major axis of the inner binary is of order ten parsecs, while the outer binary is approximately ten times wider.The simulation is terminated once the inner binary reaches a separation of sub-parsec scale, which is the typical distance for a 10 8 M⊙ SMBHB with a highly eccentric orbit to trigger an immediate merger (within 1 Myr) due to GW emission (Khan et al. 2012b).In other words, the total simulation time depends on the efficiency of the merger process.For simplicity, we do not include the general relativistic effect.Instead, we compute the timescale of the GW emission to have a sense of how strongly it may affect the actual motion of SMBHs.
Our simulations are carried out using the direct N -body code NBODY6++, the parallel version of Aarseth's NBODY6 (Spurzem 1999).The code does not use a softening to the gravitational force in order to attain high accuracy.Instead, the Kustaanheimo and Stiefel (KS) regularisation has been employed to close binaries (Kustaanheimo & Stiefel 1965).By carrying out an ensemble of 1000 sets of three-body simulations with widespread initial conditions, Amaro-Seoane et al. ( 2010) drew the conclusion that information about initial conditions will be rapidly erased by chaotic interactions between SMBH triples, leading to a result of similar statistical characteristics for different initial inclination distributions.The influence of initial conditions is erased shortly after the simulations started.In spite of this, different initial conditions with different energy levels would have different evolution with higher or lower fraction time in binding status, as well as different oscillation types and evolutionary timescales before the inner binary starts to merge.To ensure the diversity of our investigation, the initial orbital parameters and the alignment of three black holes are well spread in the parameter space (Table 1).
In order to study the effect brought by the surrounding star clusters, we take snapshots of the SMBH triples and extract their position and velocity information, then run another set of simulations with the SMBH triple only for comparison.This comparison could give us an idea of how the surrounding stars could affect the central SMBH triples.In the target region that we selected to do the comparison simulation, we obtain snapshots at an interval of 0.1 N -body time units.Detailed method descriptions are shown in § 4.

General properties of orbital evolution
To analyze the evolution of the triple SMBHs we adopt Jacobi coordinates.The masses of the inner binary are m1 and m2, respectively, the distant third black hole has a mass of m3, the semi-major axis of the inner binary is ain, with eccentricity ein, and the semi-major axis and eccentricity of the outer binary (i.e., the system containing the third black hole and the center-of-mass of the inner binary) are aout and eout.The mutual inclination angle between the orbital plane of the outer binary and the orbital plane of the inner binary is i.We adopt this notation even when the hierarchy status of the system switches, i.e., when one component of the inner binary is exchanged with the outer SMBH, or even a component is ejected.With the center-of-mass of the inner binary Min, the gravitational constant G, the angular momentum Lin, and the total energy of the inner binary Ein we compute the semi-major axis and the eccentricity (2) Similarly, we have for the outer binary and where Mout, Lout, and Eout, are the center-of-mass, the relative angular momentum of the outer binary, the total energy, and the z-component of relative angular momentum of the outer binary, respectively.The mutual inclination angle between the orbital planes of the inner binary and the outer binary is Figure 2 shows the evolution of the semi-major axis, eccentricity, and the relative inclination of the orbital planes of the inner and outer binary, for simulation D (see Table 1).
At the start of simulation D, the inner binary has a semi-major axis of order 1 kpc, while the outer SMBH is unbound to the inner binary.As a result of interactions and exchanges, the semi-major evolution track is discontinuous.Note that even for the inner binary, the data of semimajor axis evolution are sometimes discontinuous because the SMBHs are sometimes unbound to the binary system  (but still bound to the cluster), and the semi-major axis becomes temporally negative.Although these binaries appear temporally unbound, they do not escape from the system in most cases, since dynamical friction and frequent encounters with surrounding stars would pull them back towards the cluster center (see their distance to the center in panel D of Figure 3).During this journey (trying to escape and pulled back), the escaping SMBH loses part of its kinetic energy and the triple system becomes more stable.The net effect of these processes leads to the inner binary semi-major axis shrinking rapidly to roughly an order of hundred parsecs within the first 20 Myr, while the semi-major axis of the outer binary transforms its status from unbound to stably bound.When the inner binary or the outer binary is un- bound, or the outer binary is too wide to effectively interact with the inner binary, we refer to this as the unstable phase.
In the middle stage, the third black hole approaches the inner binary, and violent encounters start to play an important role.Close encounters then lead to chaotic motions and frequent exchanges of inner and outer components, as well as kicking out of one component temporally, therefore resulting in intermittent unbound of the binaries, especially the outer binary.Occasionally, after some close encounters between the three SMBHs (e.g., at t ≈ 50 Myr in Figure 2), the inner semi-major axis temporarily becomes close to that of the outer semi-major axis, which means components of inner and outer binary occasionally experience exchanges.Both of these unbound states and exchanges can be put down to the disturbance exerted by the galaxy center environment and the instability of the three-body system itself.In order to determine how often these exchanges occur, we redefine the inner and outer binary by comparing their mutual distances and show these switches with horizontal lines at the top of each panel in Figure 3.When these exchanges occur, the triple SMBHs are all relatively close to each other, and the three-body system is temporally unstable.We thus refer to this as the chaotic interacting phase.
In the final stage, both the inner binary and the outer binary become sufficiently hard, and the three SMBHs mostly retain a hierarchical triple structure.We refer to this as the hierarchical merger phase.During this stage, the inner binary approach is closer to each other efficiently under the repeated perturbations of the third SMBH and surrounding stellar components until the separation between the inner SMBH binary reaches the sub-parsec scale and merges due to GW emission.
Close interactions between stellar objects and central SMBHs lead to a number of interesting consequences.First of all, due to a short non-equilibrium stage immediately after a merger, there are full loss cones and surges of tidal disruption events (TDE) (Chen et al. 2009); this has been followed in more detail including stellar evolution and TDE Li et al. (2023).Furthermore, it has been suggested that the "tidal" disruption of hard binaries leads to hypervelocity stars escaping from the nucleus (and another member of the binary being tidally disrupted or absorbed by the SMBH, Zhang, Lu & Yu 2013).In both cases, if compact objects are involved (neutron stars, stellar mass or intermediate-mass black holes) extreme or intermediate mass inspirals (EMRI, IMRI) with interesting gravitational wave emission signature will ensue (Mazzolari et al. 2022;Naoz et al. 2022).In this paper, we do not discuss these effects, because we want to study the "pure" dynamics of the interaction of Kozai or other cycles of the triple SMBH with the stars.We do neither include stellar evolution or compact remnants, nor stellar binaries in our model.It is the subject of our future work -as a first step (Li et al. 2023) we published a study for TDE caused by different stars and used a method of scaling to extrapolate to the particle number of real galactic nuclei (which we cannot yet simulate directly).For EMRI and IMRI relativistic inspirals have already been included in the code (Cho, Minzburg & Spurzem 2024) EMRI and IMRI have received increasing also as possible sources for future space-based gravitational wave instruments (cf.e.g.Amaro-Seoane 2019; Vázquez-Aceves et al. 2022).

Stability of SMBH triples
We define four phases to describe the dynamical status of the triple SMBH.To study the stability of triple systems, Eggleton & Kiseleva (1995) provided an empirical criterion by numerical experiments, in terms of a critical ratio Ymin of the periastron distance of the outer orbit to the apastron distance of the inner orbit.Mardling & Aarseth (2001) improved the empirical criterion by using the resonance condition that the outer orbital angular velocity at periastron  is equal to the inner orbit frequency.The early study on related topics can be traced back to the stability criteria study of Harrington (Harrington 1977).The implemented new criterion yields a more stringent limit, especially for large eout by making use of KS description.Thus, we computed the Mardling parameter Ym as follows to examine the stability of the triple system where C = 2.8 is an empirical constant, and pcrit is the critical value for the periastron separation of the outer binary.This equation is used to identify the escape of one body if the ratio between periastron separation of the outer binary pout and semi-major axis of the inner binary ain is smaller than the critical value given above by Ym (i.e.pout/ain > pcrit/ain).The equation above is only for coplanar prograde motion, and their later studies include an inclination fudge factor Y fac for a wider application.Thus, we adopt in this study.The result is shown in the upper-left panel of Figure 4.The figure shows that the system frequently switches between stable and unstable states, but tends to be generally unstable at the starting stage and become stable at the late stage.When the system is stable according to the criteria above, the triple SMBHs system can be considered as a stable hierarchical triple system.The binaries are considered unbound binaries if their total energy (Ein or Eout) is positive.Using this modified Mardling criterion and the binding criterion, we classify the status of the triple SMBH system into four categories: (i) hierarchical merger phase (ain > 0; aout > 0; and 0 < ain/pout < Y −1 m ); (ii) chaotic interacting phase (ain > 0; aout > 0; and ain/pout > Y −1 m ); (iii) binary plus single phase (ain > 0; and aout < 0); (iv) unbound phase (ain < 0 and aout < 0).
For each simulation, we compute the fraction of time the SMBH triple system spends on each phase, and the results are shown in Table 2. Triples spend only a small fraction of time in the stable hierarchical merger phase, which is quite different from the results of isolated 3-body simulations shown in § 4.
We show the distance of the three SMBHs from the center-of-mass of the entire system as a function of time in Figure 3.Although both the inner and outer binaries are unbound intermittently, the three SMBH black holes are generally bound.Also, we can find that violent encounters and exchanges of members between inner and outer binaries are common.Ejections are not observed.
The center-of-mass of the stellar population (without SMBHs) is initially located at the coordinate origin.However, the center-of-mass of the entire system has a deviation after we add three SMBHs with different initial conditions.The position of center-of-mass of the system, rcom, is defined as and the velocity of center-of-mass where m * represents the mass of each star and mBH represents the mass of SMBHs, and the same for the position and velocity.Define the distance between the innermost SMBH and the position of center-of-mass to be Rin, while the distance between the outermost SMBH and the position of center-of-mass to be Rout.Assuming a sphere with a center at rcom and a radius equal to Rin/Rout, then the total mass inside this sphere is defined as the enclosed stellar mass inside the orbital radius of the innermost/ outermost SMBH.
In order to understand how the surrounding stellar population affects the dynamics of the triple SMBHs, we show the ratio between the enclosed stellar mass inside the orbital radius of the innermost/outermost SMBH and the mass of one SMBH in Figure 5.At the start of simulation D, the enclosed stellar mass for the innermost SMBH is approximately 10 times the SMBH mass.The corresponding mass ratio for the outermost SMBH is roughly 25.Both of these two quantities decrease dramatically in the first stage and the third stage and behave chaotically in the middle stage just as the evolution of their separation (see Section 3.1).Most of the time, the mass of stars inside the orbital radius of the outermost SMBH is comparable to the mass of SMBH, showing that these stars together may play a comparable role as the outermost SMBH.This indicates that the stellar components of these galaxies can not be ignored when we study the dynamics of the SMBHs.For the region inside the magenta box of Figure 2 where the example oscillation is located, the two ratios are rather low (around one percent of the SMBH mass) such that the steller components could not affect the triple effectively.This may be the prerequisite for the EKL mechanism to dominate the evolution of the triple SMBHs.

Periodic oscillations and their properties
During the second and third stages, various periodic oscillations are observed in the evolution of the orbital parameters between the inner and outer binary, which is the focus of this paper.These oscillations typically have a much longer period than the orbital periods of the inner and outer binary, so-called secular oscillations.An example of these secular oscillations is the Standard Kozai-Lidov oscillation as mentioned above.Under the SKL mechanism scheme, the z-component of the angular momentum of the inner binary is conserved, as a result, the eccentricity of the inner binary and the inclination angle oscillate together with a timescale TKozai.Antognini (2015) found a numerical factor to better estimate this timescale based on previous studies:

Periods of ein and i [Myr]
Figure 7. Statistical results within different windows for auto-correlations and cross-correlations, as a function of time.Upper-left panel: the amplitude of the first trough in the correlation function curve of inner eccentricity within each window, oscillations are present once the amplitude is below zero (see the result of the example window in Figure 6).Upper-right panel: the value of the orbital periods at the first trough in the correlation function results.Middle panel: same as the top panel but for the mutual inclination angle.Bottom panel: same as the top panel but for the cross-correlation between the inner eccentricity and the mutual inclination, note that for crosscorrelation we search for the first peak and the oscillations are present when the amplitude is above zero.where Tin and Tout denote the orbital period of inner and outer binary.
After relaxing the test-particle approximation -the high inclination (> 39.2 • ) and the near-circular outer orbit assumption, the oscillation between the eccentricity and the inclination still occurs but becomes different from the SKL mechanism, especially its timescale.In this so-called EKL mechanism, the octuple-level approximation also plays an important role, and can even be stronger than the quadruple level.For the octuple level approximation, a similar estimation of oscillating timescale as TKozai is given by Katz, Dong & Malhotra (2011) and Naoz et al. (2013) as where ǫoct represents the strength of the octupole level relative to the quadruple level of Newtonian Hamiltonian of hierarchical triples and is defined as follows: As is discussed in previous sections, the motion of the equal mass triple SMBHs residing in the center of galaxy centers is complicated.In order to simplify this entangled problem, we leave the GR corrections and octuple level of EKL mechanism for future study.We study only the case of equal-mass SMBHs.Thus, ǫoct = 0, and the octuplelevel oscillations are not included.Also, these simulations do not include GR corrections, instead, we compute the GR timescale.The effect of GW emission may lead to an immediate merger of the inner black holes once the two SMBHs get close enough.Thus, it is necessary to look at the merger timescale as a result of gravitational wave emission.This time scale Tgr is defined by the ratio between semi-major axis ain, and secular evolution rate of semi-major axis ȧin: The expression for the change in the semi-major axis ȧin, is given by Peters (1964) ȧin = − 64 5 where c is the speed of light.The adjusted time scale TP based on Peter's formula given by Zwick et al. ( 2020) is expressed as follows: The merger timescale required for the inner binary to coalesce only by emitting gravitational waves is shown in the bottom-right panel of Figure 4, for different moments in time.
The Kozai-Lidov precession rate ωKL compared to the GR precession rate ωGR can be expressed as (Chen et al. 2011;Bonetti et al. 2016) To search for such periodic oscillations of a dynamical system over a certain time interval in the simulation data, frequency analysis is a powerful technique.Unfortunately, the common method of Fourier transformation analysis is not suitable to study the periodic properties of those orbital parameters discussed above, since oscillations here are relatively short-lived (typically less than ten periodic cycles) and with dramatic attenuation on both amplitude and period of the oscillation.In this situation, the results obtained from Fourier analysis are usually too weak.Thus, we use the auto-correlation function of an orbital element or crosscorrelation between two orbital elements as an approach to identify the oscillations and determine their periods.
We illustrate our method of investigating the periodic oscillations in Figure 6 by showing an example of the oscillating region inside of the magenta box in Figure 2. The eccentricity of the inner binary ein and relative inclination angle i are both in clear oscillations within the region inside the magenta box.This is a secular oscillation according to the orbital timescale of the inner orbit shown in Figure 4 (upperright panel).In Figure 6, the upper panel shows the amplitude of auto-correlations of ein (red) and i (black).The lower panel shows the cross-correlation between ein and i.The oscillation for each quantity is clear, but their strength (represented by the amplitude) is decreasing irregularly and dramatically.The amplitude of the first minimum of the autocorrelation of both ein and i damped to a value around −0.4, while the first maximum of the cross-correlation damped to about 0.6.Note that the cross-correlation result shows ein and i are in anti-phase.The period of the oscillation also decreases over time.During this process, the eccentricity of the inner binary is effectively excited (at the peaks of this oscillation ein is close to unity, see Figure 2).When the orbits of the inner binary or the outer binary temporarily become parabolic or hyperbolic, the semi-major axis cannot be used.But their eccentricity is valid and could still depict the evolution of their angular momentum.The fraction of time that either the inner binary or the outer binary is unbound is shown in the right-most of Figure 8.
Our paper aims to find secular oscillations with periods much longer than their orbital periods.From the upper-right panel of Figure 4 we find that the orbital period of the inner binary ranges from hundreds of years to 0.1 Myr.Then the corresponding secular oscillations should be at least tens or hundreds times this orbital period.On the other hand, from the bottom-left panel of Figure 4, we find that the typical timescale for EKL cycles is from 0.5 Myr to 10 Myr (ignoring the sharp spikes) when the triple is in a hierarchical merger phase.Also, the window size can not be too large since these secular oscillations are typically short-lived.We, therefore, select a window size of about 10 Myr to cut the data and use the correlation function to analyze the oscillation inside each window.The result of an example window(the magenta box in Figure 2) is shown in Figure 6.Then we use this window to go over the whole data, with a step size equal to one-tenth of the window size, to avoid the missing discoveries caused by the selection of the window position.
The result of the correlation analysis on autocorrelation of ein and i is shown in the upper and middle panels of Figure 7, and the results for cross-correlation between ein and i are shown in the lower panel.The left three plots show the amplitude at the first minimum of the autocorrelation function (top and middle panel) or the amplitude at the first maximum of the cross-correlation function (lower panel).Periods of the oscillations are plotted respectively on the right.For most identified oscillations in ein and i, their amplitudes at the first peak or trough spread from 0 to 0.8, while the identified periods spread from 1 Myr to 4 Myr.
In the region corresponding to the magenta box, the periods of ein and i are from about 4 Myr down to 0.8 Myr.The decreasing orbital periods may be attributed to either the irregular evolution of EKL mechanism of the triple as the example shown in Figure 11 or the perturbation of the surrounding stellar population.While the typical timescale computed from Equation ( 10) is from about 0.6 Myr down to 0.1 Myr, which can also be found in the down-left panel of Figure 4.If we define the hardening rate as then the average hardening rate within the magenta box is 0.009 pc −1 Myr −1 , which is even lower than the value of the average hardening rate over the whole simulation 0.016 pc −1 Myr −1 .This result indicates that when we simulate without including general relativity effects, the EKL mechanism could not offer a dominating impact on accelerating the merger process of the inner binary in this example region, even the eccentricity is excited to the value very close to unity repeatedly.This agrees with our current understanding that EKL mechanism itself could not contribute to orbital shrinkage of the inner SMBHB without GR and tides, no matter whether we have a surrounding cluster or not.Here we record the value of the average hardening rate since it would be interesting to compare this result with simulations including GR in future studies.Also, when we compute the averaged hardening rate for all identified EKL oscillations and the averaged hardening rate for all six simulations, we get similar results as the example discussed above.
Given the fact that even for the one with the best example with clear oscillation, the amplitude of the first minimum/maximum is weak due to the irregular shape of the oscillation curve and their decreasing orbital periods.Thus, we set the critical value of the power spectrum of an identified oscillation at the first trough to be 0.1, which is a relatively low standard, in order to include those irregular oscillations.For oscillations with stronger amplitude at the first minimum than this critical value, we count them as identified oscillations.Then based on the number of these identified oscillations in each window, we plot the statistical results in Figure 8.
If we use a higher value as the critical value for the power spectrum of identified oscillation, for example, 0.5, then the result can be quite different.The total number of identified oscillations would be less than half of the current value.The strength of the oscillation for both eccentricity and inclination angle reaches their largest value at around 0.5 for eccentricity and 50 degrees for inclination.Oscillations with very small eccentricity or inclination are rare since the triples are unstable most of the time.Both eccentricity and inclination angle stay with relatively high values except in the late stage when the inner binary is hard enough shortly before they start the merging process.
Few oscillations are found in ain and aout according to Figure 8, and one of the reasons for this is that when the binaries are unbound the semi-major axis is not defined for hyperbolic orbits.The fraction of time that ein, eout, and i are in oscillation is relatively high.The fraction of time that ein and i are correlated is only slightly smaller than the fraction of time that i is in oscillation, showing a good covariation between these two quantities.This is different from the situation of the eccentricity of the outer binary eout -it has the most fraction of time in oscillation but its crosscorrelation with i is lower.When the eccentricity of the inner binary ein, the inclination angle i and the cross-correlation between ein and i are all in oscillation at the same time, we refer it as the identified EKL oscillation.In this example simulation, the SMBH triple is under EKL oscillation for about 8% of the time.The quantity with the most fraction of time in oscillation is eout (43%), and even much greater than the number of oscillations in ein (19%) and i (11%).Unlike for ein, the link between eout and i is relatively weak, but not negligible.
Also, we plot the averaged results for all six simulations similar to Figure 8 in Figure 9.The average fraction of time that the triples are in EKL oscillation is only 3%, which is relatively low.Mean motion resonances (MMR) are found and discussed mainly in planetary systems with high mass ratios to the central body.But when the ratio of the two orbital periods between our inner binary and outer binary are occasionally close to the ratio of two small integers, MMR may also occur.Thus, we try to find these MMR-like oscillations from our simulation data as well.Unfortunately, as is shown in Figures 8 and 9, the cross-correlation results between ain and aout show that no obvious MMR is observed, indicating the chance for forming this type of resonance is small.

COMPARISON THREE-BODY SIMULATIONS
To verify if the identified oscillations above can be reproduced by pure three-body interactions, especially three-body interaction under EKL mechanism, we carry out simulations with the SMBH triples only with the same initial conditions as the above N -body simulation.We take 20 snapshots with intervals of 0.17 Myr, starting from t = 36.8Myr in the part of the simulation shown in the magenta box of Figure 2, and extract the kinematical data for the three SMBHs.Then we carry out three-body simulations using the CHAIN code (Mikkola & Aarseth 1993;Aarseth 2003) to study the evolution of the isolated three SMBHs for a timescale of 80 Myr.With a focus on studying the effect of surrounding stellar components rather than the effect of general relativity -which has been studied by Bonetti et al. (2016), we use CHAIN instead of ARCHAIN in this study.
As an illustration, we show the first simulation at t = 36.8Myr out of the 20 simulations in Figure 10, and the last simulation at t = 40.0Myr in Figure 11.The SMBH triples clearly follow the EKL mechanism -the inner eccentricity and inclination oscillating with the same period at the anti-phase.Note that we have equal-mass triple SMBHs, so the oscillation would not be as perfect as those in SKL, and do not include the octuple-level oscillation.We successfully reproduced the identified EKL oscillation, but they are different from the original oscillation with surrounding stellar populations.For example, at time t = 36.8Myr, the oscillation period obtained from Figure 6 is about 1.8 Myr, while  that of three body simulations obtained from Figure 10 is about 4.5 Myr.The averaged oscillation period we obtain from these twenty simulations is 4.8 Myr.This is slightly longer but generally agrees with the oscillation periods get from Figure 7 when t = 36.8− 40.0 Myr, which are about 4 Myr down to 2 Myr.Note that the oscillation in three body simulations lasts much longer than the original oscillation shown in the magenta box of Figure 2.
Also, the periods vary dramatically or decrease with time instead of keeping the same as in SKL oscillations as is shown in Section 3.3, especially in Figure 11.The reason for this deviation may be attributed to the mildly hierarchical configuration of our SMBH triples with equal mass.Within each three-body simulation, if the maximum period is more than twice the minimum period, we define this set of simulations as an irregular oscillation.The number of irregular oscillations, such as the example shown in Figure 11, is ten out of twenty in total.Some of these irregular oscillations even break down within 80 Myr, and some of them break down in the next hundred Myr.This suggests that, although these three-body systems are in EKL oscillation, they are already on an irregular oscillation and even close to breaking down.The value of pout/ain/Ym in the upper left panel of Figure 4 shows that the triple system is very close to, but below the stability boundary at t = 36.8− 40.0 Myr.The system is therefore unstable, with a long instability timescale.This can also be verified by the variations in ain in Figure 10.If the system were stable, ain would be constant, apart from the slight variations on the orbital period timescale.Nevertheless, when the oscillation starts at t = 36.8Myr, the isolated triple is relatively more stable than at the end stage, at around t = 40.0Myr.This is further supported by the decrease of pout/ain/Ym in Figure 4.When in isolation, the SMBH triple could evolve more than 10 Gyr before their EKL oscillation breaks down, and the amplitude and the duration of each period of the oscillation evolve dramatically slower than the one in the magenta box.
Note that not every identified oscillation in the N -body simulation can be produced by the three-body simulation.In fact, even the best examples of identified oscillation besides the one in the magenta box could hardly be repeated with the three-body simulation using the same initial condition.While the other particles besides SMBHs are removed, the triple SMBHs are sometimes unbound to each other, or bound but not in oscillation, or in oscillation but on a very different timescale (e.g., on a much longer timescale; up to several Gyr).

DISCUSSION AND CONCLUSIONS
We have analyzed the evolution of several pathfinding numerical simulations of triple SMBHs in the center of surrounding star clusters with a mass one hundred times each SMBH mass.Depending on the stability status of SMBH triples, we have defined four characteristic phases and summarised the time fraction the triple systems spend in each phase.
We find that unlike pure 3-body simulations shown in Section 4, the triples spend only a small fraction of their time in the hierarchical merger phase, indicating that the surrounding stellar population plays an important role.This is supported by the evidence that most of the time, the enclosed stellar mass within the orbits of the innermost or the outermost SMBH is comparable to the mass of SMBH.In the example simulation D, at a region inside the magenta box where the enclosed stellar mass is relatively low, an oscillation similar to the EKL oscillation is observed.Although it is one of the long-lived examples observed, it lasts for only several periods, with its amplitude damping and its orbital period shrinking with time.We further confirmed this oscillation to be a disturbed EKL oscillation in two ways.First, we use the correlation function analysis to confirm that (i) ein and i are both in oscillation, (ii) the oscillation found in ein and i are correlated with an anti-phase or in-phase and (iii) the oscillation period is comparable to the typical timescale of EKL oscillation.
We extract the dynamical information of the SMBH triples and use this information as the initial conditions for the comparison of three-body simulations.The averaged oscillation period obtained from these simulations generally agrees with the oscillation periods obtained by the correlation function.We also found that for these equal-mass SMBH triples, EKL oscillations are distorted and even become irregular.We use the correlation function to find oscillations over the entire data and show the statistical results for each type of oscillation.
The overall average fraction of time that the triples are in EKL oscillation for all six simulations is only 3%.In simulation D, the SMBH triple is under EKL oscillation for about 8% of the time.The averaged hardening rate of the inner binary during the EKL oscillations is 0.009 pc −1 Myr −1 , which is even lower than the averaged hardening rate within the full-time simulation.This suggests that the role of EKL oscillations in accelerating the merger process of the inner binary may have been overestimated in previous studies.We also find that the orbital parameter with the largest fraction of time in oscillation is eout, and the link between eout and i is relatively weak, but not negligible.No MMR-like oscillation is observed, indicating the chance of forming this kind of resonance is rather low.
Several physical processes have not been included in our current study, which may be considered in future investigations.Our study is limited to major mergers with equal-mass SMBHs, and our future study on unequal-mass systems may be fruitful.The gravitational interaction of the triple SMBH system is coupled to other processes such as the mass growth of the SMBH, tidal effects, and gravitational wave emission due to the spin of the SMBH, etc.Similarly, Mikkola & Tanikawa (1998) studied the mass transfer in the secular evolution of the triple system CH Cygni, triggered by eccentricity oscillation of the inner binary.Also, tidal interaction together with mass transfer for compact xray binaries in hierarchical triples has been studied (Bailyn 1987).

Figure 1 .
Figure 1.Example of the initial condition of the hierarchical SMBH triple.The two red dots represent the inner binary, and the green dot represents the third outer black hole.The curves show the trajectories of the three SMBHs at the beginning of the simulation, red for the inner binary and green for the third SMBH.Black dots indicate the stellar population (each with a mass of 1.53 × 10 5 M ⊙ if we scale each SMBH mass to 10 8 M ⊙ ) that represents the galactic centre.

Figure 2 .
Figure2.Evolution of the semi-major axis (upper panel), eccentricity (middle panel), and relative inclination between the inner and outer orbits (bottom panel), in simulation model D. The red curve represents orbital parameters of the inner binary and the green curve represents the outer binary.The discontinuities in the semi-major indicate times in which the three-body system is unbound.The magenta box highlights a selected time span during which clear oscillations are observed.

Figure 3 .Figure 4 .
Figure3.Distance between the SMBHs and the enter-of-mass of the galactic nucleus, as a function of time.The green, red, and magenta curves represent the first, second, and third SMBH, respectively.Note that when two of the SMBH form a close binary, their distances are nearly identical and the curves overlap.The horizontal curves at the top of each panel show the switches in the order of the hierarchical triple due to the exchanges of the members of the inner and outer binary, since here (and only for this figure) we define the inner binary as the pair with the smallest separation.

Figure 5 .
Figure5.The ratio between the enclosed stellar mass inside the orbital radius of the innermost SMB (red) and outermost SMBH (green) and the mass of a single SMBH as a function of time.

Figure 6 .
Figure6.Auto-correlation of inner eccentricity (red curve) and relative inclination angle (black curve) and cross-correlation between them.

Figure 8 .Figure 9 .
Figure 8. Fraction of time that the orbital elements are in oscillation for model D. The first five bars are the identification of auto-correlation results.The subsequent four bars are the results of the cross-correlation between the two parameters.The result labeled 'EKL' represents the case in which the cross-correlation between e in and i, and the auto-correlations of e in and i are all in oscillation.The quantity α at the right represents the fraction of time that either the inner binary or the outer binary is unbound.

Figure 10 .
Figure10.Evolution of the orbital elements of an isolated hierarchical triple SMBH within the first t = 80 Myr.The initial conditions for three SMBHs are obtained from the snapshot of model D at t = 36.8Myr at the beginning of the magenta box.Unlike the triple in the clustered environment (shown in Figure2), the EKL oscillation in this isolated system persists much longer than 15 Myr, and the distortion of the oscillation curve is limited.

Figure 11 .
Figure 11.Evolution of the orbital elements of an isolated hierarchical triple SMBH.The initial conditions for three SMBHs are obtained from the snapshot of model D at t = 40.0Myr.