Merging Hierarchical Triple Black Hole Systems with Intermediate-mass Black Holes in Population III Star Clusters

Theoretical predictions suggest that very massive stars have the potential to form through multiple collisions and eventually evolve into intermediate-mass black holes (IMBHs) within Population III star clusters embedded in mini dark matter haloes. In this study, we investigate the long-term evolution of Population III star clusters, including models with a primordial binary fraction of $f_{\rm b}=0$ and 1, using the $N$-body simulation code PETAR. We comprehensively examine the phenomenon of hierarchical triple black holes in the clusters, specifically focusing on their merging inner binary black holes (BBHs), with post-Newtonian correction, by using the TSUNAMI code. Our findings suggest a high likelihood of the inner BBHs containing IMBHs with masses on the order of $\mathcal{O}(100)M_{\odot}$, and as a result, their merger rate could be up to $0.1{\rm Gpc}^{-3}{\rm yr}^{-3}$. The orbital eccentricities of some merging inner BBHs oscillate over time periodically, known as the Kozai-Lidov oscillation, due to dynamical perturbations. Detectable merging inner BBHs for mHz GW detectors LISA/TianQin/Taiji concentrate within $z<3$. More distant sources would be detectable for CE/ET/LIGO/KAGRA/DECIGO, which are sensitive from $\mathcal{O}(0.1)$Hz to $\mathcal{O}(100)$Hz. Furthermore, compared with merging isolated BBHs, merging inner BBHs affected by dynamical perturbations from tertiary BHs tend to have higher eccentricities, with a significant fraction of sources with eccentricities closing to 1 at mHz bands. GW observations would help constrain formation channels of merging BBHs, whether through isolated evolution or dynamical interaction, by examining eccentricities.


INTRODUCTION
A class of black holes (BHs) with masses ranging from 10 2  ⊙ to 10 5  ⊙ , known as intermediate-mass black holes (IMBHs), is speculated to exist in the gap between stellar-mass black holes (SBHs) and massive black holes (MBHs) in the Universe.IMBHs have attracted significant attention due to their potential role in explaining the formation of MBHs in the first hundreds of millions of years after the Big Bang (Volonteri 2010;Wu et al. 2015;Bañados et al. 2018;Greene 2012;Reines & Comastri 2016;Mezcua 2017;Inayoshi et al. 2020;Greene et al. 2020;Wang et al. 2021b), as well as the anomalies observed in dwarf galaxies (DGs), e.g., core-cusp and the number of DGs (Silk 2017;Barai & de Gouveia Dal Pino 2019).In particular, IMBHs could serve as seeds for the formation of MBHs through their coalescence and gas accretion.Moreover, the early feedback ★ E-mail: wanglong8@mail.sysu.edu.cn† E-mail: huyiming@mail.sysu.edu.cnfrom IMBHs could quench star formation (SF), reduce the number of DGs, and impact the central density profile of DGs.
While searching for IMBHs, various possible explanations for their formation have been proposed, which could be divided into three categories roughly.The first one involves the mergers of stellar-mass objects in a dense environment.For example, successive or runaway mergers of SBHs in GCs could produce IMBHs with masses greater than 10 3  ⊙ (e.g.Miller & Hamilton 2002;Gultekin et al. 2004;Giersz et al. 2015;Mapelli et al. 2021Mapelli et al. , 2022;;Rizzuto et al. 2021;Mouri & Taniguchi 2002;Giersz et al. 2015).In environments with higher density, e.g., circumnuclear giant H II region (e.g.Taniguchi et al. 2000) and nuclear star clusters (e.g.Antonini et al. 2019;Fragione & Silk 2020;Kroupa et al. 2020;Fragione et al. 2022a;Rose et al. 2022), massive IMBHs would form, because the gravitational potential from galaxies could prevent merger remnants with natal kick from escaping from clusters, where the natal kick could be from supernovae (SNe) and dynamical few-body interactions, as well as recoil kick from asymmetric GW radiation.Secondly, IMBHs could form in gaseous environments.In AGN disks around MBHs, IMBHs could form through mergers of stars due to mass migration and subsequent gas accretion (McKernan et al. 2012(McKernan et al. , 2014)).IMBHs may also be formed through the direct gravitational collapse of metal-poor giant gas without the formation of stars at the galactic center (e.g.Mayer et al. 2010Mayer et al. , 2015)).Lastly, IMBHs could be formed through the evolution of very massive stars (VMSs) with low metallicities.In regions of starburst SF with high central densities, the successive mergers of massive stars could produce VMSs, which would then evolve to IMBHs (e.g.Portegies Zwart & McMillan 2002;Portegies Zwart et al. 2004, 2006;Freitag et al. 2006;Kremer et al. 2020;González et al. 2021), or even collapse to IMBHs directly without pair-instability SNe (Spera & Mapelli 2017).
Population III (PopIII) star clusters are also potential formation sites for VMSs which could evolve to IMBHs (Sakurai et al. 2017a).Most PopIII stars are massive, and they could merge to VMSs or become IMBHs, due to their extremely low metallicities (Stacy et al. 2016;Chon et al. 2021a;Latif et al. 2022).If PopIII clusters are embedded in mini dark matter haloes (e.g.Skinner & Wise 2020), the haloes would prevent them from being disrupted by the galactic potential, allowing them to survive from the early stage of the Universe to the present (Wang et al. 2022).
In this study, we first examine triple systems comprising IMBH or SBH components dynamically formed rather than being primordial in PopIII clusters, depending on the framework proposed by Wang et al. (2022).Nevertheless, it's important to acknowledge certain limitations in the models introduced by Wang et al. (2022).Firstly, these models do not incorporate primordial binaries.Although the characteristics of primordial binaries within PopIII clusters remain uncertain at present, observations of young clusters indicate a prevalence of multiple systems among OB stars (Sana et al. 2012;Duchêne & Kraus 2013;Gaia Collaboration et al. 2023), suggesting that this trend might extend to PopIII stars as well.Furthermore, in the models presented by Wang et al. (2022), binary mergers are handled using the orbital average method (Peters & Mathews 1963;Peters 1964).However, this method may not provide a robust description of their orbital evolution, especially when they exist within multiple systems (Antonini et al. 2014).
In this investigation, we conduct additional simulations, adhering to the initial conditions outlined in (Wang et al. 2022), to enhance the statistical robustness of our results.Simultaneously, we undertake a series of simulations that incorporate primordial binaries.Furthermore, for triple BHs with merging inner BBHs (hereinafter referred to simply as "merging triple BHs") formed in the petar simulation, we evolve their orbits using the direct few-body numerical code including post-Newtonian (PN) correction tsunami (Trani et al. 2019;Trani & Spera 2020), which can accurately simulate the orbital evolution of triple systems.Additionally, we study GWs from merging inner BBHs under perturbation.Note that henceforth, when discussing merging triple BHs in the context of petar simulation, it will be necessary to clarify that the inner BBHs are evolved to merge using the orbital average method, which may not accurately describe their actual merger processes.
The rest of this paper is organized as follows: we present the initial conditions of PopIII clusters, and the -body and few-body methods (codes) in Sec. 2. The population properties of merging triple BHs and GWs emitted by their inner BBHs are investigated in Sec. 3. Finally, we draw the conclusion in Sec. 4. Throughout the paper, geometrical units ( =  = 1) are used, unless otherwise specified, and the standard ΛCDM cosmological model (Ade et al. 2016) is adopted.

METHOD
The complete process to track the long-term evolution of PopIII clusters and to evolve merging triple BHs within them is depicted in Fig. 1, where a schematic diagram of merging triple BHs is illustrated in Fig. 2. It can be divided into two main parts: simulations of PopIII clusters using an integrated approach combining the -body code petar with single and binary population synthesis code bseemp (Tanikawa et al. 2020), as well as simulations of merging triple BHs using the high-accuracy method with the PN correction tsunami based on the results obtained from the cluster simulations.These two parts will be explained in the following two subsections in detail, respectively.

petar
We simulate the evolution of PopIII clusters using the highperformance -body code petar (Wang et al. 2020b), which combines the particle-tree and particle-particle (P 3 T) algorithm (Oshino et al. 2011) with the slow-down algorithmic regularization (SDAR) method (Wang et al. 2020a).The P 3 T component of the code handles long-range gravitational interactions and is implemented within the framework of pentacle (Iwasawa et al. 2017) and the Framework for Developing Particle Simulator (Iwasawa et al. 2016;Iwasawa et al. 2020).Conversely, the SDAR method is employed to manage shortrange interactions, ensuring precise and efficient treatment of binary orbital evolution and close encounters.Additionally, the impact of the galactic potential on PopIII cluster is implemented by petar using Galpy (Bovy 2015), with a detailed description provided in the following subsection.

bseemp
The single and binary stellar evolution in PopIII clusters is simulated by the fast population synthesis method bseemp, which could trace the stellar wind mass loss and BH formation of PopIII stars with the minimum metallicity  = 2 × 10 −10 .The details of star evolution in bseemp and the definition of production of star evolution are explained in (Hurley et al. 2000(Hurley et al. , 2002;;Tanikawa et al. 2020Tanikawa et al. , 2021b;;Wang et al. 2022).For instance, the BHs within pair-instability mass gap, ranging from 60 ⊙ to 121 ⊙ , are referred to as "pair-instability BHs" (PIBHs).The BHs with mass greater than 121 ⊙ and less than 60 ⊙ are known as IMBHs and low-mass BHs (LBHs), respectively.

Initial condition
The initial conditions for the PopIII clusters in this study are derived from the long-term model "NFWden_long_w9_noms_imf1" without primordial binaries introduced in (Wang et al. 2022).We select this model as it is expected to lead to the formation of VMSs, which eventually evolve to IMBHs with masses of up to 10 3  ⊙ .In order to improve statistical accuracy, we carry on 5 times more simulations (168 simulations) for the model "NFWden_long_w9_noms_imf1" than those performed by Wang et al. (2022).Each simulation is performed using different random seeds for generating initial positions, velocities and masses of stars.
During the stellar evolution, the kick velocities of compact objects forming after SNe are also token into account.We adopt the rapid SNe model for remnant formation and material fallback (Fryer et al. 2012), and the pulsation pair-instability SNe (Belczynski et al. 2016).The kick velocity are assumed to follow a Maxwell velocity distribution with a dispersion of 265km/s, based on the measurement of the proper motions of field NSs (Hobbs et al. 2005).The fallback scenario assumes that some BHs have less or zero kick velocity due to mass fallback or failed SNe, resulting in their retention in star clusters.Meanwhile, the kick velocity of NSs formed by electron-capture supernovae are assumed to follow a Maxwell velocity distribution with a dispersion of 3 km/s.
Additionally, we include a new set of models that incorporate the presence of primordial binaries whose initial period, mass ratios, and eccentricities follow Sana's distributions (Sana et al. 2012).The specific key parameters of these models are highlighted as follows.
The initial mass of clusters  clu is set to 10 5  ⊙ , similar to that of Model A, which is a fiducial model of star clusters embedded in mini-haloes models in (Sakurai et al. 2017b).The initial half-mass radius  h = 1pc is a typical value observed in star clusters.The central density model for clusters in (Michie 1963;King 1966) is adopted, where the central concentration is determined by the ratio between the core radius  c and the tidal radius  t of clusters, denoted as .The initial value of  is set to  0 = 9.
We construct the initial mass function (IMF) of stars in PopIII clusters based on the hydrodynamic simulations (Stacy et al. 2016;Chon et al. 2021b;Latif et al. 2022).The IMF follows a single power-law profile: where  = 1,  min = 1 ⊙ , and  max = 150 ⊙ .The IMF of PopIII star clusters would be more top-heavy, favoring the formation of more VMSs and more BHs, than the canonical IMF with the power index of −2.35 (e.g., Kroupa 2001;Chabrier 2003).
As the structure of the dark matter halo where the PopIII clusters embed is unclear to date, its potential is assumed to follow the model in (Navarro et al. 1996), where the virial mass  vir = 4 × 10 7  ⊙ , the virial radius  vir = 280pc and the concentration follows  () =  (0)/(1 + ), with the assumption that clusters evolve from  = 20 to the present, and  (0) = 15.3.The properties of primordial binaries in PopIII clusters are still uncertain, and the observations of the young SF region show that OB stars are mostly in binaries or high-order multiple systems (Sana et al. 2012;Duchêne & Kraus 2013;Gaia Collaboration et al. 2023), which could be applicable to PopIII clusters as well.In this study, we introduce a set of models featuring similar initial conditions but with a primordial binary fraction of  b = 1, with the characteristics of these binaries aligning with observational constraints derived from (Sana et al. 2012).
We employ petar and to evolve 168 PopIII clusters with  b = 0 and 1 for up to 12Gyr, respectively.

tsunami
Binary mergers in petar are handled by the orbital average method described in (Peters & Mathews 1963;Peters 1964), which may formally break down if the binary is orbited by the third object on a highly inclined orbit at a moderate distance (Antonini et al. 2014).Therefore, we use tsunami to evolve orbits of merging triple BHs in the petar simulation.tsunami is a direct few-body numerical code (Trani et al. 2019;Trani & Spera 2020) used to evolve the dynamics of triple systems with high accuracy.It implements Mikkola's algorithmic regularization with a chain structure (Mikkola & Tanikawa 1999a,b), including 1PN and 2PN precession and 2.5PN GW radiation.

Initial condition
The orbital evolution of merging triple BHs formed in PopIII clusters can be traced by petar.The parameters (masses, positions and velocities) of the components of merging triple BHs are recorded by petar at various times between the first and last record times within the integrator time-step where inner BBHs are merging.The first record time refers to the moment when the triple BHs form in the petar simulation.The last record time represents the moment when the pericenters of the merging inner BBHs are ∼ 100 times the Schwarzschild radii of ( 1 +  2 ).We regard the parameters of merging triple BHs at the first and last record times in the petar simulation as the initial conditions of tsunami simulation, respectively.
Due to the dense star environment of PopIII clusters, merging triple BHs might be perturbed through encounters with field objects on timescales shorter than their evolution time, which would alter their orbital properties significantly or even disrupt them.The timescale for collision with field objects is (Binney & Tremaine 1987) where  is the number density of PopIII clusters.The collision time is estimated to be O(1) Myr, depending on the average number density of PopIII clusters ⟨⟩ = O (100) pc −3 (Wang et al. 2022), the average semimajor axis of outer binaries of merging triple BHs ⟨ 2 ⟩ = O (100) AU, the average total mass of them ⟨ 1 +  2 +  3 ⟩ = O (100)  ⊙ , and the velocity dispersion  = O (1)km/s in the petar simulations1 .Furthermore, we also estimate the timescale of dynamical processes which could disrupt merging inner BBHs, e.g., evaporation and flyby encounter , using Eq. ( 42) in (Binney & Tremaine 1987) and Eq. ( 17) in (Winter-Granic et al. 2023).Both of these timescales exceed 10Gyr, which are much longer than  coll , indicating that they would not occur within  coll .
Therefore, we evolve the merging triple BHs at the first and last record times in the petar simulation with tsunami for up to 10 Myr respectively, and only keep the triple BHs whose inner BBHs could merge eventually, i.e., merging triple BHs in the tsunami simulation.We compile statistics on the merger times of inner BBHs in the tsunami simulation and find that over 90% and 95% of them merge within 0.3 Myr and 0.5 Myr, respectively, which are significantly shorter than 10 Myr.This also indicates that setting the evolution time threshold to 10 Myr in the tsunami simulation is reasonable.

RESULT
In this section, we investigate population properties of merging triple BHs and GWs emitted by their inner BBHs within PopIII clusters, along with the impact of  b on them.

Number of merging triple BHs in petar simulation
We investigate the number of merging triple BHs and merging BBHs in the petar simulation, as shown in Fig. 3.In clusters without primordial binaries (  b = 0), all the merging BBHs naturally form through dynamical processes.On average, there are 5 merging BBHs in one cluster, with approximately one fifth of them occurring in merging triple BHs.If all the stars are paired initially, i.e.,  b = 1, the average number of merging BBHs increases to ∼ 100, with the majority formed through the evolution of primordial binaries, and ∼ 6% formed by dynamical captures.The average number of the merging triple BHs also increases to 2, with most of them having The light coral color represents the total merging BBHs, of which ones formed by the evolution of primordial binaries and dynamical captures are denoted by suffixes "-PB" and "-DY", respectively.The burly wood color represents the total merging triple BHs, the symbols representing the formation channels of their inner BBHs the same as those of merging BBHs.The white dots denote median values.The x-axis represents the cases with  b of 0 and 1, respectively.Note that, in order to better present the chart, the scales located above and below the truncation of the y-axis are set to be different.
inner BBHs originating from primordial binaries and a few formed through dynamical captures.

Comparision between merging triple BHs in petar and tsunami simulation
In this subsection, we will compare merging triple BHs in the petar and tsunami simulations.In order to express convenience hereinafter, we refer to the merging triple BHs in the petar simulation as "petar merging triple BHs".It is also important to clarify here that petar merging triple BHs are those whose inner BBHs are evolved to merge using the orbital average method, which descriptions may deviate from their actual merger processes.We use the abbreviation just for the sake of convenience in the following comparisons.The merging triple BHs whose inner BBHs could evolve from the first and last record times of petar records to merge in the tsunami simulation are referred to "tsunami merging triple BHs (the first record time)" and "tsunami merging triple BHs (the last record time)", respectively.
Fig. 4 provides a comparison among the numbers of the petar merging triple BHs ( tbh,p ), the tsunami merging triple BHs (the first record time;  tbh,t0 ), and the tsunami merging triple BHs (the last record time;  tbh,tf ).For models with  b = 0 and 1, we find that  tbh,t0 ≲  tbh,tf ≲  tbh,p .This will be explained in detail as follows.
The petar merging triple BHs used as the evolved objects in the tsunami simulation are those that have been ascertained to undergo mergers in the petar simulation.Consequently, if the tsunami results show divergent orbital evolution for these triple BHs, some of them may not ultimately merge in the tsunami simulation.When tsunami adopts the parameters of merging triple BHs at the last record time from petar as its initial conditions, the inner BBHs would already be in the advanced stages of merging in the tsunami simulation.Consequently, these results lead to  tbh,tf being greater than  tbh,t0 and comparable to  tbh,p .
In the following, we will further study the difference between the petar merging triple BHs and tsunami merging triple BHs (the first record time) providing more evolution information.Hereinafter, we only analyze tsunami merging triple BHs (the first record time) and they are simply referred to "tsunami merging triple BHs".
We investigate the evolution results of petar merging triple BHs in the tsunami simulation, as listed in Table 1, and investigate their dynamical stability, which could lead to the difference between the petar and tsunami simulations.In particular, stable triple systems have a constant semimajor axis of inner binaries on a secular timescale, unstable triple systems, in contrast, would experience chaotic energy exchange, leading to one body escaping from the systems over a short timescale.The stability could be determined by the following criterion (Mardling & Aarseth 2001) where  mut is the angle between inner and outer orbital planes.In the petar simulation, more than 50% merging triple BHs are stable in both cases with  b = 0 and 1.In the tsunami simulation, almost all the stable petar merging triple BHs could also merge.However, more than half of unstable petar merging triple BHs have different evolution results (merge, break, and neither merge nor break) in the tsunami simulation.This is an inevitable consequence of the chaotic nature of triple systems.Small differences in the initial conditions or in the integration algorithm are bound to create differences in the evolution (Hayashi et al. 2022;Portegies Zwart et al. 2023).This is especially true for meta-stable or unstable triples, where the extreme chaotic behavior can amplify small differences between two neighboring solutions, until their macroscopic outcomes diverge (Lalande & Trani 2022;Portegies Zwart et al. 2022;Trani et al. 2024).Fortunately, it has been shown that, despite different final outcomes for an individual simulation can arise from small differences in accuracy, algorithms, or machine architectures, the statistical outcome of many realizations is independent of these factors (Suto 1991;Boekholt et al. 2020).
We compare orbital evolutions of several representative events in both petar and tsunami simulations.Two stable merging triple BHs are shown in Fig. 5.In the first example (a), there are no discernible perturbations from the tertiary BHs.In the second example (b), significant perturbations from the tertiary BHs affect both the inner and outer BBHs.These two examples illustrate scenarios where either GW radiation or dynamical perturbations predominantly influence the orbital evolution of the inner binaries.To illustrate this distinction, we can employ an analytical criterion proposed by (Antonini et al. 2014): where ℓ 1 = √ 1 −  1 and ℓ GW are the defined dimensionless angular momenta of the inner BBHs and the critical angular momentum below which GW energy loss dominates its evolution, respectively.
For the example (a), ℓ 1 < ℓ GW , the inner BBH decouples from the third BH and evolves as an isolated binary approximately.The orbital evolutions using these two codes consistently agree with each other, implying that the GW effects from the orbital averaged method in petar and the direct numerical integration including the PN correction method in tsunami are practically identical for (almost) isolated binaries.
However, in the case of example (b), where ℓ 1 > ℓ GW , while the orbital evolution using both codes appears similar in the early stages, they exhibit significant differences as time progresses.These differences between the results suggest that when the dynamical influence from the third BHs on the inner BBHs cannot be ignored, the numerical three-body integration including the PN correction method and the orbital-averaged method may not agree with each other anymore, which is consistent with the conclusion from (Antonini et al. 2014).
An unstable triple BH whose inner BBH could merge in the petar simulation but break in the tsunami simulation is also shown in Fig. 6.In the early stage of evolution, the evolutions between the petar and tsunami simulations are nearly identical.However, with the evolution, the differences between these two simulations become more and more significant.In particular, at about 50 yr, when the inner BBH is near the periapsis of outer BBH, the inner BBH starts to merge in the petar simulation.In contrast, the inner BBH gradually escapes from the gravitational influence of the third BH and the outer BBH breaks apart eventually in the tsunami simulation.This divergence in outcomes is a result of the difference between the two methods and the instability of triple BHs triggered near the outer periapsis.
Triple systems may have interesting Kozai-Lidov (KL) oscillation, i.e.,  1 and  mut oscillate quasi-periodically with time (Kozai 1962;Lidov 1962;Naoz 2016).The orbital evolution of several merging triple BHs with KL oscillation in the tsunami simulation, along with them in the petar simulation, is shown in Fig. 7.The orbital elements, including  1 and  mut are almost the same at early stages of evolution in both petar and tsunami simulations.However, as the merging triple BHs evolve,  1 increases rapidly, leading to the merger of the inner BBHs within a short time in the petar simulation (note that the KL oscillation could also occur in the petar simulation, but due to the quick merger of the inner BBH and low time resolution, it is challenging to observe).On the other hand, in the tsunami simulation, the merging triple BH exhibits quasi-periodic oscillations of  1 and  imut over a very long period.For example, the last merging triple BH in the plot evolves about 7×10 4 years in total, with a period of oscillation for  1 being ∼ 7000 years, which is consistent with the KL timescale analytically obtained by (Antognini 2015;Trani et al. 2022) where is the Kepler period of the inner BBHs.Furthermore, we explore the effect of turning off the PN correction in the tsunami simulation and observe the orbital evolution of these events.While  1 and  mut undergo quasi-periodical oscillations, their behavior is distinct.Specially, at certain parameter regions,  1 in the tsunami simulation without PN correction can be either smaller or larger than those in the tsunami simulation, which consists of the prediction that the PN effect can suppress or excite  1 (Naoz et al. 2013).
Based on the above comparison, tsunami (direct three-body numerical integration including the PN correction) provides a more realistic evolution of merging triple BHs.Therefore, in the following  analysis, we focus on tsunami merging triple BHs, unless explicitly stated otherwise.

Population property of merging triple BH
We investigate the mass distribution of merging triple BHs and compare it with that of merging isolated BBHs, as shown in Fig. 8.In the case of  b = 0, both merging inner BBHs and merging isolated BBHs have primary masses  1 of O (100)  ⊙ , which are indicative of IMBHs.Their inner BBHs and isolated BBHs tend to be incomparable, meaning the majority have mass ratio  1 =  2 / 1 smaller than 0.2.In contrast, in clusters with primordial binaries (  b = 1), most primary masses  1 are typically O (10) ⊙ , and component masses tend to be comparable.We can understand the difference as follows.
In the model with  b = 0,  h and  c are smaller, as shown in Fig. 9, due to the absence of binary heating (Wang et al. 2021a).This lead to an environment with higher number density of objects.Therefore, compared to the case with  b = 1, both the number and proportion of VMSs heavier than 400 ⊙ via multiple collisions, leading to the potential evolution into heavier BHs (Wang et al. 2022), become greater in the case with  b = 0, as shown in Fig. 10.In clusters with primordial binaries, most merging inner BBHs and merging isolated BBHs originate from primordial binaries directly, as shown in Fig. 3, which are expected to evolve into lighter BHs.In addition, BHs formed from the evolution of primordial binaries inherit the property of comparable mass, which results in a higher occurrence of mass ratios close to 1.
The mass distribution of merging BBHs inferred by LVKC based on GWTC-3 is also presented2 .Compared to merging inner BBHs and merging isolated BBHs in PopIII clusters, they appear to be much lighter, which could be attributed to the observation selection effect of LIGO/Virgo.Since LIGO/Virgo are more sensitive at high frequency bands (hundreds of hertz), only merging comparable light BBHs with low redshifts are expected to be observed.Therefore, merging BBHs inferred by LVKC are biased towards lower masses.Regarding the outer BBHs of merging triple BHs, most of tertiary BHs  3 are O (10)  ⊙ in both cases with  b = 0 and 1.However, the outer BBHs prefer to be comparable in the case with  b = 1, because primordial binaries with comparable masses dominate the formation of merging triple BHs as well.
The distribution of orbital eccentricity and semimajor axis of merging triple BHs and merging isolated BBHs is shown in Fig. 11.In both of cases with  b = 0 and 1, due to the interaction with the tertiary BHs, merging inner BBHs tend to have significantly higher eccentricities, some of which could be as low as (1 −  1 ) < 10 −3 , even could reach (1 −  1 ) ∼ 10 −5 .In contrast, merging isolated BBHs tend to have steeper distributions of (1 − ), indicating that the proportion of sources with higher eccentricities is smaller than that of merging inner BBHs, as they lack the perturbation from the third BHs.Furthermore, since all the merging isolated BBHs are formed through dynamical captures in the case with  b = 0, the proportion of more eccentric sources is higher than that in the case with  b = 1, where all the sources are formed from the evolution of primordial binaries, inherently tending to have small eccentricities.In both cases with  b = 0 and 1, the outer BBHs of merging triple BHs prefer to have eccentricities 0.6 <  2 < 0.9, due to the gravitational capture process between the inner BBHs and the third BHs.The semimajor axes of the outer orbits tend to be much larger than those of the inner orbits, otherwise the latter would be disrupted by the tidal interaction.In addition, most merging triple BHs have large  mut ranging from 45 • to 135 • , as shown in Fig. 12.Note that in the above comparison, the orbital parameters of merging isolated BBHs are considered at their formation time, which is when BBHs form during the evolution of binary stars.However, for the merging triple BHs, their orbital parameters are considered at the first record time, as the formation time of triple systems can not be defined well.Nevertheless, this comparison can still highlights qualitatively the distinctions in orbital parameters between merging inner BBHs and isolated BBHs.A more fair comparison of orbital parameters at the same frequency will be further explored in the following subsection.
We also study the evolutionary dominance of inner BBHs of merging triple BHs at their first record time, as listed in Table 2.In the case with  b = 0, the orbital evolution of more than 50% of inner BBHs are dominated by the GW driven, while the evolutionary dominance is the dynamical interaction between them and the tertiary BHs in the case with  b = 1.This discrepancy can be explained as follows.The inner BBH in the case with  b = 0 are more eccentric than those  ) in the tsunami simulation, are shown from the upper to lower panels, respectively.Each pair of adjacent panels corresponds to one event.The results in the petar and tsunami simulations are plotted in red dots and orange solid lines, respectively.Additionally, the result in the tsunami simulation without the PN correction is plotted in a gray dashed line (due to the lack of PN correction, the evolution time becomes much longer.For the convenience of comparison, a part of the evolution is truncated).The x-axes are shifted based on the time recorded.in the case with  b = 1, as shown in Fig. 11.This makes that ℓ 1 is more likely to be smaller than ℓ GW in the case with  b = 0.

Merger rate of inner BH-IMBH of merging triple BH
The merger remnants and merger times of inner BBHs of merging triple BHs are shown in Fig. 13.In both cases of  b = 0 and 1, the inner BBHs begin to merge at about 10Myr .The total number of mergers increases logarithmically with time approximately, following the standard distribution of merger time ∝  −1 (e.g.Tremou et al. 2018) until the present Universe, with about 80% occurring before 3000Myr ( > 2).However, the values of  b affect the mass spectrum of merger remnants significantly.Specifically, when  b = 0, most remnants would be IMBHs with O (100)  ⊙ .As  b increases to 1, the remnants become lighter, with most of them being O (10)  ⊙ .We can understand the difference depending on Fig. 8 and the explanation on it.
For comparison, the merging isolated BBHs in cases with  b = 0 and 1 are also investigated.Their mergers follow the trend of inner BBHs of merging triple BHs roughly, but the number of them increases relatively slowly during the period of 200-5000Myr.This is because the isolated BBHs, without perturbation from other BHs, take more time to undergo mergers.For the influence of  b on the remnant masses of merging isolated BBHs, it exhibits a similar trend to that of merging inner BBHs.
The merging inner and merging isolated BBHs can be categorized based on their component masses.The average number  of such mergers and their corresponding merger rates R averaged over redshift are listed in Table 3, where R can be estimated as below where  clu = 10 5  ⊙ is the total stellar mass of PopIII clusters.The stellar mass formed per volume  ★ during a period of time could be determined using a fitting formula based on star formation rate (SFR) density (Tanikawa et al. 2022;Skinner & Wise 2020).The lower and upper limits of  ★ , derived by integrating the formula over the period from 100 to 500Myr (corresponding to a redshift range from 20 to 10) and considering re-ionization history, are 3.2 × 10 4  ⊙ Mpc −3 inner LBBH are 0.01Gpc −3 yr −1 .As  b increases to 1, the merger rates increase by several to dozens of times, as expected, except for that of inner IMBH-BH.The slight decline in the merger rate of inner IMBH-BH could be attributed to the fact that most inner BBHs consist of comparable components with dozens of solar masses in the case with  b = 1, as shown on the upper panel in Fig. 8.For merging isolated BBHs with different components, the upper merger rates are several to dozens of times those of corresponding inner BBHs in the case with the same  b .Depending on these above results, the merger rate of inner and isolated IMBBH could make a contribution to that of IMBBH constrained by LIGO/Virgo/KAGRA collaboration (LVKC) (Abbott et al. 2022(Abbott et al. , 2019)).If we regard PIBH-LBH and LBBH as GW190521-like sources and stellar-mass binary BHs (SBBHs) respectively, their merger rates could contribute to or explain those inferred by LVKC (Abbott et al. 2020a;Abbott et al. 2020b;Abbott et al. 2023).

GW from inner IMBH-BH of merging triple BH
Merging eccentric BBHs are accompanied by the emission of GWs consisting of different harmonics, whereas only the second one is included in the case of circular orbits.The characteristic amplitude of the  order (th) harmonic of inspiral GWs is given, following (Kremer et al. 2019), by where   is the frequency of the th harmonic, M and  are the chirp mass and distance of the merging BBHs, respectively.Among these harmonics, the ( peak =  peak /  orb )th order one focused in the following has the maximum radiation power, where  orb is the orbital frequency calculated by Kepler's third law.The peak frequency  peak can be obtained by (Hamers 2021) where  and  are the semimajor axes and orbital eccentricities of merging BBHs, respectively.As the merging BBHs evolve in their orbits, they undergo successive stages of merger and ringdown after the inspiral stages, rendering Eq. ( 8) unsuitable.Due to the circularization of GW radiation, the orbital eccentricities of BBHs decrease and approach zero as they enter the merger and ringdown phases, making the ( peak = 2)th harmonic become dominant.The GWs from the merger and ringdown phases could be described by PhenomD waveform (Husa et al. 2016;Khan et al. 2016).Considering the significant effect of cosmic expansion on the motion of distant merging BBHs, some physical quantities should be replaced with their redshifted counterparts:  →   (luminosity distance), M → (1 + )M,   →   /(1 + ), in both Eq. ( 8) and PhenomD waveform.Furthermore, the inspiral phase of merging BBHs can last much longer than the observation time of GW detectors, thus the fraction of mission lifetime sources spent within given frequency bins should be considered.Therefore, when calculating the characteristic strain of inspiral GWs, it should be multiplied by the square root of min[1,   ( obs /   )] (D'Orazio & Samsing 2018; Willems et al. 2007;Sesana et al. 2005), where  obs = 5 years is adopted as a fiducial value here.
We calculate the  peak and corresponding GWs from inner BBHs of merging triple BHs, and compare them with those from merging isolated BBHs depending on Eq. ( 8) and PhenomD waveform model.It should be noted that since calculating the orbital eccentricities and semimajor axes of merging inner BBHs with tsunami will be not accurate anymore, when they begin to enter the bands of ground-based GW detectors like LIGO, because in this regime the PN approximation stops holding.Thus, we stop calculating them once  reaches O (10 −4 )AU.At this stage, GW radiation starts to dominate the evolution of merging inner BBHs and the perturbation from the tertiary BHs can be neglected safely, so we continue evolving the orbital of merging inner BBHs until  approaches zero alternatively by the averaged orbital Eqs.(10), which describes the evolution of merging isolated BBHs (Peters & Mathews 1963;Peters 1964) where  = 64 1  2 (1 + 2)/5.As for the merging isolated BBHs, since they are not perturbated throughout their whole evolution, we evolve them with Eqs.(10) until  approaches to zero.The evolution of  1 with  peak of all the inner BBHs, including inner IMBH-BH, within merging triple BHs, is shown on the upper panel in Fig. 14.In the cases with both  b = 0 and 1,  1 of some sources always decreases, similar to those shown on the upper panel in Fig. 10 in (Wang et al. 2022)  4 .This occurs because these sources decoupled from the third BHs, as seen in the example (a) in Fig. 5, and their orbital evolution is dominated by GW radiation.However, the evolution of  1 of the remaining sources displays peaks, oscillations and sharp turning points at lower frequencies before they evolve as isolated binary systems approximately.This is attributed to the oscillations of their orbital elements caused by the dynamical interactions between them and the third BHs before their GW radiation becomes dominant, as observed in the example (b) in Fig. 5 and other sources in Fig. 7. Most merging inner BBHs have detectable orbital eccentricities at the mHz space-borne detectors LISA/TianQin/Taĳi sensitive bands.The distribution of  1 at  peak = 0.01Hz where these mHz detectors are sensitive is plotted on the upper panel in Fig. 15.Most of  1 are larger than 0.01, and a fraction of sources have  1 closing to 1.However, when the merging BBHs are at ∼ 10Hz where the ground-based CE/ET/LIGO/KAGRA become sensitive, they can be regarded as quasicircular due to the circularization of GW radiation, with most  1 smaller than 10 −3 , as shown on the lower panel in Fig. 15.
The characteristic strains of GWs from merging inner BBHs and those of different detector noise, with the color bar denoting their redshifts, are shown on the middle panel in Fig. 14.In both cases with  b = 0 and 1, the characteristic strains of GWs emitted by the merging inner BBHs decoupled from the third BHs are similar to those on the medium panel in Fig. 10 in (Wang et al. 2022) 5 .For the merging inner BBHs that are affected by the tertiary BHs, the characteristic strains of GWs from them also display peaks, oscillations and sharp turning points at lower frequency ranges, before they are decoupled and emit GWs as approximately isolated merging BBHs.Depending on the characteristic strains of GWs and the noise of detectors, certain sources, including merging inner IMBH-BH, are detectable, i.e., the GW strains are higher than those of detector noise within certain frequency ranges.The redshift distributions of detectable sources are shown in Fig. 16.In both cases of  b = 0 and 1, for the space-borne the characteristic strains of GWs with  peak , with the color bar representing the eccentricities at  peak = 0.01Hz of sources.The size of dots scales with the total mass of the merging inner BBHs, the characteristic strains of noise of GW detectors 5 (Kawamura et al. 2011;Robson et al. 2019;Wang et al. 2019;Michimura et al. 2020;Abbott et al. 2017;Hild et al. 2011) are plotted in different colors.The columns from left to right are cases with  b = 0 and 1, respectively.
detectors LISA/TianQin/Taĳi/DECIGO, the majority of detectable sources are located at  < 3. Particularly, the number of detectable sources diminishes for Taĳi, LISA and TianQin sequentially at high redshifts due to their reduced sensitivities at low frequency ranges.
Compared with these detectors, the Deci-Hz space-borne detector DECIGO would have the most detectable sources, as it is sensitive over wider redshift range.For the detectors CE/ET/LIGO/KAGRA, the detectable sources also tend to concentrate at  < 3. CE/ET would observe more sources at higher redshift than LIGO/KAGRA, because of their enhanced sensitivities at lower frequency ranges.The detectable sources in the case of  b = 1 are more than those in the case of  b = 0, as there would be more triple BHs formed when primordial binaries exist, as illustrated in Fig. 3.The characteristic strains of GWs from merging inner BBHs, with the color bar representing their eccentricities  1 at  peak = 0.01Hz, are shown on the lower panel in Fig. 14.In both cases of  peak = 0 and 1, the eccentricities at  peak = 0.01Hz of detectable sources for LISA/TianQin/Taĳi/DECIGO concentrate between O (10 −3 ) and O (10 −1 ), along with a significant fraction of detectable sources with  1 ∼ 1, as shown on the upper panel in Fig. 17.Among the mHz detectors LISA/TianQin/Taĳi, the more sensitive one (e.g., Taĳi) would have more detectable sources, including those with higher eccentricities.This is because detectors with better sensitivities would detect GWs with lower strains caused by higher eccentricities.For CE/ET/LIGO/KAGRA, the  1 at  peak = 10Hz of most detectable sources ranges between O (10 −5 ) and O (10 −3 ), as the sources themselves have small eccentricities, as shown on the lower panel in Fig. 15.CE/ET would have more detectable sources than LIGO/KAGRA, due to their better sensitivities at low frequency range.We also explore the sources that are not detectable by mHz detectors but detectable by at least one ground-based detectors (multiband observable sources).Their eccentricity distributions at  peak = 0.01Hz are illustrated in Fig. 18.In both cases of  b = 0 and 1, a significant fraction of multiband observable sources have eccentricities  1 ∼ 1. e 1 (f peak = 0.01Hz), e(f peak = 0.01Hz) This can be explained by the fact that high eccentricities suppress the strains of GWs so that they are lower than those of the mHz detector noise.The eccentricities of such sources at mHz band could not be identified by LISA/TianQin/Taĳi alone, but they would be constrained via multiband observation, e.g., LISA+ET.This could be performed by tracking sources in ground detected catalogs back to data streams from mHz detectors.If the a source is not a significant trigger in LISA/TianQin/Taĳi, it could have large eccentricity at mHz band.Furthermore, the eccentricities at  peak = 0.01Hz of sources detectable by both at least one mHz detector and at least one ground-based detector are also shown in Fig. 18.The majority of them have eccentricities ranging from O (10 −3 ) to O (10 −1 ).There is also a fraction of multiband detectable sources with eccentricities closing to 1, because of their low redshifts or large masses, as illustrated in the middle panel of Fig. 14.

GW from merging isolated BBH
For comparison, we also plot the evolution of  with  peak of all the merging isolated BBHs in both cases with  b = 0 and 1, as well as the corresponding characteristic strain of GWs emitted from them in Fig. 19.Due to the lack of perturbations from the tertiary BHs, there are no peaks, oscillations and sharp turning points observed in the evolution of  and GW strains with  peak .
The distributions of  at 0.01Hz and 10Hz of merging isolated BBHs are also plotted in Fig. 15.The distribution of  in the case with  b = 0 is peaked at values similar to those of merging inner BBHs, along with a larger fraction of smaller eccentricities, because all the merging isolated BBHs are formed by dynamical captures in this case, resulting in larger orbital eccentricities.However,  in the case with  b = 1 is peaked at significantly smaller values.This is attributed to the prevalence of merging isolated BBHs originating from primordial binaries, which tend to exhibit orbits with lower eccentricities.
The redshift  distributions of detectable merging isolated BBHs for various detectors are plotted in Fig. 20.These distributions are similar to those of merging inner BBHs detectable, because of the selection effect of detectors, i.e., the sources with redshifts that make their GW strains larger than detector noise become detectable.
The distributions of  of detectable merging isolated BBHs at  peak = 0.01Hz and 10Hz for various detectors are illustrated in Fig. 21.In comparison with the case of merging inner BBHs, these distributions of  show a shift towards smaller values of .Particularly at  peak = 0.01Hz, there is a notable absence of eccentricities closing to 1.The distributions of  of multiband observable and multiband detectable merging isolated BBHs at  peak = 0.01Hz and 10Hz across various detectors are also plotted in Fig. 22.These distributions of  also indicate a tendency towards smaller eccentricities, with an increased fraction of low eccentricities, compared with the case of merging inner BBHs.This trend above can be attributed to the fact that merging isolated BBHs without perturbation from tertiary BHs tend to have smaller eccentricities, as depicted in Fig. 15.

CONCLUSION
In this study, we track the long-term evolution of PopIII star clusters embedded in mini dark matter haloes, including models with a primordial binary fraction of  b = 0 and 1, using the -body code petar which binary mergers via GW radiation are dealt with the orbital average method.To obtain a more accurate evolution of triple BHs, we utilize the tsunami code with PN effect to evolve the merging triples found in the petar simulation.Then, we compare the dynamical evolution results of the merging triple BHs between these two methods, and investigate the orbital properties and GW radiation of merging triple BHs in detail.
In the petar simulation, when  b = 0, the inner BBHs of all the merging triple BHs are formed by dynamical capture, and the average number of the merging triple BHs in one PopIII cluster is 1.As  b increases to 1, almost all the merging triple BHs have inner BBHs formed by the evolution of primordial binaries, with the average number of the merging triple BHs becoming 2. In both cases with  b = 0 and 1, the tsunami simulation yields a comparable number of merging triple BHs.Specially, the number of stable merging triple BHs, which account for the vast majority of all the merging triple BHs, are almost equal in these two methods.For unstable merging triple BHs, however, the tsunami results are significantly less.This discrepancy can be attributed to the difference between these two methods and the instabilities of merging triple BHs.The mass distribution of merging triple BHs is affected by the values of  b significantly.In the case with  b = 0, most primary components  1 of inner BBHs are IMBHs with O (100) ⊙ , and the inner BBHs prefer to be unequal.The tertiary BHs  3 concentrate at O (10)  ⊙ , and the outer BBHs also tend to be unequal.When  b = 1, however,  1 becomes lighter, most of them are dozens of solar masses, and the inner BBHs tend to be comparable.For the outer BBHs, most  3 are still dozens of solar masses, but the mass ratio tends to be 1.The mass distribution of merging isolated BBHs follows trends similar to that of the inner BBHs.Unlike the mass distribution of merging triple BHs, the distribution of orbital parameters is not dependent on  b significantly.At the early evolutionary time of the merging triple BHs, both inner and outer orbits prefer to be largely eccentric (∼ 0.9), larger than those of the merging isolated BBHs.Furthermore, the orbital evolution of more than half of inner BBHs   is dominated by GW driven at early evolutionary time in the case with  b = 0, whereas the dominance of the orbital evolution of most inner BBHs is the dynamical interaction between them and the third BHs.
In both cases of  b = 0 and 1, the inner BBHs of merging triple BHs could merge from ∼ 10Myr to the present Universe, with about 80% occurring at the redshift of  > 2. In the case with  b = 0, most merger remnants of inner BBHs are IMBHs with hundreds of solar masses.The upper merger rates of inner and isolated BBH could be 0.4Gpc −3 yr −1 and 15.6Gpc −3 yr −1 , respectively.Specially, the upper merger rates of inner IMBH-BH and inner IMBBH are 0.1Gpc −3 yr −1 and 0.01Gpc −3 yr −1 respectively, which are about one-tenth those of the isolated cases.The upper merger rates of inner and isolated IMBBH could make a contribution to that constrained by GW observations.The upper merger rates of inner and isolated PIBH-LBH are 0.02Gpc −3 yr −1 and 0.1Gpc −3 yr −1 respectively, contributing to or explaining that of GW190521 inferred by LVKC.Furthermore, inner and isolated LBBH could also make a contribution to SBBHs detected by GWs.
The inner BBHs of merging triple BHs tend to have significant orbital eccentricities  1 within space-borne detector bands in both cases of  b = 0 and 1. Specially, most of them have  1 between 10 −2 and 10 −1 , with a significant fraction of them having  1 ∼ 1 at  peak = 0.01Hz where the mHz space-borne detectors LISA/TianQin/Taĳi are more sensitive.When the merging inner BBHs reach at  peak = 10Hz where the ground-based detectors CE/ET/LIGO/KAGRA become sensitive, the residual  1 concentrates between ∼ 10 −5 and ∼ 10 −3 .The merging inner BBHs, including merging inner IMBH-BH and PIBH-LBH, detectable for various space-borne and ground-based detectors concentrated at  < 3.
At  peak = 0.01Hz, the eccentricities of most detectable sources for space-borne detectors are at  1 < 0.1, with a significant fraction of them having  1 ∼ 1.Some sources with small eccentricities at  peak = 0.01Hz could be detectable by both at least one mHz space-borne detectors and at least one ground-based detectors (multiband detectable).While some sources with high eccentricities at  peak = 0.01Hz would not be detectable by mHz space-borne detectors but detectable by at least one ground-based detectors, whose eccentricities at mHz could be constrained via multiband observation, e.g, LISA+ET.
The redshift  distributions of merging isolated BBHs detectable by various GW detectors are similar to those of detectable merging inner BBHs, because of the selection effect of detectors.Compared with merging inner BBHs, the distributions of eccentricity  at  peak = 0.01Hz and 10Hz of merging isolated BBHs exhibittails with lower values.In particular, most merging isolated BBHs would have lower eccentricities when  b = 1.Consequently, the distributions of  at  peak = 0.01Hz and 10Hz of detectable merging isolated BBHs, as well as of multiband detectable and multiband observable merging isolated BBHs, shift towards smaller values, with no significant fraction of sources having  ∼ 1.The difference in eccentricities between merging inner BBHs and merging isolated BBHs, along with relevant GW detection results, are expected to help us understand and constrain the formation channels of merging BBHs in PopIII clusters.
In this work, GW detection results of merging inner BBHs and isolated BBHs are studied qualitatively by comparing characteristic strains of GWs and detector noise.Quantitative and more precise estimations through signal-to-noise ratio calculations will be explored in the future.The distributions of eccentricities at  peak = 0.01Hz of multiband observable and multiband detectable merging isolated BBHs, represented by green and orange dashed lines, respectively.The multiband observable sources denote the merging isolated BBHs are not detectable by mHz spaceborne detectors but are detectable by at least one ground-based detectors.The multiband detectable sources are those detectable by both at least one mHz space-borne detectors and at least one ground-based detectors.The upper and lower panels correspond to the cases of  b = 0 and 1, respectively.

Figure 1 .
Figure 1.The flow chart outlines the complete process of evolving merging triple BHs in PopIII clusters, which consists of two main parts: simulations of PopIII clusters and simulations of merging triple BHs, respectively.Rectangles and ellipses represent target physical systems and methods (codes), respectively.

Figure 3 .
Figure 3.The number of merging triple BHs and merging BBHs during the evolution of PopIII clusters simulated by petar, which is denoted by box charts.Each box chart represents an entire set of simulations.The width of boxes represents the number fraction of merging triple BHs or BBHs within the values of the y-axis, similar to a histogram.The light coral color represents the total merging BBHs, of which ones formed by the evolution of primordial binaries and dynamical captures are denoted by suffixes "-PB" and "-DY", respectively.The burly wood color represents the total merging triple BHs, the symbols representing the formation channels of their inner BBHs the same as those of merging BBHs.The white dots denote median values.The x-axis represents the cases with  b of 0 and 1, respectively.Note that, in order to better present the chart, the scales located above and below the truncation of the y-axis are set to be different.

Figure 4 .
Figure 4. Comparison of the number of petar merging triple BHs and merging triple BHs in the tsunami simulation.Each box chart represents a complete set of simulations.The width of boxes indicates the number fraction of merging triple BH within the values on the y-axis, similar to the way a histogram displays data.The purple color represents petar merging triple BHs.The green and orange colors denote tsunami merging triple BHs (the first record time) and tsunami merging triple BHs (the last record time), respectively.The white dots denote median values.The x-axis represents cases with different values of  b .

Figure 5 .
Figure 5.The evolution of orbital parameters of (a) one stable merging triple BH with  1 = 200.1⊙ ,  2 = 53.3⊙ ,  3 = 44.8⊙ and (b) one stable merging triple BH with  1 = 350.6⊙ ,  2 = 49.2⊙ ,  3 = 55.6⊙ in both petar and tsunami simulations.The results from petar and tsunami simulations are plotted as red dots and orange solid lines, respectively.The x-axes are shifted by the time recorded.

Figure 6 .
Figure 6.The evolution of orbital eccentricity (upper panel) and semimajor axis (lower panel) of an unstable triple BH with ( 1 = 56.1⊙ ,  2 = 55.8⊙ ,  3 = 33.7⊙ ) that the inner BBHs could merge in the petar simulation but break apart in the tsunami simulation.The results from petar and tsunami are plotted in red dots and orange solid lines, respectively.The x-axes are shifted by the time recorded.

Table 3 .Figure 8 .
Figure8.Upper: the distribution of the primary mass  1 and mass ratio  2 / 1 .Bottom: the distribution of the tertiary mass  3 and mass ratio  3 /( 1 +  2 ).The blue and orange colors denote merging triple BHs in the cases with  b = 0 and 1, respectively.In the corresponding cases, merging isolated BBHs are represented by red and gray colors, respectively.The black color denotes merging BBHs inferred by LVKC based on GWTC-3.

Figure 9 .Figure 10 .Figure 11 .
Figure 9. Upper: The evolution of  h of PopIII clusters.Lower: The evolution of  c of PopIII clusters.The cases with  b = 0 and 1 are plotted in spring green and deep pink colors, respectively.The light colors represent 168 simulations, while the deep colors are mean values of them.

Figure 12 .
Figure 12.The distribution of angles  mut between the inner and outer planes of merging triple BHs at the first record time.The lower and upper panels display the distribution of  mut and the corresponding normalized cumulative distribution function (CDF).The blue and orange colors correspond to the cases with  b = 0 and 1, respectively.

Figure 13 .
Figure 13.Upper: the merger remnant masses ( f ) of inner BBHs of merging triple BHs and merging isolated BBHs in the cases with  b = 0 and 1, which are plotted in different colors and shapes, respectively.The merging inner BBHs and merging isolated BBHs from primordial binaries are marked with dots.Lower: the normalized CDF of merger times ( merger ) of inner and isolated BBHs in the cases with different values of  b .The upper x-axis denotes the redshift  corresponding to  merger .The merger times of inner BBHs here are equal to the first record time of inner BBHs recorded by petar plus their evolution time to mergers in the tsunami simulation.The merger time of isolated BBHs refer to the points during the evolution of PopIII clusters simulated by petar when the mergers occur.

Figure 14 .
Figure14.Inner BBHs of merging triple BHs.Upper: the evolution of inner orbital eccentricity  1 with  peak , and the color bar denotes the inner mass ratio  2 / 1 .The cyan and gray shadow areas represent the frequency ranges where the space-borne GW detectors LISA/TianQin/Taĳi and the ground-based GW detectors CE/ET/LIGO/KAGRA are sensitive.Middle: the characteristic strains of GWs with  peak , and the color bar represents the redshift of sources.Lower: the characteristic strains of GWs with  peak , with the color bar representing the eccentricities at  peak = 0.01Hz of sources.The size of dots scales with the total mass of the merging inner BBHs, the characteristic strains of noise of GW detectors 5(Kawamura et al. 2011;Robson et al. 2019;Wang et al. 2019;Michimura et al. 2020;Abbott et al. 2017;Hild et al. 2011) are plotted in different colors.The columns from left to right are cases with  b = 0 and 1, respectively.

Figure 15 .
Figure 15.Upper: the distribution of orbital eccentricity  1 () of merging inner (isolated) BBHs at  peak = 0.01Hz.Lower: the distribution of  1 () of merging inner (isolated) BBHs at  peak = 10Hz.The distributions in the cases with  b = 0 and 1 are plotted in different colors, respectively.

Figure 16 .
Figure 16.The redshift distributions of detectable merging inner BBHs are depicted for various GW detectors, distinguished by different colors and line styles.The left and right columns represent cases where  b = 0 and 1, respectively.The upper panels showcase the space-borne detectors LISA/TianQin/Taĳi/DECIGO, while the lower panels feature the ground-based detectors CE/ET/LIGO/KAGRA.Detectable sources are characterized by GW characteristic strains exceeding those of detector noise within a certain frequency range.

Figure 17 .
Figure17.The distributions of orbital eccentricity at  peak = 0.01Hz and 10Hz for detectable merging inner BBHs across various GW detectors are illustrated using different colors and line styles.The left and right columns correspond to cases where  b = 0 and 1, respectively.The upper panels depict the space-borne detectors LISA/TianQin/Taĳi/DECIGO, while the lower panels represent the ground-based detectors CE/ET/LIGO/KAGRA.Detectable sources are those for which the characteristic strains of GWs exceed those of detector noise within a certain frequency range.

Figure 18 .
Figure18.The eccentricity distributions at  peak = 0.01Hz of multiband observable and multiband detectable merging inner BBHs are shown, represented by green and orange dashed lines, respectively.The multiband observable sources refer to merging inner BBHs that are undetectable by mHz space-borne detectors but are detectable by at least one ground-based detectors.The multiband detectable sources are those detectable by both at least one mHz detectors and at least one ground-based detectors.The upper and lower panels represent cases where  b = 0 and 1, respectively.

Figure 19 .Figure 20 .
Figure 19.Merging isolated BBHs.Upper: the evolutions of orbital eccentricity  with  peak , and the color bar denotes mass ratio  2 / 1 .The cyan and gray shadow areas represent the frequency ranges where the space-borne GW detectors LISA/TianQin/Taĳi and the ground-based GW detectors CE/ET/LIGO/KAGRA are sensitive.Middle: the characteristic strains of GWs with  peak , with the color bar representing the redshift of the sources.Lower: the characteristic strains of GWs with  peak , with the color bar representing the eccentricities at  peak = 0.01Hz of sources.The size of the dots is scaled with the total mass of the sources, and the characteristic strains of the noise of GW detectors are plotted in different colors.The columns from left to right represent the cases with  b = 0 and 1, respectively.

Figure 21 .
Figure21.The distribution of orbital eccentricity at  peak = 0.01Hz and 10Hz of detectable merging isolated BBHs for various GW detectors are represented in different colors and line styles.The left and right columns correspond to the cases of  b = 0 and 1, respectively.The upper panels depict the space-borne detectors LISA/TianQin/Taĳi/DECIGO, while the lower panel show the ground-based detectors CE/ET/LIGO/KAGRA.The detectable sources are those whose characteristic GW strains exceed the detector noise strains within a certain frequency range.

PopIII cluster tertiary BH merging inner BBH outer orbit inner orbit m 1 m 3 m 2
Figure 2. A sketch of merging triple BHs in PopIII clusters.The gray ball represents a PopIII cluster, in which small gray dashed circles are merging triple BHs.The enlarged version of a merging triple BH is in a gray large dashed circle, where a merging inner BBH with primary mass  1 and secondary mass  2 ( 2 ≤  1 ) orbits around a tertiary BH with mass  3 .Orange dotted and red dashed circles are inner orbit with semimajor axis  1 and outer orbit with semimajor axis  2 , respectively.

Table 1 .
Evolution results of petar merging triple BHs in the tsunami simulation and their stabilities.Code b Merge (stable) Merge (unstable) Break (stable) Break (unstable) Neither merge nor break (stable) Neither merge nor break (unstable)

Table 2 .
The evolutionary dominance of inner BBHs of merging triple BHs at the first record time.