Gravitational Waves from Hierarchical Triple Systems with Kozai-Lidov Oscillation

We study gravitational waves from a hierarchical three-body system up to first-order postNewtonian approximation. Under certain conditions, the existence of a nearby third body can cause periodic exchange between eccentricity of an inner binary and relative inclination, known as Kozai-Lidov oscillations. We analyze features of the waveform from the inner binary system undergoing such oscillations. We find that variation caused due to the tertiary companion can be observed in the gravitational waveforms and energy spectra, which should be compared with those from isolated binaries and coplanar three-body system. The detections from future space interferometers will make possible the investigation of the gravitational wave spectrum in mHz range and may fetch signals by sources addressed.


INTRODUCTION
The LIGO-VIRGO collaboration observed eleven gravitational wave signals from compact binary mergers during the first and second runs (O1 and O2) [1]. The third observation run (O3) by LIGO and VIRGO began in April 2019 and new detection alerts have started following within months of operation [2]. With current sensitivity design, ground detectors are focusing on the high-frequency range (10Hz to 1000Hz). Most expected sources lying in this range are compact binaries in the merger and ringdown phase. For ground-based detectors, currently, we have five ground-based detectors with characteristic strain of order 10 −22 ; advanced LIGO (aLIGO) detectors in Hanford and Livingston[3], advanced VIRGO (aVIRGO) in Italy [4], KAGRA in Japan [5], and GEO600 in Germany [6].
Future space observatories like Laser Interferometer Space Antenna (LISA) are expected to explore low frequency range from 10 −4 Hz to 10 −1 Hz with characteristic strain of order 10 −21 [7] whereas DECi-hertz Interferometer Gravitational wave Observatory (DECIGO) is aiming to fill the gap between LIGO and LISA with frequency band around 10 −2 Hz to 10 Hz with characteristic strain down to 10 −24 [8]. Several big projects have also been proposed in space; Big Bang Observer (BBO) [9], advanced LISA (aLISA) [10], and TianQin and TAIJI in China [11,12]. With the launch of those space detectors, we can aim for low-frequency sources. It will be possible to map a source all the way from inspiral to a merger phase by incoming gravitational waves. The observations so far, and the future GW detections may shed light on binary formation channels, enable precision tests of general relativity (GR) in the strong-field regime, and initiate new avenues in astrophysics.
In this new stage of the beginning of gravitational wave astronomy, we have to list up all possible gravitational wave sources and analyze what we expect from those in observations. So far, we have seen intensive analysis for binary systems composed of objects with various masses [13][14][15], and references therein]. However, we may expect more exotic sources in nature. One of the likely sources is a three body system, which we will study in this paper. The environment near supermassive black holes (SMBHs) in galactic nuclei comprises of a large number of stars and compact objects. Triple systems could also emerge in these surroundings [16-20, and references therein] with the likelihood of the presence of hierarchical structure. The potential astrophysical sources of gravitational waves involving three-body system have seen increasing interest over time [21][22][23][24][25]. All the detections so far (direct or indirect) came from binary inspirals and merger phase. It is good to contemplate how the waveform would look like if the binary system was accompanied by a tertiary companion in a hierarchical arrangement. The hierarchical triple body approximation has useful applications and can be applied to diverse settings, varying from planetary and stellar mass scales to SMBHs. In this approximation, the energy of each orbit is separately conserved and therefore the two semi-major axes are constants. The Kozai-Lidov (KL) mechanism [26,27], which is one of the most important phenomenon in a hierarchical triple system, is particularly exciting. KL mechanism occurs when the inner orbit inclines enough with respect to the outer orbit. The main feature of KL mechanism is the secular change (timescale much larger than the orbital periods and referred as KL timescale) of the eccentricity of the inner orbit and the relative inclination between inner and outer orbit. Both values oscillate inversely with respect to the other, that is, when the eccentricity decreases, the inclination increases, and vice versa. The eccentricity of the orbit can reach extreme values leading to various astrophysical phenomena. For instance, if the eccentricity becomes large enough because of KL oscillations, the merger rate of BHs can be enhanced by large emission of gravitational waves [28][29][30][31], the tidal disruptions rate in a stellar cusp containing SMBH binary can be several orders larger in magnitude than expected in case of single SMBH [32][33][34][35][36]. Further, the formation of hot Jupiters or ultra-short-period planets is studied through the secular gravitational interactions between two planets in eccentric and inclined orbits [37][38][39][40].
Recently, there has been extensive work focusing on gravitational wave detections from such systems. The probability of detection and frequency range from different sources has been discussed by many authors. The secular evolution of stellar-mass BH binaries near SMBHs is studied, focusing on the tidal perturbation by the SMBH [41,42]. The analysis shows that merging binaries will enter the LIGO observational window while on orbits that are still very eccentric (e ∼ 0.5). The efficient GW analysis for such systems would, therefore, require the use of eccentric templates. The constraints are put on the importance of KL induced mergers for producing gravitational wave sources detectable by aLIGO, for triples with an inner pair of stellar-mass BHs. Another interesting model with triple SMBHs is studied [43]. Numerical simulations of the dynamics of triples within galaxy cores exhibit phases of very high eccentricity (as high as e ∼ 0.99) showing the likelihood of detection of GWs by LISA.
Some authors discuss variations in GWs in the presence of the third body due to relativistic beaming, Doppler, and gravitational redshift [44][45][46][47]. The acceleration imparted by the hierarchical companions leads to potentially observable changes in the waveforms that would reflect their gravitational interactions with the surrounding matter. The GW signal would be distinguished by a net phase change or even a time-dependent Doppler shift arising from the orbital motion. This could provide direct information about the black hole binary environments and otherwise invisible ambient mass.
Recent work analyzes triple systems composed of the SMBH near the center of M87 and a pair of BHs with different masses, from stellar to intermediate-mass BHs (IMBHs) [48]. They identify frequency peaks and examine the detectability of GWs with LISA, µAres, and DE-CIGO by integrating over the observation time within the lifetime of different GW detectors.
Indirect observation of GW from a triple system is also studied by analyzing the cumulative shift of periastron time of a binary pulsar undergoing KL oscillations [49].
In this paper, we focus on the features of the gravita-tional waveform generated by the inner binary in a hierarchical triple system and draw comparisons with the isolated binary systems. We also discuss the detectable parameter range and observation of such sources by future detectors. The paper is organized as follows: We review the characteristics of the hierarchical structure and KL-mechanism in §II A. In section §II B , we discuss the constraints on KL oscillations due to relativistic effects.
In section §III, we recall the fundamentals of Quadrupole formalism and discuss GWs from highly eccentric orbits.
In section §IV , we classify models and explore parameter space for KL oscillations. We show results for GWs and discuss the observation possibility in §V. The paper is concluded in §VI. In Appendix A, we provide the waveforms from isolated binary systems and hierarchical triple coplanar systems for comparison.

II. DYNAMICS OF TRIPLE SYSTEM
In this section, we sequentially explain a hierarchical triple system, Kozai-Lidov mechanism, and their stability and constraints.

A. Hierarchical Triple System and Kozai-Lidov Mechanism
We can have several stable configurations and shapes when it comes to a three-body dynamical system. Broadly addressing, it can be concentric, coplanar orbits or inclined orbits in a hierarchical arrangement. In this paper, we treat the so-called hierarchical triple system, whose schematic picture is given in Fig. 1.
FIG. 1: The hierarchical triple system is constructed from inner and outer binaries. The inner binary consists of objects whose masses are m1 and m2, and the outer one is the pair of the inner binary and the third body with mass m3.
As depicted in the figure, we can decouple a system into the two-body inner binary orbit and the outer orbit with tertiary companion in motion around center of mass of inner binary, if the distance between the first and second bodies is much shorter than the distance to the third body. This hierarchy is supported by assuming that the gravitational effect of the third body is much smaller than the gravitational interaction between the first and second bodies.
In a two-body problem in Newtonian dynamics, a bound system is given by an elliptic orbit, which is described by six orbital elements; the semi-major axis a, the eccentricity e, the inclination i, the argument of periastron ω, the longitude of ascending node Ω, and the mean anomaly M . Although these elements are constant (except mean anomaly) in an isolated two-body system, in the hierarchical three-body system, the perturbations from the tertiary companion affect the binary motion and modify the trajectory from that of the isolated one. Such a trajectory is not closed in general, but we can introduce an osculating orbit at each time, whose trajectory is approximated by the elliptical orbit with the above six orbital elements determined by the instantaneous position and velocity [50]. As for the outer orbit, we pursue the center of mass of the inner binary rotating around the tertiary companion (see Fig. 1). It can also be described as another osculating orbit. Hence, we introduce two osculating orbits, which are called as inner and outer orbits: The masses of the inner binary are m 1 and m 2 , while the tertiary companion has the mass m 3 . We use subscriptions 'in' and 'out' to show the elements of inner and outer orbits, respectively.
In the hierarchical three-body system with large relative inclination, a secular change of orbital elements may occur. Under some conditions, there appears an oscillation between the relative inclination and the eccentricity of the orbit, which is called Kozai-Lidov (KL) oscillation [26,27]. This mechanism was independently discovered by Kozai and Lidov in 1962. The relative inclination I is defined as the argument between the inner and outer orbital planes. It is given by two orbital elements as cos I = cos i in cos i out + sin i in sin i out cos (Ω in − Ω out ) .
If we assume m 1 m 2 , m 3 , we find that the oscillation occurs in secular time-scale under conservation of energy and angular momentum. In Newtonian quadrupole approximation method [51], it results in the secular exchange of e in and I with the conserved value of Θ, which is defined by This approximation also gives the criterion of KLoscillation for circular orbit as which is equivalent to 39.2315 • I 140.7685 • . This KL type oscillation will occur even if m 1 is not much smaller than m 2 , m 3 [37]. When we include a relativistic effect, we also find the similar oscillation [30,52]. The effect of general relativity (GR) changes these KL-criterion as [53].
where the GR correction term H GR is derived from the double-averaged post-Newtonian Hamiltonian of twobody relative motion [54] as where m in = m 1 +m 2 is the total mass of the inner binary.
The KL-oscillation time-scale is also modified when the GR effect is taken into account [52].
Note that it has been pointed out [55,56] that it is imperative to take into account the 'cross terms' between the Newtonian perturbation and the Post-Newtonian precession. In our numerical evolution, we integrate the equations of motion (up to 1PN) directly. The effects of the cross terms are automatically included.

B. Stability and Constraints
In this section, we discuss stability criteria for a hierarchical triple system. We also give constraints for KL oscillation to occur, which are obtained from the conditions such that KL-oscillations are not suppressed by relativistic effects, i.e., the relativistic precession of the periastron of the inner binary and the Lense-Thirring (LT) precession effect due to the rapid rotation of heavy third body.

Typical time scales
Before discussing stability and constraints, we first present typical time scales of a hierarchical triple system. The orbital periods of the inner and outer binaries are given by where m out = m 1 + m 2 + m 3 is the total mass of a triple system. The KL oscillation time scale is given by [57] t KL ∼ P in m in m 3 a out a in 3 1 − e out 2 3 2 .

Stability of Hierarchical Arrangement
For KL oscillations to occur, the first aim is to maintain the stability of a three-body system. A necessary (not sufficient) criteria is given by [28], which predicts a minimum separation between inner and outer orbits for stability, beyond which a system may cause a chaotic instability. Hence, we have one stability condition, For simplicity, throughout this paper, we set outer orbit not to be so eccentric such that e out 1. Thus the scaled expression is written as, This expression gives a lower bound to the separation between two orbits.

Constraints from the relativistic effects
3-1. Relativistic Precession of Periastron KL oscillations are suppressed or maximum eccentricity is restricted by the relativistic precession of the periastron of the inner binary. An analytical expression can be derived using double averaged approximation up to quadrupole order to estimate parameter space such that KL oscillations are not destroyed by precession [28]. The GR precession timescale of the periastron of the inner binary can be estimated as, where r g,in ≡ Gm in /c 2 is the gravitational radius of the inner binary. Comparing this with the KL time scale t KL , we find the constraint For clarity, the scaled expression takes the form: assuming e out 1.

3-2 Lense-Thirring Precession Effect
Recent studies have shown changes in KL oscillations caused due to GR effects involving a SMBH [58][59][60]. If the third body is a rapidly rotating SMBH, then LT effect might become important changing the evolution of eccentricity excitation from the usual KL oscillation. It appears in the 1.5 PN order. For a rotating black hole, which is the third body in our setting, the spin angular momentum is given as S 3 = χ 3 Gm 2 3 /c where χ 3 is the spin parameter, and outer orbital angular momentum is given as where the reduced mass of the outer orbit is defined by The time scale of the orbit-averaged precession of L out around S 3 is given as The KL oscillations can be affected when t LoutS3 becomes comparable to the KL oscillation time scale t KL . Hence, in order to observe the usual KL oscillation, we have the constraint such that t LoutS3 > t KL , for which the scaled expression is given as: assuming χ 3 ∼ 1. If the ratio is comparable or greater than 1, then LT precession dominates resulting in chaotic evolution of eccentricity as shown in [58] and hence more dissipation of gravitational waves, which may be interesting for further study. However, in this paper our discussion is limited to 1PN effect and we focus more on the conventional KL oscillations in the inspiral phase.

III. GRAVITATIONAL WAVES
Under this section, we set the framework by recalling the famous quadrupole formula and definition of 'near coordinate zone' (NCZ). We further discuss GWs from highly eccentric orbits, which are expected in KL oscillation.

A. Quadrupole Formula
In GR, the leading order contribution of the observed gravitational waves is given by the quadrupole formula. It allows the evaluation of gravitational wave energy and waveforms emitted by a dynamical system. In a wave zone, the gravitational waves denoted by a small deviation of the metric h ij are described by where TT refers to transverse-traceless gauge and a dot denotes the time derivative. D denotes the distance of the source from the observer. The reduced quadrupole moment Q ij is defined by where i, j, · · · run from 1 to 3 (x 1 = x, x 2 = y, x 3 = z), and the mass quadrupole moment M ij for N masses is expressed as m A denotes the A-th mass at the location x i A with origin at the center of mass of a system.
The power of GW emission is given by where the bracket denotes the Brill-Hartle averaging, which is taken over several wave length [61]. The GW amplitude from a circular binary system consisting of masses m 1 and m 2 with orbital radius a is approximately given by which is described by the scaled expression as Another important quantity is energy spectra of the GWs. For an isolated circular orbit, the energy radiated in gravitational waves is concentrated in the second harmonic such that where f in is the orbital frequency of the inner binary given by We find it in the energy spectrum, whose example is given in Appendix A 1 and Fig. 20. Large eccentricities induced in the inner binary during KL oscillation and variation of eccentricity with time produce more harmonics in gravitational waves, which make the energy spectra rich as compared to the spectra from an isolated circular binary as well as an enhancement of the wave amplitude as shown in the next subsection.
In the derivation of the quadrupole formula, a system must be contained in its NCZ , i.e., in order to use the quadrupole formula, a system should be inside a region (centered on the origin of the coordinates) of radius comparable to (or smaller than) the GW wavelength λ (see some good example of what not to do given in [62]). Therefore, we must keep a check that a hierarchical triple system is in its NCZ when the formalism is applied to the center of mass of the triple system for the model in consideration. There are two NCZs for a hierarchical triple system, whose sizes correspond to two wave lengths: one is that of an inner binary λ in = a in c/v in and the other is that of an outer binary with λ out = a out c/v out , where v in and v out are typical velocities of the inner binary and the outer binary, respectively. Since we focus on GWs from the inner binary, the NCZ condition is given by a out < λ in . An inaccurate application may lead to apparent and unphysical behavior as shown in Fig. 18 of [63].

B. Gravitational waves from highly eccentric orbit
In this subsection, we discuss the GWs from a highly eccentric binary system. When the KL oscillation occurs, the inner binary system is supposed to be in an inspiral phase. Hence, it is usually expected to be difficult to observe the GWs unless the inner binary is just before a coalescence. However, the eccentricity in the KL oscillation can become very large, which enhances the amplitudes of GWs and those frequencies. As a result, we may have a chance to observe the GWs even in an inspiral phase.
Gravitational Waveform from a binary system depends strongly on the eccentricity of the orbit. In an eccentric orbit, the emission of gravitational waves is higher and distributed over many harmonics. The GWs from an eccentric binary was studied by Peter and Mathews [61]. The total radiation power P averaged over one period is given by 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 when the eccentricity becomes large. The total radiation power P is the sum of the power radiated in the n-th harmonics P n as The radiation power P n in the n-th harmonic is given by where J n (x) is the Bessel function of the first kind. The spectrum P n is peaked at the GW frequency f peak , where the magnification factor n peak is well approximated by for the range of 10 −6 < 1 − e 2 < 1 [30]. For e 1, n peak is just 2, while considering high eccentricity, say e ∼ 0.98 , we find n peak ∼ 574. As the eccentricity increases, change in f peak becomes significant and may lie in the detectable frequency band 10 −4 Hz ∼ 10 kHz. It is worthwhile evaluating the peak frequency as

Hz
Note that although the n-th harmonics power P n has a peak at f = f peak , the spectrum is rather broad. As a result, the observable range of the frequency can be wide as shown in [47].
In observation, another important quantity is the wave amplitude. The GW amplitude h n of the n-th harmonic is related to P n by where the bracket denotes the time average over a few orbital period, and f n is the frequency of the n-th harmonic wave.
When we evaluate at the peak frequency f n = f peak , we find Suppose that m 1 = m 2 = 10M , a in = 0.01 AU, and D = 10kpc. When the eccentricity gets large enough via KL oscillations, we obtain for e in = 0.98. It may be observable.
To answer the question whether it is really observable or not, we have to perform a numerical simulation and evaluate the dimensionless characteristic strain, which is defined as [64], To calculate this quantity, we first find gravitational waveforms h + (t) and h × (t) and transform them into the Fourier components h + (f ) and h × (f ). In our analysis, we apply Fast Fourier Transform (FFT) to perform the Fourier transformation, and evaluate the characteristic strain.
The total GW energy E = dt P consists of many harmonics, which can be described as The energy spectrum dE df is obtained as

IV. CLASSIFICATION AND CONSTRAINTS ON PARAMETERS
Before we perform our analysis, we classify the possible models with the KL oscillations by the masses of a triple system. We focus on stable hierarchical triple systems and list the models in three different mass ranges. This general discussion helps in visualizing several realistic situations varying on the scale of mass.
We have three possible cases: Case (A): m in ∼ m 3 , in which three bodies have comparable masses, or one of the inner binary system can have a smaller mass.
Case (B): m in m 3 , which mimics the state of a less massive inner binary around a supermassive black hole (SMBH) in the galactic nuclei. The inner binary can be a stellar mass binary system such as a binary black hole (BBH), a binary neutron star (BNS) or a neutron star-black hole (NSBH) binary. It can also be a BBH consisting of an intermediate BH (IMBH).
Case (C): m in m 3 , which is the reverse of Case (B), hence, instead of SMBH, we may have a planetary object or an asteroid as a tertiary companion. For instance, the inner binary consists of intermediate BHs, while the third companion is stellar mass size.
In Case (C), although the KL oscillation is possible for appropriate parameters, the oscillation occurs between the outer eccentricity and the relative inclination. The eccentricity of the inner binary may not be affected so much. As a result, the gravitational waves from the inner binary system, which are much more important than those from the outer orbit, has almost no effect due to KL oscillation. We then do not analyze Case (C) in this paper. There is another reason we do not discuss Case (C), i.e. KL timescale is much longer and the cycle may not be captured in one observation run.
However, our previous work [49] discusses the indirect detection of KL oscillations. Therefore, if a third body is a pulsar with the inner binary of intermediate-mass black holes then the presence of KL oscillations can be captured in radio signals from the pulsar. Now, we show the parameter range in which the KL oscillation appears and is stable.
As discussed in §II B, we have stability conditions and constraints in order to find the KL oscillations. In this case, the relation between the typical time scales are given as We present some examples in Fig. 2, in which we show the allowed region in the parameter space (a in , a out ) for fixing the the masses of a triple system as m 1 = m 2 = m 3 = 10M . In §II B, we show the condition with which the orbit does not become chaotic. This gives the maximum value of a in , which is given by green line [Eq. (2.2)]. The other two constraints come from the relativistic effects. One is the precession of the periastron, whose condition is shown by the blue line [Eq. (2.3)]. Beyond this line, such relativistic effect will suppress the KL oscillation. The other is the LT precession effect, whose critical line is out of the present parameter range. Taking into account the above conditions, in Fig. 2, we present the possible range in the parameter space of (a in , a out ) for the stable KL oscillation, which is shown by the lightblue shaded region. For the present parameter choice,  In this case, the relation between the typical time scales are given as In Fig. 3, we present the possible range in the parameter space of (a in , a out ) for fixing the masses of a triple system as m 1 = m 2 = 10 3 M , and m 3 = 10 6 M . The constraints are determined only by the relativistic precession effect of the periastron and chaotic boundary just as Since the range of mass parameter in Case (B) is very wide, in Fig. 4, we also present the possible range in the parameter space of (a in , m in ) for the stable KL oscillation, which is shown by the light-blue shaded region. We fix the mass of tertiary component and the semi-major axis of the outer orbit as m 3 = 10 6 M and a out = 10AU. For the present parameter choice, the constraints are determined by the relativistic precession effect of the periastron and chaotic boundary condition just as Case (A).
However, note that the LT precession may become significant in Case (B) when the third body is much heavier. See Fig. 5, in which we set m 3 = 10 9 M . In this case, a hierarchical triple system will evolve in the left direction horizontally as a in decreases because of the emission of GWs. If m in 10 3 M , the evolution curve will hit on the LT critical line before the relativistic periastron shift becomes important. Hence the system will evolve into a chaotic KL oscillation phase, which may be very interesting to study, but it is out of our present analysis because the LT precession appears in 1.5 PN order.   We then choose the semi-major axis a in and a out in the stable regions for the KL oscillation, which are given in Figs. 2-4. In addition, to confer the results of our analysis, we pick Case (A) to contrast the key features in waveforms from an isolated binary system (Models ABC and ABE) and in waveforms from a binary in a coplanar hierarchical triple system (Models ACC and ACE). The parameters are summarized in Table I. The observer's location (ι) is defined as the angle between the detector and the "axis" of the reference plane (i.e., the "normal" to the initial outer orbital plane). As we are focused on the inspiral phase of the inner binary, there are special features that can be seen in the waveform.

B. Method of Orbit Evolution
The first order post-Newtonian equations of motion also known as the Einstein -Infeld -Hoffmann equations of motion [65], describes the dynamics of a system of point-like masses due to their mutual gravitational interactions, including GR effects. We employ Eq. (5.1) in order to solve the three-body system. This equation could also be derived from the Lagrangian given by Lorentz and Droste [66].
where m k , v k , x k (k = 1, 2 and 3) are the mass, velocity and position of the k-th component of the system, G is the gravitational constant, and c is the speed of light. Eq. (5.1) has been numerically integrated by using 6-th order implicit Runge-Kutta method, whose coefficients are obtained from Butcher [67].
On integrating, we obtain the numerical data of positions and velocities of the triple system. The initial outer orbital plane is considered to be the reference frame. In order to set up initial conditions, we convert initial orbital elements of inner and outer orbits into the variables x k and v k in Cartesian coordinates, with its origin in the center of mass of the whole system [50]. We integrate the EIH equations (5.1) numerically and evaluate the osculating orbital elements at each step from the numerical data of positions and velocities of the triple system (see e.g. Murray and Dermott [50]). Since the inner orbit is not exactly an ellipse, the obtained osculating elements are oscillating with small amplitudes in the cycle of inner orbit. Hence, we take an average of the osculating elements for each cycle of inner orbit and then obtain the averaged semi-major axesā in ,ā out and eccentricitiesē in , e out , which may give the effective values of the orbital elements. Those elements will evolve secularly in time because of the effect of the tertiary body.

C. Numerical Results
We first summarize the numerical results for Model AKL1 in next subsections 1, 2 and 3, and then mention about Model BKL1 in subsection 4.

Evolution of the eccentricity
We first show the result of the evolution of eccentricity and relative inclination of the inner binary for Model AKL1, in which the masses of a triplet are m 1 = m 2 = m 3 = 10M and initial semi-major axis are a in = 0.01AU and a out = 0.1AU. The typical time scales are given by P in ∼ 0.08 days, P out ∼ 2.11 days and t KL ∼ 164 days.
The initial eccentricities are e in = e out = 0, but the relative inclination is chosen to be I = 90 • . We show the time evolution of the eccentricity e in and the relative inclination I for one KL oscillation cycle in Fig 6. We have run our code for several values of inclination varying from 40 • to 90 • , for which we obtain stable KL oscillations.
For I = 90 • , we find the maximum value of the eccentricity as e in,max ≈ 0.987, while e in,max ≈ 0.02 for I = 40 • . For lower values for relative inclination, we find lower maximum eccentricities as expected from Eq. (2.1).
FIG. 6: The time evolution of the osculating orbital elements (the averaged eccentricity of the inner binaryēin,and averaged relative inclinationĪ) in Model AKL1 (m1 = m2 = m3 = 10M , ain = 0.01AU, and aout = 0.1AU) for one KL oscillation cycle. The red and blue lines show the evolution ofēin andĪ, respectively. When the inclination drops to Imin ≈ 55 • , the eccentricity reaches the maximum value emax ≈ 0.987. The small oscillations are caused by the outer companions motion.

GW Waveform
Next we show the waveform of the GWs for Model AKL1 (m 1 = m 2 = m 3 = 10M , a in = 0.01AU, and a out = 0.1AU). In Fig. 7, we present the waveform of two polarization modes for one KL oscillation cycle. For the small eccentricity until day 98, the amplitudes are low as h ∼ 10 −21 , but when the inner binary is in high eccentricity regime around day 100, we find a significant rise in amplitude as h ∼ 10 −19 . To see more detail, we enlarge the part of waveforms. In Fig. 8, we show the waveform of + -polarization mode for two days. The top figure depicts the data for the first two days, when the eccentricity is small, while the bottom one shows the data from day 106 through day 108, when the eccentricity becomes very large. We also show the similar enlarged waveforms for ×-mode in Fig. 9 FIG. 8: The zoomed in parts of waveform of + polarization in low and high eccentricity regime. The top figure shows the wave form in low eccentricity period (the first two days), while the bottom one depicts the wave form in high eccentricity period (two days from day 106 through day 108). These features can be easily understood when we analyze more simple systems such as an isolated binary or a coplanar motion of a hierarchical triple system, for which the results are summarized in Appendix for Models ACC-ABE (see Table I).
For an isolated circular binary, the waveform is sinusoidal, but if the eccentricity is large as e = 0.9, the sharp peak in the waveform appears near the periastron point. The amplitude also increases compared with the circular case. For a hierarchical triple system, if the orbit is coplanar, we find similar features with those in an isolated binary. When the eccentricity e in is small, we find the sinusoidal wave form, while if e in is large, a sharp peak appears near the periastron point. In the case of a triple system, however, there appears small modulation in the amplitudes because of the effect from a tertiary component. In Model AKL1, when the eccentricity is small, the wave form is almost the same as that in a coplanar circular orbit of a triple system. When the eccentricity becomes large via KL oscillation, the local feature of the wave form is similar to that of an eccentric binary or in a coplanar eccentric orbit of a triple system. As shown in Figs. 10 and 11, the amplitude is modulated randomly because of the oscillation of the eccentricity, which is different from the coplanar eccentric case.
For the observer at a different position, we also find the same features. The difference appears in the global shape of the wave form. We show the examples for the observers at ι = 45 • and ι = 90 • measured from the initial outer orbital plane in Figs. 10 and 11.

GW Spectra
Large eccentricities during KL oscillation and variation of eccentricity with time introduce more harmonics in GWs which make the energy spectra rich as compared to the spectra from an isolated binary (shown in Appendix A 1). In Fig. 12, we show energy spectra for two stages in one KL cycle. During low eccentricity regime, the spectra is similar to that of circular coplanar hierarchical triple (Model ACC).
FIG. 12: The GW energy spectra from the inner binary with KL oscillation for low and high eccentricity regime. The top figure shows the spectrum at day 1-5 (ein ∼ 0), and the bottom one gives that at day 106-108 (ein = 0.983 ∼ 0.987). We find a broad band spectrum for a high eccentric orbit.
In the top figure in Fig 12, which is obtained by the Fourier transformation for the time interval from day 1 to day 5, we find two sharp peaks which correspond to the two periods of the inner and outer orbits. On the other hand, during high eccentricity regime, as we see from the bottom one of Fig 12, which is obtained for the time interval from day 106 to day 108, the spectra is quite similar to the eccentric coplanar model ACE (see Fig.  23). Because of high eccentricity, many higher harmonics appear. As a result, we find a broad band spectrum for a high eccentricity regime.

Model BKL
We also perform numerical calculation for Model BKL1, i.e., m 1 = m 2 = 10 3 M , a in = 0.15AU, and a out = 10AU. We find that the properties of GWs are very much similar to Model AKL1. There appears KL oscillation, in which the maximum eccentricity becomes 0.99. As for the waveform, it is a sinusoidal shape with modulation caused by a tertiary object when the eccentricity is small, while it is deformed during the large eccentricity phase and shows a sharp peak by passing near a periastron point. The GW spectra also change in one KL cycle from ones with two sharp frequencies, corresponding to two orbital frequencies, to a broad band spectra. The difference appears only in the time scales given in Table I. The observability is also different from Model AKL as we will show in the next subsection.

D. Observability: Frequency, Strain, and SNR
To examine the possibility of observation of the present triple models by future detectors, we have to check whether the frequency and strain are in the observable range and evaluate the signal to noise ratio (SNR).
We first plot f peak for Model AKL1. We show in Fig.  13 that due to KL oscillations the frequency sweeps up to the detectable range of proposed future space interferometers. However, as the inner binary is in inspiral phase, the GW strain may lie below the sensitivity curve of detectors. We then plot strain versus frequency curve for Model AKL1 in Fig. 14.
The black, brown, blue, green curves correspond to strain sensitivity curves of LISA, aLISA, DECIGO and BBO respectively [10, 64,68]. We plot GW characteristic strain curve for Model AKL1 at different days. The orange curve is plotted for day 102-103 when the eccentricity is e in ∼ 0.944. The red curve is plotted for day 103-105 with e in ∼ 0.972. The violet curve is plotted for day 106-108 with e in ∼ 0.987, which is the maximum value in one KL cycle. After day 108, the eccentricity decreases and then the strain curves decrease to the orange one through the red one. These curves show the evolution of the strain in one KL cycle. At high eccentricity regime of KL cycle, the GWs may be observable for at least two days every 164 days (one KL cycle) by DECIGO or BBO.
In Fig. 15, we also show the similar plot for Model AKL2 with m 1 = m 2 = m 3 = 30M and a in = 0.01AU. When the observer distance is again kept to be at 10kpc, the signal is louder due to increase in mass.
One may wonder that the GW signal may be louder when a in becomes smaller via an evolution of the inner binary by the emission of GWs. For example, one may choose a in = 0.006AU, which is still in the stable range for KL oscillations (see Fig. 2). However it does not give a louder signal. It is because the maximum eccentricity e in,max = 0.87 for a in = 0.006AU is not high enough to give an observable strain compared with Model AKL1 or AKL2. We show the change of the maximum eccentricity e in,max in terms of a in in Fig. 16.
We find that the maximum eccentricity is very close to unity for a in > ∼ 0.01AU, but it decreases rapidly past a in ∼ 0.01AU. This is because GR effect suppresses the KL mechanism [69][70][71]. As a result, the case with a in = 0.01AU, which shows the shortest semi-major axis as well as the largest eccentricity, gives the loudest characteristic strain.
We also show the characteristic strain for Model BKL1 in Fig. 17. In this case, the strain is not large enough to give the observable value by DECIGO or BBO. We may need aLISA.
The plot of the frequency vs the characteristic strain gives helpful assessment of detectability provided high signal to noise ratio (SNR) of the waveform exists. Our hierarchical triple system may be one of the many source populations captured by future GW space detectors with following merger phase captured by ground detectors.
The area between the source and detector curves is related to SNR by the following relation [72], where h c (f ) is the strain amplitude of the source and h d (f ) denotes noise amplitude of detector. For the models shown in our analysis, the SNR for DECIGO is about 8 at e in,max ∼ 0.987 for Model AKL1. Similarly the SNR is about 16 for Model AKL2 at e in,max ∼ 0.964. Since the design sensitivity of BBO is better than DECIGO, the SNR for BBO would be larger. For Model BKL1, the signal is not loud enough for BBO, but that for aLISA SNR is ∼ 17, which could be observable.

VI. CONCLUDING REMARKS
We have shown features of gravitational waves from an inner binary in a hierarchical triple system with KL oscillations. The waveform changes its shape in time because of oscillation of the eccentricity. When the eccentricity is small, the waveform is a sinusoidal shape modulated by a tertiary companion, while it becomes one with sharp peaks near periastron point when the eccentricity gets large.
We have also examined the time variation of the characteristic strain curve, which may appear in sensitivity range of detectors when the eccentricity becomes large via the KL oscillations. This can be the first direct observation of KL oscillation which may further shed light on binary formation channels and surrounding environment of the binary system.
The interesting remaining questions are what happens when the system evolves beyond the KL stable region.
Because of the GW emission, the semi-major axis of the inner binary will decrease. Since other parameters such as m 1 , m 2 , m 3 and a out are conserved, the system will evolve vertically in Figs. 2 and 3 or horizontally in Figs. 4 and 5. Hence, we reach at the boundary of relativistic periastron shift. Beyond this boundary, KL oscillation may be suppressed, and then the effect of the tertiary component may become weak. As a result, we expect that GWs may not be observed before the coalescence of the inner binary. The question is whether the orbit of the inner binary is circularized via the GW emission at the coalescence or not. If the eccentricity remains then, we find the difference from a simple binary system. We may need numerical relativity or some other approach such as "effective one body system" to clarify it.
Another interesting issue is what happens when the LT effect becomes important. It can be the case if we consider SMBH. In Fig. 5, if the total mass of the inner binary is smaller than 10 3 M , when the semi-major axis decreases, the system will evolve into the region where the LT precession becomes important. In this region, we may find a chaotic KL oscillation [58], which may provide us strong emission of GWs. In order to study such a process, we have to include higher PN effects, i.e., at least 1.5 PN terms.
The work on these two issues is in progress. As for a triple system, we are also interested in more general situations, in which we expect chaotic behavior of a system. The GWs from such a chaotic system should also be studied. Since it may be difficult to make templates, we should find some other ways to clarify typical features of such a system in GW observations [73][74][75].
One of the most important issues about the GWs from a hierarchical triple system is the event rate. As for a formation of three body system, there are many works including numerical N -body simulation [18,20,22]. The event rate depends on many uncertain factors like star formation rate, binary fraction and distribution of initial parameters. Several papers discuss the increase in the merger rate of compact binaries around SMBHs due to KL oscillations [42,76,77]. The binaries in these triple systems undergo several KL cycles finally leading to the merger. In Ref. [42], the authors estimate the rate to be 0.56 Gpc −3 yr −1 . Similar analysis is done by Refs. [76,77], estimating the rate 1-3 Gpc −3 yr −1 . The difference is caused by different distance regimes from the SMBH and different initial conditions. For isolated triples, the merger rate to be 6 Gpc −3 yr −1 in the absence of natal kicks [78]. Globular clusters may harbor IMBHs at their centers. Ref. [79] studies the triple systems composed of the Stellar-mass binary with IMBH as a third body. The merger event rate induced by KL oscillations is estimated to be 0.06-0.46 Gpc −3 yr −1 .
We are interested in the models where merger timescale is much longer than Kozai timescale. Hence, if there are many KL oscillations before the merger, the signal about KL oscillations may lie in the observable range. In Ref. [47], using double-averaged equations, they analyze the fraction of KL binaries in two channels (isolated triples and galactic center) lying in the LISA frequency range. They show that of all merging binaries considered in their simulation, ∼ 39% of the isolated-triple channel, and ∼ 22% of the galactic-center channel, display significant KL oscillations in LISA frequency band.
However, the present models of a hierarchical triple system with KL oscillations, which gives observable GWs in an inspiral phase, may be difficult to be formed by evolution of a triple star system because both semi-major axes are very short. One possible way would be a capture process. The third body is bounded due to a close encounter with a BH binary in globular clusters or near galactic centers. The question is how likely such an interaction is. In [18,80], binary-single and binary-binary BH interactions are numerically simulated and several possible end states are discussed. In Ref. [80], they have shown that stable triple formation end state has higher cross section compared to some other possible end states.
We are planning to study in detail the formation process and probability to evaluate the event rate of compact hierarchical triples such that KL oscillations are not quenched by relativistic effects.

GWs from an isolated binary
First, we consider an isolated binary system for the illustration of waveform from a two-body configuration. In Fig. 18, we show gravitational waveform from circular binary with frequency twice the orbital frequency, 2f orb . When a binary has an elliptical orbit, GW signal is modulated due to different speeds in a single orbit. It results in 'spikes' due to higher emission near the periastron and hence higher amplitude as shown in Fig. 19. The energy spectra for Models ABC and ABE are shown in Fig. 20. As expected, the spectra from an isolated circular binary shows one peak corresponding to twice of the orbital frequency.
FIG. 20: The energy spectra for Model ABC (a circular binary) and ABE (an eccentric binary) with m1 = m2 = 10M and ain = 0.01AU. Many higher harmonics are produced because of high eccentricity.

GWs from a coplanar triple system
Here, we present the results from coplanar systems for both circular and eccentric inner binary, i.e., Models ACC and ACE (a coplanar circular binary in a hierarchical triple system with m 1 = m 2 = m 3 = 10M and a in = 0.01AU , a out = 0.1AU). The waveform from a circular binary in a coplanar hierarchical triple shows two harmonics corresponding to twice of inner and outer orbital frequency (see Fig. 21). In fact, just as the case of a circular isolated binary, we find a sinusoidal oscillation with one frequency 2f in , which is slightly modulated with another frequency 2f out because of the effect of the tertiary component. Fig. 22 shows waveform from an eccentric inner orbit with e in = 0.9 in a hierarchical coplanar system. In this case, the amplitude becomes large near the periastron point just as the case of an eccentric isolated binary, but the modulation by the tertiary component is also found. The energy spectra for Models ACC and ACE is shown in Fig. 23. Just as the case of an isolated circular binary, the spectrum from a circular inner binary in a coplanar hierarchical system shows two peaks corresponding to twice the orbital frequency of outer and inner orbit. Unlike circular orbits, many higher harmonics appear in the spectra from an eccentric inner binary in a coplanar hierarchical system. As a result, we find a broad band spectrum of GWs, which is a little different from the spectrum of the eccentric isolated binary. This is because the eccentricity in a hierarchcal triple system is oscillating, which gives a rich structure in a spectrum.