Post-Newtonian Kozai-Lidov Mechanism and its Effect on Cumulative Shift of Periastron Time of Binary Pulsar

We study the Kozai-Lidov mechanism in a hierarchical triple system in detail by the direct integration of the first-order post Newtonian equations of motion. We analyse a variety of models with a pulsar to evaluate the cumulative shift of the periastron time of a binary pulsar caused by the gravitational wave emission in a hierarchical triple system with Kozai-Lidov mechanism. We compare our results with those by the double-averaging method. The deviation in the eccentricity, even if small, is important in the evaluation of the emission of the gravitational waves. We also calculate the cumulative shift of the periastron time by using obtained osculating orbital elements. If Kozai-Lidov oscillations occur, the cumulative shift curve will bend differently from that of the isolated binary. If such a bending is detected through the radio observation, it will be the first indirect observation of gravitational waves from a triple system.


INTRODUCTION
Gravitational wave (GW) is one of the most interesting phenomena predicted by general relativity. It is the ripple on space-time caused by motions of massive objects like black holes. Orbital motions of close binaries emit GW which extracts orbital energy and gradually shrinks the orbit. The shrinking binary orbits can be observed through radio signals if the binary includes a pulsar as its component (Weisberg & Taylor 2005). Such a binary system with a pulsar is called a binary pulsar. A pulsar is a neutron star rotating fast and emitting radio signals with peaks whose period is quite precise. Due to this feature, it is possible to obtain various types of information from the observation of the radio signals from the pulsar; for example, we can know pulsar's rotational period, binary orbital period, and the information of binary orbital elements like semi-major axis and eccentricity (Smarr & Blandford 1976). Hence, if it is observed for a long term, the time evolution of orbital shape due to GW emission can be followed.
nary. The triple systems that can be divided into an inner binary and outer orbiting companion are called as hierarchical triple systems. Triple systems sometimes exhibit completely different orbital motions even if they have hierarchical structures. One of the most remarkable phenomena in hierarchical triple systems is the Kozai-Lidov (KL) mechanism (Kozai 1962;Lidov 1962). It is one of the most important orbital resonances that is mainly characterised by the secular changes of the eccentricity of the inner binary and the relative inclination between inner and outer orbits. These values oscillate exchanging their values with each other in secular timescale, that is, when the eccentricity increases, the inclination decreases, and vice versa, with timescale longer than both orbital periods. The eccentricity excitation in the inner binary is quite important for various astrophysical phenomena. For example, the large eccentricity can enhance GW emission in the binary and finally cause the merger of black holes (Blaes et al. 2002;Miller & Hamilton 2002;Liu & Lai 2017). In addition, the tidal force can also be enhanced with the excited eccentricity and the tidal disruptions of stars by supermassive black holes can be caused (Ivanov et al. 2005;Chen et al. 2009Chen et al. , 2011Wegg & Bode 2011;Li et al. 2005). In the context of planetary science, the formation of hot Jupiters (Naoz et al. 2012;Petrovich 2015;Anderson et al. 2016) or ultra-short-period planets (Oberst et al. 2017) are also said to be caused by KL mechanism. Recently, the GW emission from the hierarchical triple systems with KL mechanism has attracted attention of researchers. Some authors discussed about the waveform of GW from a binary in a hierarchical triple system and its observability (Hoang et al. 2019;Randall & Xianyu 2019;Gupta et al. 2020). If such systems exist and include pulsars as components of the binaries, the radio signal from the pulsar should also be detected. The CSPT curve described from the signal will tell how the third companion and GW emission affect the evolution of the binary.
In this paper, we first analyse the KL mechanism in relativistic systems in detail and compare the orbital evolution by the direct integration of the equations of motion and that by the well-known double-averaging method. We then investigate how CSPT curve changes with GW emission in hierarchical triple systems with KL mechanism. We treat general hierarchical triple systems in this paper, expanding the discussion in our previous letter (Suzuki et al. 2019), which treated only one example. If the CSPT curves predicted in this paper are detected through radio observation, it will be the first indirect observation of GW from a triple system. The paper is organised as follows: we summarise the important features of KL mechanism in §2. We describe our models in §3 and explain our methods in §4. The results and discussions are in §5. The conclusion follows in §6.

KOZAI-LIDOV MECHANISM
Hierarchical triple systems are three-body systems in which the motions of components can be divided into two Keplerian elliptic orbits called inner and outer orbits due to highly hierarchical configuration such that the outer semi-major axis is much longer than the inner one (see Fig. 1). We denote the masses of the components of inner binary by m1 and m2, and that of the tertiary companion by m3. Each orbit Figure 1. The hierarchical triple system is constructed from inner and outer binaries. The inner binary consists of objects whose masses are m 1 and m 2 , and the outer one is the pair of the inner binary and the third body with mass m 3 . The outer semi-major axis aout is much larger than the inner one a in .
in the hierarchical triple system is described with six orbital elements. In this paper, so called Kepler elements are used as the 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. It is well-known that these elements are constant in a twobody system, except the mean anomaly, which corresponds to the phase in an elliptic orbit. In the system that consists of three or more objects, in general, the trajectory of each component is not a closed elliptical orbit even in Newtonian dynamics. However, when the Hamiltonian of the total system is given by the sum of two-body Hamiltonians with perturbative interactions like a hierarchical triple system, each trajectory can be approximated by an elliptical orbit but its shape gradually changes in time. In such a case, the orbital elements of the osculating orbit, which is obtained by the instantaneous position and velocity, are used to describe the trajectory (see e.g. Murray & Dermott (2000)). In this paper, the osculating orbital elements of inner and outer orbits are represented with the subscripts 'in' and 'out', respectively. As for the outer orbit, we pursue the centre of mass of the inner binary rotating around the tertiary companion (see Fig. 1). Kozai-Lidov (KL) mechanism is one of the orbital resonances seen in hierarchical triple systems, which is discovered by Kozai (1962) and Lidov (1962) 1 . In the system where KL-mechanism occurs, the eccentricity of inner orbit ein and relative inclination I between inner and outer orbits oscillate in secular timescale. In this section, we shortly summarise some important features of KL-mechanism in Newtonian and post-Newtonian dynamics. The basic features of KL-mechanism are well described with quadrupolelevel approximation for a restricted triple system (see e.g. Shevchenko (2017)), in which one of the components of the inner binary is assumed as a test particle. We keep the lowest quadrupole order of the perturbed interaction terms in the Hamiltonian expanded in terms of the ratio of the semimajor axes. The detailed explanation of this treatment is given in Appendix A.
Not all of our models in this paper are the case of this restricted triple system. For example, some models have the inner binary constructed with two neutron stars. As shown in §4, we will not use the double averaging method in our analysis but we directly integrate the equations of motion. Hence the deviation from the test-particle limit is automatically taken into account. Here we just introduce the basic features of KL-mechanism obtained from the test-particle treatment in order to analyse our results. Note that the detailed analysis for non-restricted hierarchical triple system was given in Naoz et al. (2013a,b). In §5, we will revisit this point and will discuss the deviation seen in our simulation results from theoretical prediction with test-particle limit approximation.

KL oscillations in Newtonian Dynamics
First we discuss a restricted hierarchical triple system in Newtonian mechanics. KL-mechanism is an orbital resonance in hierarchical triple systems characterised by the oscillation of the eccentricity of inner orbit ein and the relative inclination I between inner and outer orbits on a secular timescale. We call this characteristic oscillation of ein and I as KL-oscillation. The relative inclination I is defined by cos I = cos iin cos iout + sin iin sin iout cos (Ωin − Ωout) . (1) The amplitude and timescale of KL-oscillation are determined by the conserved quantities in the restricted hierarchical triple system. From the quadrupole-order restricted triple treatment, two conserved quantities are obtained: CKL ≡ e 2 in 1 − 5 2 sin 2 I sin 2 ωin .
The detailed derivation is given in Appendix A1. When these values satisfy appropriate conditions, the KL-oscillation occurs. The KL-oscillations are classified into two types depending on the sign of CKL. When CKL ≥ 0, the condition for KL-oscillation is given as This KL-oscillation is called the "rotation" type because the periastron of the inner orbit rotates when the KL-oscillation proceeds, that is, the argument of periastron ωin increases monotonically. When CKL ≤ 0, on the other hand, the condition is described as This KL-oscillation is called the "libration" type because the argument of periastron ωin oscillates (librates) around π/2 or 3π/2 when the KL-oscillation proceeds. The ranges of conserved values (θ 2 , CKL) for both rotation and libration types are depicted in Fig. 1 in Antognini (2015). When the orbit is initially circular (CKL = 0 in Eq. (5)), the wellknown inclination range for KL-oscillation is obtained as The amplitude and timescale of the KL-oscillation depend on the type of oscillations even if the system size (masses and semi-major axes) is the same. For the amplitude, the difference is clearly seen in the exact formulae of maximum and minimum eccentricities shown in Appendix A1. The timescale of the KL-oscillation TKL is roughly estimated as where G is the gravitational constant and min = m1 + m2 is the total mass of the inner binary. This timescale depends only on the system size and the eccentricity of the outer orbit, but the exact oscillation period also depends on the conserved quantities of the system (see Appendix A1 for the reason). In §5, we confirm it by comparing our simulation results with different conserved quantities.

Post-Newtonian Correction
In the restricted hierarchical triple system with quadrupolelevel approximation, the GR correction is usually discussed by adding a simple correction term to the perturbation potential, which is derived by double-averaging of the first order post-Newtonian (1PN) Hamiltonian of two-body relative motion (the detail is given in Appendix A2). Note that Will (2014a,b) pointed out that this approach for the GR corrections is not always appropriate. Strictly speaking, for secular calculation due to the risk of the violation of energy conservation, we have to consider the effect of "cross terms" between the Newtonian perturbations and the post-Newtonian precession effect. In this section, however, we consider the GR correction without cross terms for interpretation of our numerical results (see also Appendix A2). In our simulation, as shown in §4, the equations of motion are directly integrated. Hence the effect of the cross terms is automatically taken into account. The restricted triple systems with the GR correction have two conserved values as in the Newtonian dynamics. θ does not change from Newtonian one, but CKL is modified as which is a dimensionless constant describing the strength of GR effect with rg,in = Gmin/c 2 . Note that C (GR) KL is the same as CKL for circular orbit.
The classification conditions of KL-oscillations are C (GR) KL ≥ 0 for "rotation" type while C (GR) KL ≤ 0 for the "libration" type, respectively. The amplitude and timescale of KL-oscillation with the GR correction depend on the conserved quantities and vary from those in Newtonian analysis. In §5, we compare the Newtonian and GR results.

MODELS
We study GW emission effects on CSPT (cumulative shift of periastron time) of binary pulsars in hierarchical triple systems with the KL-oscillations. As discussed in our previous letter paper (Suzuki et al. 2019), this effect could be found in long-time observation of radio pulses from the pulsar. We have shown only one model with initially circular inner binary as an example. In this paper, we analyse a broad range of parameters. We first obtain constraints on parameters by imposing stability of the system and observable timescale and we then analyse several models in the allowed parameter range. Before discussing the constraints, we first classify hierarchical triple systems into three classes according to their mass ratio: , we may also see the KL-oscillations (Blaes et al. 2002;Wen 2003;Thompson 2011;Liu & Lai 2018) as in Class [1] as long as aout ≫ ain. If aout is not large enough as compared to ain, such a system does not have a sufficient hierarchy and then the interaction between the inner and outer orbits becomes strong. As a result, both orbital elements will change extremely with time and the orbit will become chaotic. It may become unstable.
In Class [3], when the outer object can be treated as a test particle (aout ≫ ain), the inner orbit is not affected so much by the tertiary object, while the orbital elements of the outer orbit may change with time. However, it is known that the eccentricity of the outer orbit does not change with time at least in the quadrupole order approximation. Instead we may expect the oscillation between the relative inclination I and the longitude of ascending node of the outer orbit Ωout in secular timescale (Naoz et al. 2017). Since we are interested in CSPT with the KL oscillations, i.e., CSPT via the time change of the pulsar's eccentricity, we discuss only Classes [1] and [2].
In order to see CSPT through radio signals, each model should contain a pulsar as a component of the inner binary. As a companion of the pulsar in the inner binary, in order to find large GW emissions from the inner binary and to neglect the tidal dissipation effect, we may choose a compact object with a similar or larger mass than that of the pulsar, i.e., a neutron star (NS) or a black hole (BH). If the companion is a non-compact object like a main sequence star, a strong tidal force from the pulsar deforms the companion star and . In the blue thin-stripe region, a hierarchical triple system is stable. The condition for the KL-oscillations not to be suppressed by the post-Newtonian relativistic effect is given by the magentacoloured region. The overlapped region gives a stable KL oscillations. The dark-green lines show the timescales of KL-oscillations (T KL = 10, 10 2 , and 10 3 yrs), which should be shorter than our lifetime (< 100 yrs) for observation. Our models given in Table 1 are shown by the black dots.
the orbital energy is dissipated by friction in the star. Since such dissipation by the tidal force may affect the periastron shift in addition to the GW emission, CSPT becomes more complicated, which is beyond the scope of this paper. Hence we analyse three types of model for inner binaries; P-NS binary (pulsar + NS), P-BH binary (pulsar + BH) P-IMBH binary (pulsar + intermediate mass BH) . m1 and m2 are the masses of the companion and pulsar in the inner binary, respectively. We choose those concrete values given in Table 1. There exist some conditions for the parameters of the outer orbit in order for the inner binary to exhibit stable KLoscillations. We show those constraints in Fig. 2-4 in terms of the semi-major axis of the outer orbit aout and the mass of the third body m3 by fixing parameters of the inner binary. The dashed black line shows the constraint for the outer binary mass m3, which should almost be the same or larger than the mass of the inner binary min. The second condition is stability of the hierarchical triple systems, i.e., the socalled "chaotic boundary". As given in Mardling & Aarseth (2001), the following condition should be satisfied so that the hierarchical structure of the system does not break at least in the initial state: The stability condition (11) is shown by the blue thin-stripe region. The third condition is given by Eq. (10), which ensures KL-oscillation occurs even in a relativistic system. We  (0) and aout(0) are the initial values of the semi-major axes of the inner and outer orbits, respectively. ǫ (1PN) is the strength of the relativistic effect defined by Eq. (9) for a restricted hierarchical triple system. P, NS, BH, IMBH and SMBH mean a pulsar, neutron star, black hole, intermediate mass black hole and supermassive black hole, respectively.  depict this condition by setting ein = eout = 0 because it does not change so much even for non-zero eccentricities. This relativistic constraint is given by the magentacoloured region. In order to observe the effect of KLoscillation on CSPT, the timescale of KL-oscillation should be short enough, compared with our lifetime. As mentioned in §2.1, the timescale of KL-oscillation is roughly estimated by Eq. (7). We show some contour lines of TKL by the darkgreen lines (TKL = 10, 10 2 and 10 3 years). When the tertiary companion has the parameters both in the blue thin-stripe region and the magenta-coloured region in Figs. 2-4, the KL-oscillation will occur with appropriate timescale. We also show our model parameters by the black dots with the model names in Fig. 2-4. We analyse nine models: for P-NS inner binary, we discuss four models; PNN, PNB, PNIB and PNSB, in which the tertiary companion is a neutron star (NS), black hole (BH), intermediate mass black hole (IMBH), and supermassive black hole (SMBH), respectively. For P-BH inner binary, we consider three cases: PBB, PBIB and PBSB, in which the tertiary companion is a BH, IMBH and SMBH, respectively. We also analyse the model PIBIB and PIBSB ; both systems have a P-IMBH inner binary, and an IMBH or SMBH as a tertiary companion. We choose the masses of a pulsar (or NS), BH, IMBH and SMBH as 1.4M⊙, 30M⊙, 10 3 M⊙ and 10 6 M⊙, respectively. The model parameters are summarised in Table 1.
Here we remark the Lense-Thirring precession effect. This is one of the spin-orbit coupling effects appearing in 1.5 post-Newtonian order correction (Barker & O'Connell 1975). Recent studies have shown that the Lense-Thirring precession caused by the rapid rotation of an outer supermassive black hole in a hierarchical triple system changes the evolution of the KL-oscillation Liu et al. 2019). As in Liu et al. (2019), TLT is evaluated by where χ3 ≤ 1 is the rotation parameter of the third object in the hierarchical triple system. By using Eq. (7), TLT ≫ TKL gives the condition to neglect the Lense-Thirring effect, i.e., ain au We imposed χ3 = 1 in above estimation. Since all models in Table 1 satisfy this condition, we can neglect the Lense-Thirring effect in our calculation.

BASIC EQUATIONS
For the models explained in §3, we directly integrate the equations of motion for their orbital evolution. Then we analyse the behaviour of KL oscillations and evaluate the cumulative shift of periastron time (CSPT) of the inner binary.

Equations of motion for three body system
In order to solve relativistic motions of our three-body system composed of compact objects, we use the first-order post-Newtonian equations of motion, which are called as the Einstein-Infeld-Hoffmann (EIH) equations (Einstein et al. (1938)): where m k , v k , x k are the mass, velocity and position of the k-th component of the system (k = 1, 2 and 3). Note that this equation could be derived from the Lagrangian given by Lorentz & Droste (1917). In our study, Eq. (14) is numerically integrated by using the 6-th order implicit Runge-Kutta method. The coefficients of 6th-order Runge-Kutta are obtained from Butcher (1964). The back reaction of GW emission to the orbital evolution can be treated by including the 2.5 order post-Newtonian terms. However, since the back reaction in a few KL-oscillation timescale is so small, it does not change our result. Hence we consider only the first order of the post-Newtonian equations for the orbital evolution.

Initial Conditions
In order to set initial conditions for our simulation, we not only need the semi-major axis a but also other parameters like the eccentricity e and the inclination i. These parameters fix the conserved quantities θ and C GR KL , which classify the type of KL-oscillation as "libration" or "rotation" (see §2). Hence we prepare four sets of initial parameters named as "initially circular libration (ICL)", "initially circular rotation (ICR)", "initially eccentric libration (IEL) " and "initially eccentric rotation (IER)". For "Initially circular", we set ein = 0.01, while for "initially eccentric" we choose ein = 0.6. The other parameters are determined to find C GR KL < 0 for libration and C GR KL > 0 for rotation. The parameters of each type are summarised in Table 2 and are used for post-Newtonian calculations.
To study the relativistic effect, we also perform the Newtonian calculation. We choose two conserved quantities as CKL = C (GR) KL and the same value of θ 2 as the post-Newtonian one, which are obtained by setting the initial periastron argument as ωin given in the last column in Table 2.
These initial orbital elements are converted into the position and velocity vectors, x k and v k , in Cartesian coordinates, whose origin is the centre of mass of whole system. The x-y plane of our coordinate system is chosen to be the initial outer orbital plane. The detailed conversion formula is given in Appendix B1 (See also e.g. Murray & Dermott (2000)). By using Cartesian initial variables, the above EIH equations (14) are integrated numerically and the osculating orbital elements are evaluated at each time step. The procedure to evaluate orbital elements from positions and velocities at each time step is also explained in Appendix B2.
The integrated inner orbit is not exactly a closed ellipse, but it fluctuates with small amplitudes because of the effect of the tertiary component. As a result, the orbital parameters of the osculating orbit evaluated at each step are oscillating, which seem to be artificial. Hence we take an average of these elements for each inner cycle to extract the effective values at each cycle. We describe such averaged orbital elements with a bar, e.g.,āin andēin. Those elements evolve in secular timescale due to the effect of the third body.

Cumulative Shift of Periastron Time (CSPT)
The orbital energy of inner binary, if it is close enough, is extracted little by little via the GW emission. The energy dissipation makes the semi-major axis of the orbit shrink and then the period of the orbit becomes shorter and shorter. As derived in Peters & Mathews (1963), the period change for each orbital cycle iṡ where Pin is the orbital period of the inner binary given by When the energy dissipation is evaluated for one binary cycle, the orbital elements can be treated as constant because the back reaction of energy dissipation is small enough in such a timescale. Here we use the averaged values,ē andā, instead of the osculating orbital elements, e and a, to reflect the effective shape of the orbit for one cycle. Whenē and a evolve with secular timescale such as the KL-oscillation timescale,Ṗin also changes with time. This period shift can be seen by observing the cumulative shift of periastron time (CSPT) through radio signals  Table 2. The important parameters in initial conditions for KL-oscillations for post-Newtonian calculations. We analyze four sets of initial parameters; " initially circular libration" (ICL), "initially circular rotation" (ICR), "initially eccentric libration" (IEL) and "initially eccentric rotation (IER). e, i, ω are the eccentricity, the inclination, and the argument of the periastron, respectively. We also show two conserved quantities, C (GR) KL and θ 2 , in post-Newtonian dynamics. For "Initially circular", we set e in = 0.01, while for "initially eccentric" we choose e in = 0.6. The other parameters are determined to find C (GR) KL < 0 for libration and C (GR) KL > 0 for rotation. For the outer orbit, eout = 0 and iout = 0 • are used and ωout cannot be defined. About the parameters other than those shown in the table, the longitude of the ascending node Ω is set as 0 for both inner and outer orbits, and the mean anomaly M is set as 0 • and 20 • for inner and outer orbits. To study the relativistic effect, we also perform the Newtonian calculation. We choose two conserved quantities as and the same value of θ 2 as the post-Newtonian one, which are obtained by setting the initial periastron argument as ω in given in the last column.
from a binary pulsar just as the observation of the Hulse-Taylor binary (Weisberg & Taylor 2005). In this paper, we expand the analysis to hierarchical three-body systems. The CSPT of the inner binary ∆P is defined as where TN is the N -th periastron passage time and Pin (0) is the initial orbital period of the inner binary. From the definition of TN , we obtain where Pin(t) is the binary period at time t, which changes in time by the GW emission as By substituting Eqs. (18) and (19) into Eq. (17), the CSPT ∆P is described as Since the emission energy of GWs is quite small, we usually expect In fact, for Hulse-Taylor binary pulsar (Weisberg & Taylor 2005), since we have the condition (21) is true if t ≪ 3.7 × 10 8 yrs. Hence, when we are interested in the time-scale such that TN ≪ 10 8 yrs, we approximate ∆P as Note that if we assumeṖin(t) is almost constant, that is, Pin(t) ≈Ṗin(0), ∆P is given by which was used in Weisberg & Taylor (2005). However, in a hierarchical triple system with the KL oscillation,Ṗin(t) is not constant but may change in time with the KL-oscillation timescale. Hence, in this study, we evaluate ∆P by Eq. (24) with Eq. (15).
Our analysis can be applied to a general stable threebody (or N -body) system with a binary pulsar as long as the condition Eq. (21) is satisfied. Here we stress that the CSPT could be observed through radio signals from the pulsar as the accumulated effect. Highly accurate observation of radio pulsars enables us to see this CSPT even for such weak GW emission that the back reaction of GW emission on the orbital elements is negligibly small. The CSPT observation through the radio signals from a binary pulsar in a triple system may be the precursor of detection of gravitational waves from a triple system with the KL-oscillation (Gupta et al. 2020).

Orbital Evolutions
In our simulation results, the stable orbital evolutions are observed in all the models shown in Table 1. We show the results of PNN model and PNIB model as representative.
The mass hierarchy in PNN model is the smallest in all the models and it is expected that the deviation from testparticle approximation used in §2 is the largest. In PNIB model, on the other hand, ǫ (1PN) is the second largest as seen in Table 1 and relativistic effect in this model may become important.

Evolution of Orbital Parameters
Figs. 5 and 6 show the time evolution of the averaged inner eccentricityēin, relative inclinationĪ and KL-conserved valueθ 2 of the PNN model. Fig. 5 shows the result of libration type KL-oscillations, while Fig. 6 exhibits those of rotation type KL-oscillations. In each figure, top and bottom panels correspond to the results of the initially circular and eccentric cases, respectively, whose parameters are given in Table 2. In Figs. 5 and 6, the KL-oscillation is observed in all panels with different amplitudes and timescales: The initially eccentric cases (bottom panels) have smaller amplitude and shorter timescale than those of initially circular cases (top panels). This is because in the initially circular case, the eccentricity oscillates between zero and some finite value, while in the initially eccentric case, it oscillates between two finite values around the initial value. The same behaviour is also found from the figures of the eccentricity in the double-averaging method. CKL (or C KL ) is very small in the initially circular case, while it is not so small in the initially eccentric case (see Fig. A1 in Appendix A for the Newtonian case).
However, in all panels of both figures, θ 2 is not exactly constant but oscillates with the same period as that of the KL-oscillation although it should be constant in the analysis with test-particle quadrupole approximation in §2. This is because all masses of the components in the system are the same in PNN model and the hierarchy assumed in §2 is not enough in this model, that is, the test-particle approximation does not work exactly in this model. This small deviation from the test-particle limit is consistent with the discussion given in Naoz et al.  the quadrupole approximation 3 . Each panel in Figs. 7 and 8 is the same evolution as that shown in the corresponding panel in Figs. 5 and 6. We find the difference between two oscillation curves in all panels. The timescale of KLoscillation obtained from direct integration is smaller than that calculated in double-averaging method. The deviation in timescale is much more obvious in initially eccentric results (bottom panels of Figs. 7 and 8). For the amplitude, the tendency of the difference is not the same in all panels. In the results of ICL and ICR types (top panels in Fig. 7 and  8), the amplitude of KL-oscillation is smaller in our direct simulation than that obtained with double-averaged calculation. Both curves in these panels have the same minimum values but the maximum values are enhanced in light-green lines. In the result of IER type (bottom panel of Fig. 8), the enhancement of the amplitude in the double-averaging method is observed as seen in ICL and ICR types, but both maximum and minimum values in light-green line are a little different from those of dark-green line: the maximum value is larger and the minimum value is smaller in lightgreen line. The result of IEL type (top panel of Fig. 7), on the other hand, has an opposite feature from the the other three types: The amplitude obtained from double-averaged calculation is smaller than that obtained from direct simulation. In this case, the maximum values in both curves are the same but the minimum value in the dashed lightgreen curve is larger than that of solid dark-green curve. We remark that those differences between eccentricity evolution obtained from direct integration and that by doubleaveraging method may be crucial when we evaluate the GW emission for the systems with finite masses, that is, one may overestimate or underestimate the maximum or minimum value of the eccentricity when we use the double-averaging method. The amplitude and frequency of the gravitational waves are strongly sensitive to the eccentricity, especially for the highly eccentric orbit like e > 0.9. It may be important to calculate the evolution of such an orbit by direct integration 4 .

Evolution of Orbital Parameters
Figs. 9 and 10 are the same figures as Figs. 5 and 6 but for PNIB model. Figs. 9 and 10 reflect the features of libration and rotation types of KL-oscillations, respectively. In each figure, top and bottom panels are the results of initially circular and eccentric types, respectively. As seen in Figs. 5 and 6, initially eccentric cases (bottom panel) have smaller amplitude and shorter timescale than those of initially circular cases (top panels), which is similar to the PNN model. As for the KL oscillation period, it does not seem to depend on the oscillation types in the initially circular case, while in the initially eccentric case, the rotation type (the bottom   Table 3. Comparison between the results by the direct integration and those by the double averaging method for the PNN model with ǫ (1PN) = 0.177. We show the maximum and minimum eccentricities, ∆e = emax − e min , which gives the oscillation amplitude, and the KL oscillation period T KL . The first rows give the results by the direct integration, while the second rows show the results by the double averaging method. panel in Fig. 10) gives shorter oscillation time than that in the libration type (the bottom panel in Fig. 9). θ 2 is almost constant in PNIB model unlike that in PNN model. It is because the test-particle approximation is valid in PNIB model. In fact, the deviation from the double-averaging method is smaller than that of PNN model.
Newtonian v.s. post-Newtonian PNIB model has the second largest value of ǫ (1PN) in Table 1 and its relativistic effect is the strongest in our models except PIBIB model. Since the main features are the same in both models, we shall discuss the PNIB model as a representative of relativistic ones. In Figs. 11 and 12, we show the evolution of the eccentricities obtained by Newtonian and post-Newtonian direct simulations. Each figure exhibits the results of libration and rotation types of KL-oscillations. The top and bottom panels in each figure correspond to the results of initially circular and eccentric types. The Newtonian and post-Newtonian results are described by the light-and dark-green curves, respectively. The tendency of the difference between two curves is not the same in all panels. In the results of ICL and ICR types (top panels of Figs. 11 and 12), the amplitude of KL-oscillation is smaller in post-Newtonian simulation than that obtained from Newtonian calculation. Both curves in those results have the same minimum values (about zero), but the maximum value is suppressed in post-Newtonian curve. The KL-timescale is a little longer in post-Newtonian result in the initially circular types. In the results of IEL and IER types (bottom panels of Figs. 11 and 12), on the other hand, the KL-timescale obtained in post-Newtonian calculation is shorter than that obtained from Newtonian one. Interestingly, IEL (bottom panel of Fig. 11) and IER (bottom panel of Fig. 12) have different features in the amplitude. In the result of IEL type, the amplitude obtained by post-Newtonian simulation are smaller than those of Newtonian result; unlike results of the ICL and ICR types, both maximum and minimum values are suppressed in this case. Newtonian simulation post-Newtonian simulation Figure 11.
Comparison between Newtonian and post-Newtonian evolution curves of the averaged inner eccentricityē in for the "liblation" type of KL oscillations in the PNIB model. Top and bottom panels correspond to the results of ICL and IEL types, respectively. The light-and dark-green curves describe the results obtained from Newtonian and post-Newtonian direct simulations.
Type On the other hand, in the result of IER type, the amplitudes of Newtonian and post-Newtonian results are almost the same but both maximum and minimum values of post-Newtonian result are shifted downward. These complicated features can be understood basically by using the double-averaging method, which is given in Appendix A2. As shown in Fig. A4, the curve of the maximumminimum eccentricity in terms of CKL in Newtonian dynamics is shifted to the right when the post-Newtonian correction term is taken into account. Here we have used C (GR) KL instead of CKL as the horizontal axis because it is conserved and classifies the oscillation types, libration or rotation. Hence when we include the post-Newtonian correction term, fixing two conserved quantities (θ 2 and C (GR) KL = CKL), we find that the maximum value decreases and the minimum value increases for the libration type, while both maximum and minimum values decrease for the rotation type. As for the KL oscillation, the analysis by the double-averaging method explains the results by the direct integration (compare Figs. 11 and 12 with Table A1. ).

Irregularity of KL oscillation period
As we showed above, the amplitude of KL oscillation and its period can be understood basically by the double-averaging method. However we find that there appears an irregularity of the period in some models. For example, the KL-oscillations in ICR type of the PNB and PBB models show irregular periods (see Fig. 13). This irregular behaviour of the KL-oscillation period was already found in Antonini & Perets (2012). They calculated orbital evolutions of BH binaries around SMBH by using N-body integrator and found the irregular periods and amplitudes in the KL-oscillation (Fig. 3 in their paper) .
Since the calculations in Antonini & Perets (2012) and ours are performed by the direct integration, one may naturally expect some deviation from the double averaging method, in which the KL-oscillation period is regular. However, since the deviation in our calculation is very small, the double averaging method may provide almost correct results.   3.2 × 10 −5 5.6 × 10 −5 3.1 × 10 −5 −1.1 × 10 −4 9.34 yrs 8.88 yrs 7.94 yrs 7.09 yrs Table 5. The period of KL oscillations. The period n (n = 1, 2, 3) denotes the period from the n-th peak of the eccentricity to the (n + 1)-th peak. C KL /C (GR) KL is the "conserved" value after the eccentricity passes through the maximum value. The periods calculated by the double-averaging method with the same values of C KL /C (GR) KL are given in the third rows of each period.
Note that the amplitude and timescale of KL-oscillation are strongly dependent on two conserved quantities θ 2 and C (GR) KL , but not so much on ǫ (1PN) except for the ICL oscillation type, in which the relativistic effect is large because it changes the existence range of KL oscillation. Hence we analyse the behaviour of the "conserved" quantities in our simulations. As for θ 2 , although it oscillates with the outer orbit period, the averaged value is almost constant except around the time when the eccentricity reaches the maximum value. We then show the time evolution of CKL and C is not conserved when the eccentricity reaches the maximum value. However it becomes almost constant again when the eccentricity decreases.
In order to see the detail, in Table 5, we show the numerical values of the oscillation periods. The period n (n = 1, 2, 3) denotes the period from the n-th peak of the eccentricity to the (n+1)-th peak. We also show the constant "conserved" values after the eccentricity passes through the maximum value in Table 5. We evaluate the KL oscillation periods by the double-averaging method with those values of CKL/C (GR) KL , which are given in the third row of each period in Table 5. We find that those periods are consistent with the numerical ones by the direct integration. We believe that these small deviations of the "conserved" values in each period causes small irregularity of the KL oscillation period. We still have a small difference from the numerical simulation, which may be because of large deviation of CKL/C (GR) KL near the maximum eccentricity.

Cumulative Shift of Periastron Time (CSPT)
The KL-oscillations shown in §5.1 affect the evolution of the CSPT ∆P of binary pulsar in the hierarchical triple system. As we showed in the previous letter Suzuki et al. (2019), if a hierarchical triple system shows the KL oscillations in ob- KL-ICL isolated binary KL-IEL isolated binary Figure 15. The CSPT curve for libration type of PNIB model is shown. Top and bottom panels are the results integrated from the time of maximum and minimum eccentricities, respectively. The blue and red solid curves correspond to ICL and IEL types, respectively. The dashed curves are those of isolated binaries whose parameters are the same as the initial values of the inner binaries of corresponding types. servation period, we expect the bending of CSPT curve. It is because when the eccentricity becomes large, the amount of GW emission increases, and then the change of orbital period gets large. Here we shall discuss how the bending of CSPT curve depends on the models or types of KL oscillations.
For each model in Table 1, we have calculated the timeevolution of CSPT as explained in §4. Since the behaviour of the CSPT curve does not depend so much on the models except for the timescales, we show the results only for the PNIB model. Figs. 15 and 16 show the results of libration and rotation types KL-oscillations, respectively. In each panel, the red and blue solid curves show the results of initially circular and eccentric types, respectively. The top panels show the CSPT curves calculated from the time when the maximum eccentricity is found in each KL-oscillation type (at t = 15.21yr, t = 0yr, t = 14.96 yr and t = 1.32 yr for ICL, IEL, ICR and IER types, respectively), while the bottom panels exhibit those calculated from the time when minimum eccentricity is reached (at t = 0yr, t = 7.92yr, t = 0 yr and t = 6.10 yr for ICL, IEL, ICR and IER types, respectively). It shows that the CSPT curves become completely different depending on the choice of the initial time of integration TN = 0 even for the same model. For reference, we also show the CSPT curves of the isolated binary whose parameters are the same as the initial parameters of the inner binary in corresponding hierarchical triple models, by the red and blue dashed curves. The CSPT curves of isolated binaries are approximated by the quadratic functions as Eq. (25). At first the CSPT curves of KL triple system coincide with the quadratic curves of corresponding isolated binaries, but when the eccentricity changes with KL-mechanism, the curves of the triple-system bend and the discrepancy from the binary curves becomes large as already shown in Suzuki et al. (2019). This is because the period change of the inner binary due to GW emission (Ṗin) depends on the orbital eccentricity as given by Eq. (15). Hence when the orbital eccentricity changes,Ṗin also changes, and then the CSPT curve deviates largely from the quadratic curve.
In the top panels of Figs. 15 and 16, the solid curves at first coincide with the quadratic curves with eccentric orbits, but they switch to the less steeper curves as the eccentricities become smaller by KL-mechanism. This feature results in the slower decrease of ∆P in the triple system compared with that of the isolated eccentric binary. The slope and bending timescale of red and blue solid curves are different from each other depending on the amplitude and KL-timescale. While, in the bottom panels of Figs. 15 and 16, the switch from the circular curves to the eccentric steeper curves causes rapid decrease of ∆P in the triple system curves than those of isolated circular binaries.
This bending feature may be useful to see KL-oscillation from pulsar observation. The shape of the CSPT curve has the information of the eccentricity and the KL-oscillation timescale in its slope change. The bending of the CSPT curve is clear when the curve is integrated from minimum eccentricity, but the curve from the maximum eccentricity does not show clear bending. However, the change of the CSPT curve becomes clearer if the time-derivative of ∆P is In each panel, the red and blue curves show the results of initially circular and eccentric types. The top and bottom panels in those figures show d∆P /dTN curves calculated from the time when the maximum and minimum eccentricities are obtained, respectively. By comparing the values of slope and bending timescale obtained from these figures with those of the isolated binary, the shape of eccentricity oscillation curve and corresponding parameters can also be estimated. It has already been pointed out that the KL-oscillation should be observed through the long-period radio observation of the orbital elements of the binary pulsar (Gopakumar et al. 2009;Zwart et al. 2011). In real observation, however, the observational data is sometimes missed due to some reasons; for example, in the observation of the Hulse-Taylor binary, the data was not obtained for a decade of 1990s because of the major upgrades of Arecibo telescope (Hulse 1994). If this unseen period is completely overlapped with the time when eccentricity is changed from the initial value with KL-oscillation, it is difficult to recognise whether KL-oscillation occurs or not only from orbital element data. Even in such case, we can conclude that KL-oscillation occurs in the system if the CSPT curve deviates from that of isolated binary in late phase.
Some readers may worry about the spin evolution of the pulsar caused by the spin-orbit coupling in 1.5 order post-Newtonian terms (Barker & O'Connell 1975) because it may change the direction of the pulsar rotation axis and affect the radio observation, that is, the change of beaming direction of pulse signal may cause the disappearance of the pulsar. Following Liu & Lai (2017, the evolution of spin in relativistic KL-oscillation can be characterised with the "adiabaticity parameter" A defined as the ratio of the de Sitter spin precession rate ΩSL to the orbital precession rate by KL-oscillation ΩL. The adiabaticity parameter A is described as where µin = m1m2/(m1 + m2) is the reduced mass of the inner binary. This parameter is quite similar to ǫ (1PN) defined as Eq. (9). Hence for the system with KL-oscillation, which satisfies the condition Eq. (10), we find that the adiabaticity parameter A satisfies The adiabaticity parameters of our models are summarised in Table 6. In case of A ≪ 1, the spin evolution is classified as "nonadiabatic", that is, the orbital precession by KL-oscillation is much faster than the relativistic spin precession, and then the spin axis cannot 'catch up' with the precession of angular-momentum axis. In such a situation, the spin axis of the pulsar is expected to be parallelly transported just as in the Newtonian case, and then the beaming direction of the radio signal is expected not to change so much even when the inclination changes by KL-oscillation. The PBB model corresponds to this case. For the other models, A is still smaller than unity, but not so much. The spin axis of the pulsar in the system with such mid-range of A is perturbed around its initial direction as shown in Liu & Lai (2018). If the perturbation of the spin direction is large enough so that the beaming angle of the pulsar goes out from the observable range, the radio signal from the pulsar will disappear and will rarely re-appear due to its complicated evolution. If the disappearance of a pulsar in triple system is observed, it will be an important example of the 1.5 post-Newtonian effect on the KL-oscillation. The critical value of A that causes the disappearance of the signal should depend on the emission mechanism of the pulsar, the intensity of the radio signal, the distance to the system, and the opening angle of the radio telescope. If the CSPT is observed for a whole period of KL-oscillation despite the precession of the spin direction of the pulsar, it means the pulsar is successively observed from some different directions and such observation may give new information about the pulsar. The bending of the CSPT curve may be a rare event because such compact hierarchical three-body systems with high inclination may need to be formed by the dynamical interaction in dense environments like the globular clusters and the galactic nuclei (Kulkarni et al. 1993;Samsing et al. 2014;Zevin et al. 2019). However, as discussed in Suzuki et al. (2019), this interesting signal is important not only to confirm the existence of the third body but to provide a first indirect evidence of GW emission from the triple system with KL-oscillations. GW emission makes the inner binary more compact and GW waveform from such compact triple system with KL-oscillation can be observed by future GW detectors (Gupta et al. 2020) like LISA (Amaro-Seoane et al. 2017), DECIGO (Sato et al. 2017), and Big Bang Observer (Harry et al. 2006).

CONCLUSIONS
In this paper, taking the 1st post-Newtonian relativistic correction into account, we have studied the KL-oscillations in hierarchical triple systems with a pulsar and calculated the cumulative shift of periastron time (CSPT). The KLmechanism is one of the orbital resonances that appear in the hierarchical triple systems characterised as the exchanging oscillation with the inner eccentricity and the relative inclination. When the eccentricity of the binary pulsar is excited by KL-oscillation, it enhances GW emission from the binary and it changes the shape of the CSPT curve. We have analysed the KL-oscillations in several models with a pulsar, and those effects on the CSPT curves.
We have first analysed the KL-oscillations for the models with different initial parameters. We have classified those models into four types (ICL, IEL, ICR, and IER). We have calculated their orbital evolutions by the direct integration of 1st post-Newtonian equations of motion. The four KLtypes have different amplitudes and timescales and, in addition, the non-test particle limit effect and the relativistic effect appear differently. In the result of the model with weak mass hierarchy (e.g. PNN model), we find that KL-"conserved" value θ 2 is not conserved but oscillating whereas it should be constant in double-averaged method with testparticle limit approximation. It has also been found that the amplitudes and timescales obtained in direct integration do not coincide with those in double-averaged method. The tendency of these discrepancies is different in the four types of KL-oscillations. The amplitudes and frequencies of the emitted gravitational waves are quite sensitive to the eccentricity, and these differences between eccentricity evolution in direct integration and that obtained from double-averaged method may be crucial when we evaluate the GW emission for the systems with finite masses, that is, one may overestimate or underestimate the maximum or minimum value of the eccentricity when we use the double-averaged method.
In the model with large ǫ (1PN) (e.g. PNIB model), we could observe clear differences between the results obtained by Newtonian and post-Newtonian direct integrations. The post-Newtonian effects appear differently in the four types of KL-oscillations. The complicated behaviours can be understood theoretically by using the double-averaging method with 1st-order post-Newtonian corrections. However, in some models (e.g. PNB and PBB models), we have observed KL-oscillation with irregular periods, which cannot be explained by double-averaging method with quadrupole-order approximation. This may be because the KL-conserved quantities are not exactly constant in the direct integration.
The KL-oscillation effect appears in the CSPT curve as the bending of the curve. The slope of the curve at each phase reflects the maximum or minimum eccentricity and the time between two bending points corresponds to the timescale of KL-oscillation. The CSPT curves become completely different depending on the choice of the initial time of integration even for the same model. The bending of the CSPT curve is clear when the curve is integrated from minimum eccentricity, but the curve from the maximum eccentricity does not show clear bending. In such case, the time derivative of the CSPT can be a good indicator for the bending of the CSPT curve.
The system that causes this interesting signal may be rare because such compact hierarchical triple systems with high inclination need to be formed by dynamical interaction in a dense environment like a globular cluster or the galactic center. However, once such systems are observed with the pulsar signal, it is very important because it is the first indirect observation of GW from triple systems. In addition, it will be the precursor of the direct detection of the waveform by the future gravitational detectors like LISA, DECIGO and Big Bang Observer. Some highly relativistic triple systems should show the spin precession of the pulsar caused by the 1.5 post-Newtonian effect from the outer orbit and it will change the beaming angle of the pulsar. If the beaming angle of the pulsar is perturbed and goes out of the observable range, the radio signal from the pulsar will disappear and rarely appear again. The disappearance of the signal from a pulsar in triple system will provide one of the important examples of the 1.5 post-Newtonian effect on the KL-oscillation. On the other hand, if the CSPT is observed for a whole period of KL-oscillation despite the precession of the spin direction of the pulsar, it corresponds to the successive observation of a pulsar from different directions and such observation may give new information about a pulsar.

A1 Newtonian Dynamics
Here we discuss the restricted hierarchical triple system. We choose our reference plane to define the inclinations as the initial orbital plane of the outer orbit. Since the outer inclination is conserved in the restricted triple system, we find that iout = 0 and then the inner inclination iin is the same as the relative inclination I between inner and outer orbits 1 .
The secular time evolution of the osculating orbital elements of the inner orbit is described by the Lagrange planetary equations, which is decoupled from the orbital motion of the outer orbit in the restricted hierarchical triple system as where n is the mean motion of the inner orbit, which is defined by and VS is the double-averaged perturbation potential in the Hamiltonian of the motion of a test-particle in the triple system. "Double-averaged" means that the corresponding term is averaged for both periods of inner and outer orbits. In this section, we drop the subscript "in" for the inner orbit variables just for brevity. VS is obtained by expanding the perturbative interaction potential term in the Hamiltonian with a/aout up to the quadrupole moment and performing its double-averaging procedure. It is described by the orbital elements as VS = V0vS(e, i, ω), where V0 = Gm3a 2 16a 3 out (1 − eout) 3/2 , vS = (2 + 3e 2 )(3 cos 2 i − 1) + 15e 2 cos 2ω sin 2 i.
Introducing the following three variables: µ ≡ cos i, τ ≡ V0 na 2 t , (A12) Figure A1. The maximum and minimum values of eccentricity in terms of C KL . The red solid and blue dotted curves denote the maximum and minimum values of the eccentricity, respectively. We choose θ 2 = 0.01, 0.2, 0.4, 0.6, and 0.8. The libration type exists only for θ 2 < 0.6.
From the condition of emin ≤ emax, we have the constraints for θ and CKL: θ 2 ≤ −CKL + 1 (rotation type), θ 2 ≤ 1 5 (−2CKL + 3 − 2 √ −6CKL) (libration type), which are the same as Eqs. (4) and (5). In Fig. A1, we show some examples of emin and emax for four types of KL-oscillations. We find that the eccentricity oscillates between zero and the maximum value for the initially circular types, while it changes between two finite values (finite minimum and finite maximum values). For the libration types, there is no KL-oscillation beyond some critical value of θ 2 , while for rotation types, θ 2 reaches almost unity although the oscillation amplitude becomes smaller for larger θ 2 .
The exact half-period of the KL-oscillation TKL is defined by the time such that the eccentricity changes from the minimum value to the maximum value (Antognini 2015). It is evaluated as where τKL = ηmax η min dη dτ −1 dη.
Since τKL has order of unity, the dimensionful factor na 2 /V0 is used for rough estimation of the KL-timescale, which cor- Figure A2. Normalized KL-oscillation period τ KL in terms of θ 2 . The cyan and magenta curves denote τ KL for the libration and rotation types, respectively. responds to Eq. (7). We find where K(k) is the complete elliptic integral of the first kind with the modulus k. In Fig. A2, we show τKL.

A2 Post-Newtonian Correction
In the restricted triple system, the first order post-Newtonian (1PN) GR correction can be included by adding the correction term to the interaction potential, that is, where This correction term is derived by double-averaging the 1PN Hamiltonian of two-body relative motion (See e.g. Migaszewski & Goździewski (2011). The original Hamiltonian is obtained in Richardson & Kelly (1988)). When the corrected potential V (GR) S is used instead of VS, dimensionless potential vS is also replaced by where ǫ (1PN) is the dimensionless constant that describes the 1PN GR correction defined by Eq. (9). The basic equations for the orbital elements are the same as Eqs.(A13), (A14) and (A15) by replacing the potential vS with v (GR) S . Hence we find two conserved quantities again: C (GR) KL = CKL(η, µ, ω) + ǫ (1PN) 1 − η η .
Note that CKL(η, µ, ω) is not conserved in this case because of 1PN corrections. C

(GR) KL
coincides with the Newtonian value CKL if the orbit is circular (η = 1).

Model
Type ǫ (GR) C  Table A1. The maximum and minimum eccentricities for the PNIB models. ∆e = emax − e min gives the oscillation amplitude. τ KL and T KL are the reduced KL oscillation timescale and the real period, respectively, which are calculated based on the double-averaging method. The first row gives the Newtonian result, while the second row shows the result with post-Newtonian correction.