Second post-Newtonian Lagrangian dynamics of spinning compact binaries

The leading-order spin-orbit coupling is included in a post-Newtonian Lagrangian formulation of spinning compact binaries, which consists of the Newtonian term, first post-Newtonian (1PN) and 2PN non-spin terms and 2PN spin-spin coupling. This makes a 3PN spin-spin coupling occur in the derived Hamiltonian. The spin-spin couplings are mainly responsible for chaos in the Hamiltonians. However, the 3PN spin-spin Hamiltonian is small and has different signs, compared with the 2PN spin-spin Hamiltonian equivalent to the 2PN spin-spin Lagrangian. As a result, the probability of the occurrence of chaos in the Lagrangian formulation without the spin-orbit coupling is larger than that in the Lagrangian formulation with the spin-orbit coupling. Numerical evidences support the claim.


I. INTRODUCTION
On February 11, 2016, the LIGO Scientific Collaboration and Virgo Collaboration announced the detection of gravitational wave signals (GW150914), sent out from the inspiral and merger of a pair of black holes with masses 36M ⊙ and 29M ⊙ [1]. The gravitational wave discovery directly confirmed a major prediction of Albert Einstein's 1915 general theory of relativity. The LIGO project was originally proposed in the 1980s and its initial funding was approved in 1992. Since then, the dynamics of spinning compact binaries has received more attention. A precise theoretical waveform template is necessary to match with gravitational wave data. As a kind of description of the waveforms and dynamical evolution equations, the post-Newtonian (PN) approximation to general relativity was used. Up to now, the PN expansion of the relativistic spinning two-body problem has provided the dynamical non-spin evolution equations and the spin evolution equations to fourth post-Newtonian (4PN) order [2][3][4][5][6][7].
A key feature of the gravitational waveforms from a chaotic system is the extreme sensitivity to initial conditions [8][9][10], and therefore the chaos is a possible obstacle to the method of matched filtering. For this reason, several authors were interested in the presence or absence of chaotic behavior in the orbital dynamics of spinning black hole pairs [11][12][13][14][15][16][17]. There were three debates about this topic.
Sixteen years ago fractal methods were used to show that the 2PN harmonic-coordinates Lagrangian formulation of two spinning black holes admits chaotic behavior [11]. Here the Lagrangian includes contributions from the Newtonian, 1PN and 2PN non-spin terms and the effects of spin-orbit and spin-spin couplings. [41] However, the authors of [12] made an opposite claim on ruling out chaos in compact binary systems by finding no positive Lyapunov exponents along the fractal of [11]. The au- * Electronic address: xwu@ncu.edu.cn thors of [18,19] refuted this claim by finding positive Lyapunov exponents, and pointed out that the reason for the false Lyapunov exponents obtained in Ref. [12] lies in using the Cartesian distance between two nearby trajectories by continually rescaling the shadow trajectory. This is a debate with respect to Lyapunov exponents resulting in two different claims on the chaotic behavior of comparable mass spinning binaries. In fact, the true reason for the discrepancy was found in Ref. [20] and should be that different treatments of the Lyapunov exponents were given in the three papers [12,18,19]. The authors of [12] used the stabilizing limit values as the values of Lyapunov exponents, whereas those of [18,19] used the slopes of the fit lines. It is clear that obtaining the limit values requires more CPU times than obtaining the slopes.
Second debate is different descriptions of chaotic regions and parameter spaces. It was reported in Ref. [21] that increasing the spin magnitudes and misalignments leads to the transition to chaos, and the strength of chaos is the largest for the spins perpendicular to the orbital angular momentum. However, an entirely different description from Ref. [22] is that chaos occurs mainly when the initial spins are nearly antialigned with the orbital angular momentum for the binary configuration of masses 10M ⊙ and 10M ⊙ . These descriptions seem to be apparently conflicting, but they can all be correct, as mentioned in Ref. [23]. This is because a complicated combination of all parameters and initial conditions rather than single physical parameter or initial condition is responsible for yielding chaos. No universal rule can be given to dependence of chaos on each parameter or initial condition.
Third debate relates to different dynamical behaviors of PN Lagrangian and Hamiltonian conservative systems at the same order. In other words, the debate corresponds to a question whether the two formulations are equivalent. The equivalence of Arnowitt-Deser-Misner (ADM) Hamiltonian and harmonic coordinate Lagrangian approaches at 3PN order were proved by two groups [24,25]. This result is still correct for approximation accuracy to next-to-next-to-leading (4PN) order spin1-spin2 couplings [4]. Recently, the authors of [26] resurveyed this question. In their opinion, this equivalence means that all the known results of the ADM Hamiltonian approach can be transferred to the harmonic-coordinates Lagrangian one, or those of the harmonic-coordinates Lagrangian approach can be transferred to the ADM Hamiltonian one. In fact, the two approaches are not exactly equal but are approximately related. In this case, dynamical differences between the two approaches are not avoided. For a two-black hole system with two bodies spinning, the PN Lagrangian formulation of the Newtonian term and the 1.5PN spin-orbit coupling is non-integrable and possibly chaotic, whereas the PN Hamiltonian formulation of the Newtonian term and the 1.5PN spin-orbit coupling is integrable and nonchaotic [26]. This is due to the difference between the two formulations, 3PN spin-spin coupling. That is, the Lagrangian is exactly equivalent to the Hamiltonian plus the 3PN spin-spin term. As a result, the 3PN spin-spin effect leads to the non-integrability of the Hamiltonian equivalent to the Lagrangian. Ten years ago, similar claims were given to the two PN formulations of the two-black hole system with one bodies spinning [21,[27][28][29]. An acute debate has arisen as to whether the compact binaries of one bodies spinning exhibit chaotic behaviour. A key to this question in Ref. [30] is as follows. The construction of canonical, conjugate spin variables in Ref. [31] showed directly that four integrals of the total energy and total angular momentum in an eightdimensional phase space of a conservative Hamiltonian equivalent to the Lagrangian determine the integrability of the Lagrangian. Precisely speaking, no chaos occurs in the PN conservative Lagrangian and Hamiltonian formulations of comparable mass compact binaries with one body spinning.
One of the main results of [26,32] is that a PN Lagrangian approach at a certain order always exists a formal equivalent PN Hamiltonian. This is helpful for us to study the Lagrangian dynamics using the Hamiltonian dynamics. Following this direction, we shall revisit the 2PN ADM Lagrangian dynamics of two spinning black holes, in which the Newtonian, 1PN and 2PN non-spin terms and the 1.5PN spin-orbit and 2PN spinspin contributions are included. A comparison between the Lagrangian and related Hamiltonian dynamics will be made, and a question of how the orbit-spin coupling exerts an influence on chaos resulted from the spin-spin coupling in the Lagrangian will be particularly discussed. The present investigation is unlike the work [11], where the onset of chaos in the 2PN Lagrangian formulation was mainly shown. It is also unlike the paper [33], where the effect of the orbit-spin coupling on the strength of chaos caused by the spin-spin coupling was considered but the 1PN and 2PN non-spin terms were not included in the Lagrangian.
In our numerical computations, the velocity of light c and the constant of gravitation G are taken as geometrized units, c = G = 1.

II. POST-NEWTONIAN APPROACHES
Suppose that compact binaries have masses M 1 and M 2 , and then their total mass is M = M 1 + M 2 . Other parameters are mass ratio β = M 1 /M 2 , reduced mass µ = M 1 M 2 /M and mass parameter η = µ/M = β/(1 + β) 2 . In ADM center-of-mass coordinate system, r is a relative position of body 1 to body 2, and v is a velocity of body 1 relative to the center of mass. Let unit radial vector be n = r/r with radius r = |r|. The evolution of spinless compact binaries can be described by the following PN Lagrangian formulation The above three parts are the Newtonian term L N , first order PN contribution L 1P N and second order PN contribution L 2P N . They are expressed in Ref. [34] as In fact, these dimensionless equations deal with the use of scale transformation: r → GM r, t → GM t and L 0 → µL 0 . In terms of Legendre transformation H 0 = p · v − L 0 with momenta p = ∂L/∂v, we have the following Hamiltonian where these sub-Hamiltonians are The two Hamiltonians H 1P N and H 2P N are the results of [34]. Besides them, other higher-order PN terms can be derived from the Lagrangian L 0 . For example, a third order PN sub-Hamiltonian was given in [26] by It is obtained from coupling of the 1PN term L 1P N and the 2PN term L 2P N . Since the difference between the 2PN Lagrangian L 0 and the 2PN Hamiltonian H 0 is at least 3PN order, the two PN approaches are not exactly equivalent. Clearly, a Hamiltonian that is equivalent to the 2PN Lagrangian [42] cannot be at second order but should be at a higher enough order or an infinite order. This is one of the main results of [26]. When the two bodies spin, some spin effects should be considered. Now, the leading-order spin-spin coupling interaction L 2ss [35], as one kind of spin effect, is included in the Lagrangian L 0 . In this sense, the obtained Lagrangian becomes where Note that each spin variable is dimensionless, namely, It is only converted to a spin-spin coupling Hamiltonian Adding this term to the Hamiltonian H 0 , we have the following Hamiltonian When another kind of spin effect, the leading-order spin-orbit coupling L 1.5so [35], is further included in the above-mentioned Lagrangian L 1 , we obtain a Lagrangian formulation as follows: Unlike L 2ss , L 1.5so depends on the velocity. Therefore, the Legendre transformation results in the appearance of the leading-order spin-orbit Hamiltonian H 1.5so obtained from the coupling of the Newtonian term L N and the spin-orbit term L 1.5so , the next-order spin-orbit Hamiltonian H 2.5so obtained from the coupling of the 1PN term L 1P N and the spin-orbit term L 1.5so and the next-order spin-spin Hamiltonian H 3ss obtained from the coupling of the leading-order spin-orbit term L 1.5so and itself. These Hamiltonians are written in [26] as We take three Hamiltonians: For the three Hamiltonians, H 4 is the best approximation to the Lagrangian L 2 although the former is not exactly equivalent to the latter. As an important point to note, L 2 exhibits chaos when the two objects spin and the spin effects are L 1.5so but are not L 2ss . This is because many higher-order spin-spin couplings (such as H 3ss ), included in a higher enough order Hamiltonian equivalent to L 2 (with its equations of motion to another higher enough order), make this equivalent Hamiltonian non-integrable [26]. However, any PN conservative Lagrangian approach is always integrable and non-chaotic when only one body spins and the spin effects are not restricted to the spinorbit couplings. This is due to integrability of the equivalent Hamiltonian [30,31]. Clearly, the above-mentioned same PN order Lagrangian and Hamiltonian approaches (e.g. L 2 and H 2 ) are typically nonequivalent, and therefore they both have different dynamics to a large extent. Now, we are mainly interested in knowing how the spin-orbit coupling L 1.5so affects the chaotic behavior yielded by the spin-spin coupling L 2ss in the above Lagrangians. In other words, we should compare differences between the L 1 and L 2 dynamics. Considering that the 2PN spin-spin coupling H 2ss equivalent to L 2ss and the 3PN spin-spin coupling H 3ss associated to L 1.5so play an important role in the onset of chaos in the Hamiltonians, we should also focus on differences between the H 3 and H 4 dynamics. These discussions await the following numerical simulations. An eighth-and ninth-order Runge-Kutta-Fehlberg algorithm of variable step sizes [RKF8 (9)] is used to solve the related Lagrangian and Hamiltonian systems. This algorithm is highly precise. In fact, it gives an order of 10 −11 to the accuracy of the Lagrangian L 2 and an order of 10 −12 to the accuracy of the Hamiltonian H 2 , as shown in Fig. 1. Here, the errors of L 2 and H 2 were estimated according to two different integration paths. The error of L 2 , i.e., the error of H 2 , was obtained by applying the RKF8(9) to solve the 2PN Lagrangian equations of L 2 , while the error of H 2 was given by applying the RKF8(9) to solve the 2PN Hamiltonian equations of H 2 . Although the errors have secular changes, they are so small that the obtained numerical results should be reliable.
The evolution of the orbit in Fig. 2 demonstrates that the same order Lagrangian and Hamiltonian formulations L 2 and H 2 diverge quickly from the same starting point. This supports again the general result of [26] on the nonequivalence of the PN Lagrangian and Hamiltonian approaches at the same order.
For the given orbit, we investigate dynamical differences among some Lagrangian and Hamiltonian formulations using several methods to find chaos.

A. Chaos indicators
A power spectral analysis method is based on Fourier transformation and gives the distribution of frequencies to a certain time series. It can roughly detect chaos from order. Discrete frequencies are usually regarded as power spectra of regular orbits, whereas continuous frequencies are generally referred to as power spectra of chaotic orbits. In light of this criterion, distributions of continuous frequency spectra in Figs. 3(a) and (d)-(f) seem to show the chaoticity of the Lagrangian L 1 and the Hamiltonians H 2 , H 3 and H 4 . On the other hand, distributions of discrete frequency spectra in Figs. 3(b) and (c) describe the regularity of the Hamiltonians H 1 and the Lagrangian L 2 . It is sufficiently proved that the same PN order Lagrangian and Hamiltonian approaches (L 1 and H 1 , L 2 and H 2 ) have different dynamics. It is worth emphasizing that the method of power spectra is ambiguous to differentiate among complicated periodic orbits, quasiperiodic orbits and weakly chaotic orbits. Therefore, more reliable qualitative methods are necessarily used.
As a common method to distinguish chaos from order, a Lyapunov exponent characterizes the average exponential deviation of two nearby orbits. The varia-tional method and the two-particle one are two algorithms for calculating the Lyapunov exponent [36,37]. Although the latter method is less rigorous than the former method, its application to a complicated system is more convenient. For the use of the two-particle method, the principal Lyapunov exponent is defined as where d(0) and d(t) are separations between two neighboring orbits at times 0 and t, respectively. The best choice of the initial separation d(0) is an order of 10 −8 under the circumstance of double precision [36]. In addition, avoiding saturation of orbits needs renormalization. A global stable system [43] is chaotic if λ reaches a stabilizing positive value, whereas it is ordered when λ tends to zero. In terms of this, it can be seen clearly from Fig.  4 that the four approaches L 1 , H 2 , H 3 and H 4 with positive Lyapunov exponents are chaotic, and the two formulations H 1 and L 2 with zero Lyapunov exponents are regular. Note that each of these values of Lyapunov exponents is given after integration time 2 × 10 5 . Obtaining a reliable stabilizing limit value of Lyapunov exponent usually costs an extremely expensive computation [20].
A quicker method to find chaos is a fast Lyapunov indicator (FLI) [38,39]. It was originally calculated using the length of a tangential vector and renormalization is unnecessary. Similar to the Lyapunov exponent, this indicator can be further developed with the separation between two nearby trajectories [40]. The modified indicator is of the form An appropriate choice for renormalization within a reasonable amount of time span is important. See Ref. [40] for more details of this indicator. The FLI increases exponentially with time for a chaotic orbit, but algebraically for a regular orbit. The completely different time rates are used to distinguish between the two cases. Based on this point, the dynamical behaviors of the six approaches L 1 , L 2 , H 1 , H 2 , H 3 and H 4 can be described by the FLIs in Fig. 5. These results are the same as those given by the Lyapunov exponents in Fig. 4. Here each FLI was obtained after t = 5×10 4 . In this sense, the FLI is indeed a faster method to identify chaos than the Lyapunov exponent. Because of this advantage, the method of FLIs is widely used to sketch out the global structure of phase space or to provide some insight into dependence of chaos on single physical parameter or initial condition [23].

B. Effects of varying the mass ratio on chaos
Fixing the above initial conditions, initial spin configurations and spin parameters (χ 1 and χ 2 ), we let the mass ratio β run from 0 to 1 in increments of 0.01. For each value of β, FLI is obtained after integration time t = 5 × 10 4 . In this way, we plot Fig. 6 in which dependence of FLI for each PN approach on β is described.
It is found through a number of numerical tests that 7.5 is a threshold of FLI to distinguish between regular and chaotic orbits. A global stable orbit is chaotic if its FLI is larger than the threshold, but ordered if its FLI is smaller than the threshold. In this sense, this figure shows clearly the correspondence between the mass ratio and the orbital dynamics.
It is shown sufficiently that the PN Lagrangian and Hamiltonian approximations at the same order, (L 1 , H 1 ) and (L 2 , H 2 ), have different dynamical behaviors in significant measure. More details on the related differences are listed in Table 1. It can also be seen that the chaotic parameter space of L 1 is larger than that of L 2 . That means that the spin-orbit coupling L 1.5so makes many chaotic orbits in the system L 1 evolve into ordered orbits. [44] In other words, the probability of the occurrence of chaos in the system L 1 without L 1.5so is large, but that in the system L 2 with L 1.5so is small. Thus, L 1.5so plays an important role in weakening or suppressing the chaotic behaviors yielded by L 2ss . Here, the chaotic behaviors weakened (or suppressed) do not mean that an individual orbit must become from strongly chaotic to weakly chaotic (or from chaotic to nonchaotic) when L 1.5so is included in L 1 . This orbit may become more strongly chaotic, or may vary from order to chaos. For example, β = 0.03 in Table 1 corresponds to the regularity of L 1 but the chaoticity of L 2 . In fact, the chaotic behaviors weakened or suppressed mean decreasing the chaotic parameter space. The result obtained in the general case with the PN terms L 1P N and L 2P N is an extension to the special case without the PN terms L 1P N and L 2P N in Ref. [33]. As the authors of [33] claimed, this result is due to different signs of H 2ss (equivalent to L 2ss ) and H 3ss (associated to L 1.5so ). The two spin-spin terms are responsible for causing chaos in the Hamiltonian H 4 , and H 2ss has a more primary contribution to chaos. It is further shown in Fig. 6 and Table 1 that H 4 with the inclusion of H 3ss has weak chaos and a small chaotic parameter space, compared with H 3 with the absence of H 3ss . This is helpful for us to explain why L 1.5so can somewhat weaken or suppress the chaoticity caused by L 2ss in the PN Lagrangian system L 2 .

IV. CONCLUSIONS
When a Lagrangian function at a certain PN order is transformed into a Hamiltonian function, many additional higher-order PN terms usually occur. In this sense, the PN Lagrangian and Hamiltonian approaches at the same order are generally nonequivalent. The equivalence between the Lagrangian formulation and a Hamiltonian system often requires that the Euler-Lagrangian equations and the Hamiltonian should be up to higher enough orders or an infinite order.
For the Lagrangian formulation of spinning compact binaries, which includes the Newtonian term, 1PN and 2PN non-spin terms and 2PN spin-spin coupling, the Legendre transformation gives not only the same order PN Hamiltonian but also many additional higher-order PN terms, such as the 3PN non-spin term. Therefore, the same order Lagrangian and Hamiltonian approaches have some different dynamics. This result is confirmed by numerical simulations. This Lagrangian is non-integrable and can be chaotic under an appropriate circumstance due to the absence of a fifth integral of motion in the equivalent Hamiltonian. When the 1.5PN spin-orbit coupling is added to the Lagrangian, the 3PN spin-spin coupling appear in the derived Hamiltonian. The 3PN spinspin Hamiltonian is small and has different signs compared with the 2PN spin-spin Hamiltonian. In this sense, the probability of the occurrence of chaos in the Lagrangian formulation without the spin-orbit coupling is large, whereas that in the Lagrangian formulation with the spin-orbit coupling is small. That means that the leading-order spin-orbit coupling can somewhat weaken or suppress the chaos yielded by the leading-order spinspin coupling in the PN Lagrangian formulation. Numerical results also support this fact.