Dynamical Evolution of Hierarchical Triple Systems and their Gravitational Wave Detectability

1. Introduction

The existence of gravitational waves (GWs) was first predicted theoretically by Einstein using his theory of General Relativity (GR) [1]. Unlike electromagnetic radiation, GWs are produced by quadrupole sources rather than dipoles, allowing Einstein to describe the rate of energy loss due to GW emission. GWs arise from perturbations in the spacetime metric caused by asymmetries such as binary systems, deformations in compact objects like neutron stars (NSs), or density fluctuations during inflation. The GW signal is typically classified into three phases: inspiral, merger, and ringdown. During inspiral, a compact binary slowly loses energy via GW emission, increasing in frequency and amplitude. The merger produces the peak GW signal as the objects coalesce, followed by the ringdown, dominated by the relaxation of the remnant and emission of decaying quasinormal modes [2].

The first direct GW detection occurred in 2015 when LIGO observed GW150914 from a binary black hole merger [3], opening a new era of GW astronomy. Ground-based interferometers like LIGO, Virgo, and KAGRA (forming the LVK collaboration) currently operate in the 10 Hz–10 kHz range and have detected over 100 compact object mergers. However, signals below 10 Hz remain inaccessible, motivating larger ground-based detectors such as the 40 km Cosmic Explorer and Einstein Telescope, as well as space-based observatories like the Laser Interferometer Space Antenna (LISA) [4] and the Lunar Gravitational Wave Antenna (LGWA) [5]. These extended sensitivities are crucial for detecting eccentric binaries, whose GW emission is distributed over multiple harmonics and is more efficient than that from circular orbits [6,7]. The total radiation power P averaged over one orbital period is given by

$$\mathcal{P}={\displaystyle \frac{32G^4{m}_1^2{m}_2^2{m}_{\mathrm{in}}}{c^5{a}_{\mathrm{in}}^5{(1-e_{\mathrm{in}}^2)}^{7/2}}}\left(1+{\displaystyle \frac{73}{24}}{e}_{\mathrm{in}}^2+{\displaystyle \frac{37}{96}}{e}_{\mathrm{in}}^4\right)$$

(1.1)

where $a_{in}$ and $e_{in}$ are the semi-major axis and the eccentricity of the inner binary, respectively. It shows that the radiation power increases rapidly as the eccentricity increases. The total radiation power P is the sum of the power radiated in the nth harmonic ${\mathcal{P}}_n$, given as

$$\mathcal{P}=\sum _{n=1}^{\infty }{\mathcal{P}}_n$$

(1.2)

where the ${\mathcal{P}}_n$ is the power due to the nth harmonic. The peak of the spectrum $P_n$ is given by

$$f_{peak}=n_{peak}{f}_{in}$$

(1.3)

where $f_{in}$ is the orbital frequency of the inner binary and $n_{peak}$ is the magnification factor, which represents the number of harmonics. This factor can be approximated, in the range ${10}^{-6}<1-e^2<1$ [6,8], to

$$n_{peak}={\displaystyle \frac{2{(1+e)}^{1.1954}}{{(1-e^2)}^{1.5}}}$$

(1.4)

We can see that if eccentricity $e<<1$, $n_{peak}\sim 2$, describing the peak frequency of GWs from a circular orbit, whereas for high eccentricities, $n_{peak}$ is very high and would be in the detectable range of the detector (${10}^{-4}$Hz to 10kHz). Also, since the signal from a highly eccentric orbit consists many harmonics, the signal has a broad spectrum of frequencies. We also check the strain of the wave at the nth harmonic, given by

$$\langle h_n^2\rangle ={\left(2\pi f_n\right)}^2\langle h_n^2\rangle \approx {\displaystyle \frac{4G}{c^3{D}^2}}{P}_n(a,e)$$

(1.5)

where $f_n$ is the frequency of nth harmonic, D is the distance of the source. At $f_n=f_{peak}$, we can approximate the equation and express in scaled form as

$$\sqrt{\langle h_n^2\rangle }\approx 5\times {10}^{-25}{g}^{1/2}{\displaystyle \frac{(n_{\mathrm{peak}},e)}{n_{\mathrm{peak}}}}\left({\displaystyle \frac{m_1}{M_{\odot }}}\right)\left({\displaystyle \frac{m_2}{M_{\odot }}}\right)\times {\left({\displaystyle \frac{a_{\mathrm{in}}}{\mathrm{AU}}}\right)}^{-1}{\left({\displaystyle \frac{D}{10\mathrm{kpc}}}\right)}^{-1}$$

(1.6)

For $m_1=m_2=10M\odot ,a_{in}=0.01AU,D=10kpc$,

$$\sqrt{\langle h_n^2\rangle }\approx 1.9\times {10}^{-22}$$

(1.7)

which is highly detectable in our present and future detectors. Eccentric binaries are highly efficient GW emitters and hence appear at lower frequencies compared to circular orbits. These kinds of eccentric systems must dominate the detections compared to circular orbit system but due to the circularization of the eccentric orbit by GW emission, by the time of merger, the orbit has almost completely circularized, causing it to appear in the LVK band as a binary in a circular orbit. Some residual eccentricity may remain, but has a very negligible effect on the signal.

The formation of eccentric orbits can result from dynamical interactions with multiple bodies, particularly through the Kozai-Lidov (KL) mechanism in hierarchical triple systems (HTSs). The KL mechanism [9,10], independently identified by Yoshihide Kozai and Mikhail Lidov in the 1960s, occurs when a tertiary companion perturbs an inner binary, causing periodic exchanges between its eccentricity and inclination over timescales much longer than the orbital period. This interplay produces long-term modulations in both the inner binary’s eccentricity and the angular displacement between the inner and outer binaries.

Analytically, the Hamiltonian of a hierarchical triple system can be decomposed into contributions from the inner and outer binaries, along with a term representing their mutual interaction [11]:

$${\mathcal{H}}_{total}={\mathcal{H}}_{inner}+{\mathcal{H}}_{outer}+{\mathcal{H}}_{perturbation}$$

(1.8)

This decomposition allows separate treatment of the inner and outer binary evolutions while capturing the essential perturbative effects that drive KL oscillations.

In the context of gravitational-wave (GW) astrophysics, the KL mechanism is significant because it can maintain high eccentricities ($e>0.1$) for binaries entering the LIGO detection range, which may accelerate merger times. Simulations suggest that a substantial fraction of binaries in dense stellar environments are influenced by third-body perturbations, making accurate modeling of their dynamics essential. This requires integrating third-body interactions along with post-Newtonian relativistic corrections [11], especially when predicting merger rates of eccentric binaries.

A fundamental assumption in this configuration is that the tertiary’s gravitational influence on the inner binary is weaker than the interaction between the two inner bodies. In a simple two-body Newtonian system, a bound orbit is elliptical and characterized by six orbital elements ($a,e,i,\omega ,\mathsf{\Omega },M$). In hierarchical triples, perturbations from the tertiary lead to deviations from these idealized ellipses, affecting both the inner binary’s orbital shape and orientation over time [6].

To account for these perturbations, the concept of an osculating orbit [12] is introduced, providing an approximate trajectory that can be described by an elliptical path defined by the aforementioned six orbital elements, derived from the system’s instantaneous position and velocity at any given moment. As for the outer orbit, this is considered as the motion of the inner binary’s center of mass as it rotates around the tertiary companion, which can also be represented as another osculating orbit. In this framework, the masses of the inner binary components are denoted $m_1$ and $m_2$, while the tertiary companion has a mass $m_3$. To distinguish between the inner and outer orbits, subscripts "in" and "out" are used respectively. In such systems, the orbital periods of the inner and outer binary is given by

$$P_{in}=2\pi \sqrt{{\displaystyle \frac{a_{in}^3}{G{m}_{in}}}},P_{out}=2\pi \sqrt{{\displaystyle \frac{a_{out}^3}{G{m}_{out}}}}$$

(1.9)

At a particular relative inclination I, a secular change of orbital elements may occur where under certain conditions, an oscillation between the values of the eccentricity and relative inclination occurs with the relative equation defined as the argument between the inner and outer orbital planes [9,10] and is given by

$$cosI=cos{i}_{in}cos{i}_{out}+sin{i}_{in}sin{i}_{out}cos\left({\mathsf{\Omega }}_{in}-{\mathsf{\Omega }}_{out}\right)$$

(1.10)

which results in the secular exchange of $e_{in}$ and I with the conserved value [13] of

$$\mathsf{\Theta }=(1-e_{in}^2)co{s}^{2}I$$

(1.11)

with the condition for the KL oscillation, provided that $m_1<<m_2,m_3$ and the outer orbit is circular.

$$0\le \left|cosI\right|\le \sqrt{{\displaystyle \frac{3}{5}}}$$

(1.12)

Note that the KL oscillation will occur even if $m_1$ is not much smaller than $m_2$, $m_3$. The KL oscillation timescale [14] is given by

$$t_{\mathrm{KL}}\sim P_{\mathrm{in}}{\displaystyle \frac{m_{\mathrm{in}}}{m_3}}{\left({\displaystyle \frac{a_{\mathrm{out}}}{a_{\mathrm{in}}}}\right)}^3{\left(1-e_{\mathrm{out}}^2\right)}^{{\displaystyle \frac{3}{2}}}$$

(1.13)

When modeling an HTS system, we must also consider the general relativistic precession in the system and the Lense-Thirring, (LT) effect, which is a general relativistic correction to the precession of a gyroscope in the presence of a large rotating mass. This causes the relative orbital nodes of the gyroscope to precess[15]. Recent studies have shown changes in KL oscillations caused due to GR effects involving a SMBH. If the third body is a rapidly rotating SMBH, then LT effect might become important in changing the evolution of eccentricity excitation from the usual KL oscillation. The effect appears in the 1.5PN order [16,17,18]. For a rotating black hole, the spin angular momentum [6,16,18] is given as $S_3={\chi }_{3}G{m}_3^2/c$, where ${\chi }_3$ is the spin parameter, and the outer orbital angular momentum [16,18] is given as

$$L_{out}={\mu }_{out}\left[G(1-e_{out}^2)a_{out}M\right]$$

(1.14)

where M is the total mass of the gyroscope and the reduced mass of the outer orbit is defined by

$${\mu }_{out}={\displaystyle \frac{m_{in}{m}_3}{m_{out}}}$$

(1.15)

The timescale of the orbit-averaged precession of $L_{out}$ around $S_3$ [16,18] is given by

$$t_{L_{\mathrm{out}}{S}_3}={\displaystyle \frac{2c^2{(1-e_{\mathrm{out}}^2)}^{\frac{3}{2}}{a}_{\mathrm{out}}^3}{G{S}_3\left(4+3m_{\mathrm{in}}/m_3\right)}}$$

(1.16)

This study focuses on developing a more accurate model for GW emission from HTSs by conducting a preliminary study on the detectability of such systems in our present and future detectors and providing a theoretical framework for the evolution of the system and modeling the GWs.

2. HTS Simulation

2.1. Constraints

2.1.1. Stability of Hierarchical System

For KL oscillations to occur, the stability of the three body system must be maintained. For this reason, a necessary criteria is mentioned in [19], which predicts a minimum separation between the inner and outer orbits for stability, beyond which the system might become chaotic. This condition is given by

$${\displaystyle \frac{a_{out}}{a_{in}}}>{\displaystyle \frac{2.8}{1-e_{out}}}{\left[{\displaystyle \frac{m_{out}}{m_{in}}}{\displaystyle \frac{1+e_{out}}{{(1-e_{out})}^{\frac{1}{2}}}}\right]}^{{\displaystyle \frac{2}{5}}}$$

(2.1.1)

This gives us a lower bound for the separation between the two orbits.

2.1.2. Relativistic Precession of Periastron

KL oscillations are suppressed or maximum eccentricity is restricted by the periastron precession of inner binary. To constrain this effect, we consider the GR precession timescale of periastron [6,20] of the inner binary and compare it with the KL timescale (see Equation (1.13)) to get the constraint whose scaled expression is of the form

$${\left({\displaystyle \frac{a_{\mathrm{out}}}{\mathrm{AU}}}\right)}^3<{10}^8{\left({\displaystyle \frac{1-e_{\mathrm{in}}^2}{1-e_{\mathrm{out}}^2}}\right)}^{{\displaystyle \frac{3}{2}}}{\left({\displaystyle \frac{m_{\mathrm{in}}}{M_{\odot }}}\right)}^{-2}{\left({\displaystyle \frac{a_{\mathrm{in}}}{\mathrm{AU}}}\right)}^4\left({\displaystyle \frac{m_3}{M_{\odot }}}\right)$$

(2.1.2)

2.1.3. Lense-Thirring Precession Effect

For the Lense - Thirring effect to appear, we need to have a constraint such that $t_{KL}\ge t_{L_{out}{S}_3}$ (see Equation (1.16)), for which the scaled expression [6] is given as

$${\left({\displaystyle \frac{a_{\mathrm{in}}}{\mathrm{AU}}}\right)}^{{\displaystyle \frac{3}{2}}}\le {10}^{-12}\left({\displaystyle \frac{m_3}{M_{\odot }}}\right){\left({\displaystyle \frac{m_{\mathrm{in}}}{M_{\odot }}}\right)}^{{\displaystyle \frac{1}{2}}}$$

(2.1.3)

assuming ${\chi }_3\sim 1$

2.2. Parameter Space and Final Models

We set up 17 mass combinations (see Table 1) and then for each combination, using the constraints, we find the region of stability (shown in Figure 1) and select the $a_{in}$ and $a_{out}$ of each model. Based on these initial models, we find all the initial parameters for these models (see Table 2). Note that we haven’t considered triple stellar mass HTS due to the possibility of their formation being much lower.

We shall use these parameters in the simulation to verify its working process and confirm the effects of KL dynamics occurring in these systems.

2.3. Modeling a Simple 3-Body System Using REBOUND

To simulate the dynamical evolution of HTSs, we used the open-source N-body integrator REBOUND [21]. REBOUND is widely adopted in astrophysical research for its flexibility, modular architecture, and high-precision integrators suitable for simulating gravitational interactions over a wide range of dynamical regimes. In particular, it supports various integration schemes, collision detection algorithms, and allows for custom initial conditions and external forces. All the integrators work on the basis of the Drift-Kick-Drift (DKD) scheme. The Hamiltonian of the system is given by $H=H_1+H_2$, where $H_1$ and $H_2$ represent the kinetic $\left(p\right)$ and potential $\left(q\right)$ parts, respectively. The DKD scheme integrates over a timestep, say dt, where half the timestep is used to integrate the kinetic part and the other, the potential part. For this study, we used the built-in integrator, IAS15 [22]. It is a 15th-order non-symplectic integrator (Newtonian) that can handle both conservative and non-conservative forces.

Our HTS are setup by initializing our masses alongside their respective parameter using the code. We then choose our time period over which to integrate the system and also the timestep to be used by the REBOUND integrator. In this case, we have kept the default timestep, $dt=0$. The code uses a loop to capture the positions and velocities of the masses in the system over the given time period over intervals of time determined by the required resolution for analysis. In all cases, resolution has been set to 500000. Hence, the number of time intervals is 500000. The extracted positions and velocities are used to calculate the eccentricity and inclination of the inner binary. This code makes up the Newt3 code.

We shall use REBOUNDx [23] to add additional physical effects, such as 1PN corrections, tidal forces, which can be incorporated using the REBOUNDx extension library. The 1PN correction to the Newtonian N-body Hamiltonian used by REBOUND can be easily augmented onto Newt3 by attaching the gr_full module to the same code. This module, during a simulation run, calculates the corrections using the instantaneous values of positions, velocities and modifies the accelerations of the masses at each timestep. This allows us to see the effect of GR precession on our system and how it affects our KL effect. This addition of the gr_full module to the Newt3 code creates the GR3 code. We shall discuss the results of our simulations in Section 4.

3. Detectability of HTS Through GWs

The noise curves of detectors play a crucial role in determining the detectability of various sources. These curves illustrate the sensitivity of the detector across different frequency ranges, providing valuable insights into its performance. To obtain these noise curves, we analyze noise measurements collected by sensors in the time domain. By performing a Fast Fourier Transform (FFT) on this data, we can calculate the Power Spectral Density (PSD) of the noise. Alternatively, an analytical equation for the noise specific to the detector can also be employed. For our purposes, the analytical noise $S_n\left(f\right)$ proves to be sufficient, allowing us to effectively assess the detector’s capabilities.

3.1. LISA

The Laser Interferometer Space Antenna (LISA) [24] is an ambitious gravitational wave observatory set to be launched by the European Space Agency (ESA), with participation from NASA, targeting an operational period in the 2030s. Unlike ground-based detectors like LIGO and Virgo, LISA will operate in space, allowing it to detect GWs at much lower frequencies. It comprises three spacecraft arranged in an equilateral triangle with 2.5 million kilometers between each pair [25]. LISA is designed to capture GWs produced by massive celestial events, such as mergers of supermassive black holes, inspirals of stellar-mass compact objects into supermassive black holes, and potentially signals from the early universe.

LISA’s low-frequency sensitivity opens a window to events beyond the reach of current ground-based interferometers. These events can reveal crucial insights into galaxy formation, cosmology, and fundamental physics. Orbiting the Sun in tandem with Earth, LISA’s constellation will effectively remove terrestrial noise sources, providing a pristine environment for GW detection. The data from LISA will complement observations from high-frequency detectors on Earth, expanding our understanding of the universe through multi-frequency GW astronomy. The sensitivity curve for LISA [25] is given by

$$S_n\left(f\right)={\displaystyle \frac{P_n\left(f\right)}{\mathcal{R}\left(f\right)}}$$

(3.1.1)

where $P_n\left(f\right)$ is the power spectral density of the noise, and $\mathcal{R}\left(f\right)$ is the response function of the detector.

Here, $P_n\left(f\right)$ is given by,

$$P_n\left(f\right)={\displaystyle \frac{P_{\mathrm{OMS}}}{L^2}}+2\left(1+{cos}^2\left({\displaystyle \frac{f}{f_{*}}}\right)\right){\displaystyle \frac{P_{\mathrm{acc}}}{{\left(2\pi f\right)}^4{L}^2}}$$

(3.1.2)

where $P_{\mathrm{OMS}}$ is the single optical metrological noise, and $P_{\mathrm{acc}}$ is the single test mass acceleration noise.

And, the response function $\mathcal{R}\left(f\right)$ is given by,

$$\mathcal{R}\left(f\right)={\displaystyle \frac{3}{10}}{\displaystyle \frac{1}{\left(1+0.6{\left(\frac{f}{f_{*}}\right)}^2\right)}}$$

(3.1.3)

Thus, the sensitivity curve $S_n\left(f\right)$ is

$$S_n\left(f\right)={\displaystyle \frac{10}{3L^2}}\left(P_{\mathrm{OMS}}+2\left(1+{cos}^2\left({\displaystyle \frac{f}{f_{*}}}\right)\right){\displaystyle \frac{P_{\mathrm{acc}}}{{\left(2\pi f\right)}^4}}\right)\left(1+{\displaystyle \frac{6}{10}}{\left({\displaystyle \frac{f}{f_{*}}}\right)}^2\right)$$

(3.1.4)

where $L=2.5\mathrm{Gm}$ and $f_{*}=19.09\mathrm{mHz}$.

Figure 2. Sensitivity curve of the LISA detector.

Figure 2. Sensitivity curve of the LISA detector.

3.2. LGWA

The Lunar Gravitational-wave Antenna (LGWA) is a pioneering gravitational-wave observatory proposed for deployment on the Moon. Unlike space-based interferometers such as LISA, LGWA utilizes the Moon itself as a giant solid-state antenna for detecting GWs (GWs). The mission concept involves an array of high-precision inertial sensors—Lunar Inertial Gravitational-wave Sensors (LIGS)—deployed in permanently shadowed regions (PSRs) at the lunar poles, where seismic noise is minimal and thermal conditions are extremely stable [5].

LGWA is designed to be sensitive in the decihertz frequency band, spanning from about 1 mHz to 1 Hz. This band bridges the gap between the millihertz sensitivity range of space-based detectors like LISA and the hertz–kilohertz range of ground-based detectors such as the Einstein Telescope and Cosmic Explorer. The lunar environment—with its cryogenic temperatures, lack of atmosphere, and seismic quietness—provides a uniquely favorable platform for low-frequency GW detection that cannot be achieved on Earth.

Through its novel detection mechanism, LGWA is expected to observe GWs from a wide range of astrophysical and cosmological sources. These include inspirals and mergers of intermediate-mass black hole binaries, binary white dwarf systems, and potentially the early inspiral phases of binary neutron stars and stellar-origin black holes, enabling early warnings for electromagnetic follow-up. LGWA also holds promise for detecting unmodeled transients such as tidal disruptions of white dwarfs, and constraining the cosmological stochastic GW background.

Unlike traditional interferometers, LGWA measures the Moon’s elastic response to passing GWs. Its sensitivity is characterized by the displacement noise of its inertial sensors and the Moon’s GW response function.

Figure 3. Sensitivity curves of the LGWA detector.

Figure 3. Sensitivity curves of the LGWA detector.

3.3. EccentricFD

EccentricFD is a frequency-domain waveform model implemented in the PyCBC library for simulating gravitational wave (GW) signals from eccentric binary systems. Unlike traditional quasi-circular approximants, EccentricFD includes the effects of non-zero eccentricity at leading-order post-Newtonian (PN) accuracy, enabling the modeling of binaries formed through dynamical capture, globular cluster interactions, or hierarchical triple-induced Kozai–Lidov oscillations. These systems may retain measurable eccentricities within the LIGO–Virgo–KAGRA frequency band, and waveform models such as EccentricFD are crucial for their detection and characterization [26,27,28]. The model is computationally efficient due to its frequency-domain implementation and is widely used for injection studies and parameter estimation in the analysis of non-circular compact binary coalescences.

The construction of the EccentricFD model is based on the stationary phase approximation, where the waveform is expressed as a sum over harmonics of the orbital frequency. The amplitude and phase of each harmonic are modified by the orbital eccentricity, which is specified at a reference frequency, typically 10 Hz. The model incorporates post-Newtonian corrections up to 3PN order in the binary’s orbital energy and gravitational wave flux. It takes as input the component masses, reference eccentricity, and luminosity distance, and returns the frequency-domain plus and cross polarizations [27,28].

This waveform model is used to get the frequency-domain waveform of the inner binaries from our various models and are then plotted onto the sensitivity curves of LISA and LGWA to identify whether the signal will be detectable in the detector. However, since we do not have access to the waveform response of LISA and LGWA, these plots just provide us the rough estimates of their detectability.

4. Results and Discussion

4.1. Simulations

We have run the initial parameters of our models (see Section 2.2, Table 2) through our Newtonian three-body code, Newt3, and our first-order Post Newtonian (1PN) GR three-body code, GR3. In most cases, we see that inclusion of GR at 1PN order increases the number of cycles observed over a certain period. In some cases, the system is stable but no KL oscillations occur and in others, the system is chaotic. We divide our simulation results based on their stability and the occurrence of KL oscillation in the HTS in order to understand the effects occurring in these systems much more clearly.

4.1.1. Stable HTS with KL Oscillation

From our simulations, we see that systems IM1, IM2, IM3 (see Figure 4) show negligible difference in the number of KL cycles occurring when simulated at both Newtonian and 1PN level through the Newt3 and GR3 codes, respectively.

But as soon as the tertiary mass range goes into the SMBH range, we see large changes in the simulation results given at Newtonian and 1PN order, as seen in systems SM1, SM2, SM3, SM6, SM7 (see Figure 5). In SM6 and SM7, we also see that in the GR3 simulation result, there appears to be a lot more fluctuation in the eccentricity of inner binary, as it approaches the maximum eccentricity allowed by the KL mechanism.

In the case of system SM4, we see a drastic difference in the number of KL cycles (see Figure 6). This could be due to the initial mutual inclination, $I={45}^{\circ }$, being quite close to the lower limit of the effective limit range of mutual inclination for KL oscillations to occur.

4.1.2. Stable System But No KL Oscillation

We observe a curious case in system IM4 (see Figure 7). When simulated (over 20 years), we see that in both simulations, the inner binary does not exhibit any KL oscillations but the system is stable. The results from GR3 show that the inner binary is precessing at a wider angle about the tertiary compared to the Newt3 results, clearly seen in the simulation run for 200 years

4.1.3. Instability in Triple System

Systems IM4, IM5 and SM9 exhibit unstable, chaotic behaviour despite us considering parameters that are supposed to describe stable HTSs. In system I5, we see that all three bodies are unbounded in under 5 years whereas in systems SM9, the inner binary is ejected from the system instead of the creation of an unbounded system (see Figure 8).

4.2. Eccentric System Detectability

From our GW frequency domain plots (see Figure 9), we see that all our models come under the LISA band when considering their eccentricity to be very high. We also see that as eccentricity reduces, the signal comes towards the LIGO band. And since we know that binaries at high eccentricities circularize faster, we would be able to see their complete evolution over the LISA, LGWA and LIGO bands.

From Figure 10, we see that as eccentricity reduces, so does the number of harmonics observed in the band as a consequence of equation (1.1) and as a result, decreasing power loss from the system with decreasing eccentricity.

We can also consider future ground based detectors such as the Einstein Telescope (ET) [29] into the picture to create a continuous detector array for analyzing the complete evolution of the inner binary (see Figure 11). These detectors will allow us to capture the residual harmonics from the binary as its orbit circularizes.

5. Conclusion

This study investigates the dynamical evolution of hierarchical triple systems (HTSs) and evaluates their gravitational wave (GW) detectability through numerical simulations. By selecting models from astrophysically relevant parameter spaces and implementing them using REBOUND and REBOUNDx, we confirm the occurrence of KL oscillations and examine the competing influence of general relativistic (GR) precession. Our results highlight the importance of including post-Newtonian corrections, as GR precession can significantly modulate both the timescales and amplitude of KL-driven eccentricity growth. Interestingly, despite choosing initial parameters expected to yield stable HTSs displaying KL oscillations, our simulations reveal both stable and unstable systems.

Within the stable systems, two distinct behaviors emerge: systems exhibiting KL oscillations and systems that do not. For systems showing KL oscillations, the impact of including 1PN dynamics increases with the tertiary mass, particularly in the supermassive black hole (SMBH) regime, where KL oscillations become more frequent. This effect arises because GR precession can partially restart the KL mechanism even under imposed constraints, though its influence is weak [11].

We also identify notable outliers in our simulations. A stable system (IM4) does not display KL oscillations, while two unstable systems (IM5 and SM9) show unexpected chaotic behavior. In these cases, one inner binary mass is an order of magnitude larger than the other, but the tertiary mass differs: in the stable system, it is two orders of magnitude larger, whereas in the unstable systems, it is only one order larger. This suggests the existence of an additional stability constraint dependent on the mass ratios within the HTS. The stable system may lie within this constraint, explaining its stability and absence of KL oscillations, whereas the chaotic systems exceed it. Future work should aim to quantify this constraint to deepen our understanding of HTS stability.

From frequency-domain gravitational waveforms obtained using EccentricFD, we find that highly eccentric binaries driven by KL oscillations can emit GWs detectable by both current and future detectors. Specifically, LISA is well-suited to capture signals during the high-eccentricity inspiral phase, while aLIGO is more sensitive to the circularized final mergers. Proposed detectors such as LGWA, with a frequency band bridging LISA and LIGO, alongside the Einstein Telescope (ET), could identify residual eccentricity in these systems. This multi-band detectability underscores the potential for HTS-originating binaries to serve as complementary probes across a wide range of GW frequencies.

We acknowledge that our simulations of hierarchical triple systems (HTSs) are idealized and do not capture all physical effects present in realistic astrophysical environments. Additionally, the gravitational waveform curves presented here do not incorporate the detector response functions, and thus do not reflect how these signals would appear in actual observational data. Nevertheless, this study provides a foundational framework for constructing more sophisticated gravitational wave (GW) templates that include complex three-body dynamics and relativistic corrections. Ongoing development of our HTS code is aimed at incorporating higher-order PN terms, up to at least 2.5PN, to better understand the interplay between KL oscillations, GW emission, and other secular effects. These improvements will ultimately support the generation of more accurate, time-evolving GW waveforms for inner binaries embedded in hierarchical triple systems.

Data sharing is not applicable to this article as no datasets were stored or separately generated during the study. All data were produced and utilized in real-time within the simulation code for analysis and were not saved externally.

The underlying code for this study is not publicly available but may be made available to qualified researchers on reasonable request from the corresponding author.

All authors declare no financial or non-financial competing interests.