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 leads to a 3PN spin–spin coupling occurring 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 this claim.


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 e-mail: xwu@ncu.edu.cn 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-coordinate 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. 1 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 authors 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 Refs. [12,18,19]. The authors of [12] used the stabilizing limit values as the values of the 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 time than obtaining the slopes.
A second debate is concerned with 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 a 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.
A third debate is related with the different dynamical behaviors of PN Lagrangian and Hamiltonian conservative systems at the same order. Ten years ago, the 2PN harmonic coordinate Lagrangian formulation of the two-black hole system with one body spinning was numerically proved to be chaotic in [21,24], but the 2PN Arnowitt-Deser-Misner (ADM) Hamiltonian formulation of the two-black hole system with one body spinning was analytically shown to be integrable and non-chaotic in [25,26]. In fact, this debate deals with a question whether the Lagrangian and Hamiltonian formulations are equivalent. The equivalence of the ADM Hamiltonian and the harmonic coordinate Lagrangian approaches at 3PN order were proved by two groups [27][28][29]. This result is still correct for approximation accuracy to nextto-next-to-leading (4PN) order spin1-spin2 couplings [4]. In spite of this, Levin [24] found on the basis of the two opposite claims that the PN Lagrangian and Hamiltonian approaches are not exactly equal but are approximately related. It was further pointed out that the dynamical difference between the two approaches is due to replacing the acceleration in higherorder terms of the Euler-Lagrange equations with lowerorder terms but exactly deriving the equations of motion from the Hamiltonian formalism. Approximately conserving the constants of motion in the Lagrangian system but exactly conserving them in the Hamiltonian system is also a similar reason. Recently, Wu et al. [30] revisited this question. In their opinion, besides the two differences, two other differences exist. Substituting lower-order terms for the accelerations in higher-order terms of the PN harmonic coordinate generalized Lagrangian formalism depending on coordinates, velocities, and accelerations leads to the occurrence of approximations, but no such substitution is required for the ADM Hamiltonian formalism of coordinates and momenta. In addition, approximations must arise when the transformation between the ADM and harmonic coordinates is made. For the sake of a true description of the relation between the PN Lagrangian and Hamiltonian systems, the four difference sources should be avoided. An example satisfying these requirements is a special two-black hole system with two bodies spinning, whose ADM Lagrangian formulation includes the Newtonian term and the 1.5PN spin-orbit cou-pling. This Lagrangian is exactly equivalent to a full ADM Hamiltonian 2 with the Newtonian term, the 1.5PN spin-orbit coupling, and the 3PN spin-spin coupling, but it is not identical to the truncated Hamiltonian with the Newtonian term and the 1.5PN spin-orbit coupling. It is clear that the 3PN spin-spin coupling is a difference between the Lagrangian and the truncated Hamiltonian. This difference leads to the non-integrability of the Lagrangian and the integrability of the truncated Hamiltonian. Here are some details as regards this result. The truncated Hamiltonian contains five independent integrals of motion consisting of the total energy, the total angular momentum vector and the amplitude of the orbital angular momentum. As to the conserved amplitudes of the spins, they are used to construct the canonical, conjugate spin variables [31]. The five independent integrals in a ten-dimensional phase space can sufficiently support the integrability of the truncated Hamiltonian. However, the 3PN spin-spin term gives rise to the loss of the constant amplitude of the orbital angular momentum in the full Hamiltonian. Therefore, the full Hamiltonian is no longer integrable, and its equivalent Lagrangian can be chaotic under appropriate circumstances. It should be emphasized that the presence of the five independent integrals in the truncated Hamiltonian and the existence of the four independent integrals in the full Hamiltonian (i.e. the Lagrangian) are independent of the choice of coordinates or gauges. Even if the two Hamiltonians are expressed in terms of action-angle variables rather than coordinates and are coordinate (or gauge) invariant [32][33][34], these integrals still exist. The integrability or non-integrability of these systems does not depend on the choice of coordinates. As an important result given in [30], truncating higher-order PN terms are responsible for causing dynamical differences between the PN Lagrangian and Hamiltonian formulations at the same order. In particular, these differences cannot always be avoided regardless of the Hamiltonians expressed in terms of the action-angle variables if the truncation of higher-order terms occurs. The equivalence claimed in Refs. [4,25,26] is unlike that of [30] and means that all the known results of the ADM Hamiltonian approach can be transferred to those of the harmoniccoordinate Lagrangian approach, or all the known results of the harmonic-coordinate Lagrangian approach can be transferred to those of the ADM Hamiltonian approach. Noticing that the Euler-Lagrange equations and the constants of motion are generally approximate in a PN Lagrangian system but exact in its related Hamiltonian, perhaps one would believe in the preference of the Hamiltonian formalism over the Lagrangian formalism in the study of chaos. This point of view is problematic because it is unclear which of the Lagrangian and Hamiltonian formalisms can truly describe a real physical system. For example, is the Hamiltonian formalism truer in describing the two-black hole system including the Newtonian term and the 1.5PN spin-orbit coupling than the Lagrangian formalism? No one has answered this question. In fact, a Lagrangian formalism is approximate in general at the same PN order if it is derived from a Hamiltonian formalism; on the contrary, a Hamiltonian formalism is approximate in most cases if it is derived from a Lagrangian formalism. When a Lagrangian formalism is exactly equivalent to a Hamiltonian formalism, it is of course better to apply the Hamiltonian formalism to investigate the chaotic dynamics. This is because the Hamiltonian formalism can exactly provide the equations of motion and some constants of motion, and it has many advantages on the properties of a canonical system. 3 Fortunately, another important result of [30] is that for a lower-order PN Lagrangian formulation with Euler-Lagrange equations to an infinite PN order there always exists a formally equivalent PN Hamiltonian at the infinite order from a theoretical point of view or a certain finite order from a numerical point of view. This result is supported analytically and numerically via a special 1PN Lagrangian formulation of a relativistic circular restricted three-body problem with the Euler-Lagrange equations and the equivalent Hamiltonian as a converging Taylor series [36]. Thus, the integrability or non-integrability of a PN Lagrangian formalism can be shown in terms of that of its formal equivalent PN Hamiltonian. Since a conservative PN Hamiltonian of comparable mass compact binaries with one body spinning has four integrals of the total energy and total angular momentum in an eight-dimensional phase space, it is integrable and nonchaotic. This also shows the integrability of any conservative PN Lagrangian of comparable mass compact binaries with one body spinning. Stated succinctly, neither the conservative PN Lagrangian dynamics of comparable mass compact binaries with one body spinning nor the conservative PN Hamiltonian dynamics of comparable mass compact binaries with one body spinning is chaotic [31]. The debate on the two different claims of the chaotic behavior between the PN Lagrangian and Hamiltonian formulations of comparable mass compact binaries with one body spinning was ended.
As was mentioned above, one of the main results of [30,36] is that for a PN Lagrangian approach at a certain order there always exists a formally equivalent PN Hamiltonian. This is helpful for us in our study of the Lagrangian dynamics using the Hamiltonian dynamics. Following this direction, we shall revisit the 2PN ADM Lagrangian dynamics of two spinning black holes [11], in which the Newtonian, 1PN and 2PN non-spin terms and the 1.5PN spinorbit and 2PN spin-spin contributions are included. A comparison between the Lagrangian and related Hamiltonian dynamics will be made, and the question of how the orbitspin coupling exerts influence on chaos resulting from the spin-spin coupling in the Lagrangian will be particularly discussed. The present investigation is unlike [11], where the onset of chaos in the 2PN Lagrangian formulation was mainly shown. It is also unlike Ref. [38], 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.

Post-Newtonian approaches
Suppose that compact binaries have masses M 1 and M 2 ; then their total mass is M = M 1 + M 2 . Other parameters are the mass ratio β = M 1 /M 2 , the reduced mass μ = M 1 M 2 /M and the mass parameter η = μ/M = β/(1+β) 2 . In the ADM center-of-mass coordinate system, r is a relative position of body 1 to body 2, and v is a relative velocity. Let the 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 , the firstorder PN contribution L 1PN , and the second-order PN contribution L 2PN . They are expressed in [39] as In fact, these dimensionless equations deal with the use of a scale transformation: r → G Mr, t → G Mt, and L 0 → μL 0 . In terms of a 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 1PN and H 2PN are the results of Refs.
[32] and [39]. Besides them, other higher-order PN terms can be derived from the Lagrangian L 0 . For example, a thirdorder PN sub-Hamiltonian was given in [30] by It is due to the coupling of the 1PN term L 1PN and the 2PN term L 2PN in the Lagrangian L 0 . In fact, this Hamiltonian plus the 3PN Hamiltonian from the 3PN Lagrangian of [39] is in perfect agreement with the full 3PN Hamiltonian derived in [39][40][41][42][43]. 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 4 cannot be 4 Here, the equations of motion derived from this Lagrangian are required to reach a high enough order or an infinite order.
at second order but should be at a high enough order or an infinite order. This is one of the main results of [30]. When the two bodies spin, some spin effects should be considered. Now, the leading-order spin-spin coupling interaction L 2ss , as one kind of spin effect, is included in the Lagrangian L 0 . In this sense, the Lagrangian obtained becomes where Note that each spin variable is dimensionless, namely, S i = S iŜi with spin magnitude S i = χ i M i 2 /M 2 (0 ≤ χ i ≤ 1) and three-dimensional unit vectorŜ i . Because L 2ss is independent of the velocity, it does not couple the 1PN term L 1PN or the 2PN term L 2PN via the Legendre transformation. It is only converted to a spin-spin coupling Hamiltonian, The leading-order spin-spin Hamiltonian H 2ss is the same as that in Ref. [44], and it includes the leading-order spin(a)spin(b) part and the leading-order spin(a)-spin(a) part. It is also consistent with that in Refs. [45,46], but the scale factors are different. Adding this term to the Hamiltonian H 0 , we have the following Hamiltonian: When another kind of spin effect, the leading-order spinorbit coupling L 1.5so , is further included in the abovementioned Lagrangian L 1 , we obtain a Lagrangian formulation [11] 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 1PN 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 [30] as Notice that H 1.5so is the result of [44]. Additionally, the nextto-leading (3PN) order spin-spin coupling Hamiltonian H 3ss is only a part of that in Refs. [45,46], but different scale factors are used. If the full 3PN spin-spin Hamiltonian in Refs. [45,46] is derived from a PN Lagrangian formulation, it is from the 3PN order spin-spin couplings L 3ss in the Lagrangian as well as the coupling of the leading-order spinorbit term L 1.5so and itself in the Lagrangian. That is to say, this Lagrangian should be equal to L 2 plus L 3ss . Similarly, H 2.5so is not a full next-to-leading-order spin-orbit coupling Hamiltonian. In fact, the full 2.5PN spin-orbit Hamiltonian should be equal to the Hamiltonian (17) plus the 2.5PN spinorbit Hamiltonian from the 2.5PN spin-orbit Lagrangian. Now, we take three Hamiltonians: Of the three Hamiltonians, H 4 is the best approximation to the Lagrangian L 2 although the two are not exactly equivalent. Let us investigate the integrability or non-integrability of these PN Lagrangian and Hamiltonian systems. First, we consider the case without the terms L 1PN , L 2PN and L 2ss in the Lagrangian L 2 , i.e.L A = L N + L 1.5so . As was claimed in [30],L A is not equivalent toH A1 = H N + H 1.5so but is exactly equivalent toH A2 = H N + H 1.5so + H 3ss . The systemH A1 holds five exact constants of motion: these are the total energy E =H A1 , the total angular momentum vector J = L + S 1 + S 2 , and the amplitude |L| of the orbital angular momentum L = r × p. However, the amplitude |L| is no longer invariant due to the presence of H 3ss in the systemH A2 . When the conserved amplitudes of the two spins are used to construct the canonical, conjugate spin variables [31], each ofH A1 andH A2 has a ten-dimensional phase space. Consequently, the dynamics ofH A1 is integrable and regular but that ofH A2 (i.e. L A ) is non-integrable and possibly chaotic. Numerical evidence showed the chaoticity of the systemL A and the regularity of the systemH A1 [30]. 5 Second, when the terms L 1PN and L 2PN are included in the systemL A , i.e.L B = L A + L 1PN + L 2PN , no exactly equivalent Hamiltonians but approximately related Hamiltonians can be given, e.g., For the case of two bodies spinning, the five exact integrals (E, J and |L|) still show the integrability ofH B1 andH B2 , but the four exact integrals (E and J) determine the nonintegrability ofH B3 . As far as the dynamics ofL B is concerned, the Euler-Lagrange equations are accurate to second order.L B should have four exact integrals of the total energy and the total angular momentum vector from the physical point of view, but give only four approximate integrals from the mathematical point of view. The approximate energy integral may be one ofH B1 ,H B2 andH B3 , and the approximate angular momentum integral is J with p = v. Regardless of these integrals being approximate or exact, the inclusion of L 1PN and L 2PN does not destroy the non-integrability of L A . That means thatL B is not integrable. Third, when L 2ss is added toL B , L 2 is obtained and is similar to the second case. In this case, the presence of the 2PN spin-spin coupling or the 3PN spin-spin coupling in the systems L 1 , L 2 , H 1 , H 2 , H 3 , and H 4 leads to the non-integrability of these systems without doubt. As was mentioned in the Introduction, the existence of the exact (or approximate) integrals of motion in these systems is physical and does not depend on the choice of coordinates or gauges. That is to say, the integrability or non-integrability of these systems is completely independent of the choice. If k first integrals exist in an autonomous canonical Hamiltonian system, then k cyclic coordinates appear in principle in this Hamiltonian expressed by action-angle variables that render the Hamiltonian coordinate (or gauge) invariant. It is easy to apply the action-angle variables for the circular orbits of some integrable Hamiltonian systems without the spin-spin couplings in [32][33][34], but it is difficult to do so for generic orbits or the inclusion of the spin-spin couplings in the present case.
Clearly, the above-mentioned same PN order Lagrangian and Hamiltonian approaches (e.g. L 2 and H 2 ) have to a large extent dynamical differences because of many higherorder terms having been truncated. Only if these higherorder terms truncated would exist, would the differences between the PN Lagrangian and Hamiltonian approaches be present, and do approximations also occur even if actionangle variables are used. Now, we are mainly interested in 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 and see the differences between the L 1 and L 2 dynamics. Considering that the 2PN spin-spin coupling H 2ss is equivalent to L 2ss and the 3PN spin-spin coupling H 3ss associated to L 1.5so plays 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 will properly be presented in the context of 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 [30] on the non-equivalence 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.

Chaos indicators
A power spectral analysis method is based on a 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, the distributions of continuous frequency spectra in Fig. 3a 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, the distributions of discrete frequency spectra in 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 in differentiating among complicated periodic orbits, quasiperiodic orbits, and weakly chaotic orbits. Therefore, more reliable qualitative methods should necessarily be used. As a common method to distinguish chaos from order, a Lyapunov exponent characterizes the average exponential deviation of two nearby orbits. The variational method and the two-particle one are two algorithms for calculating the Lyapunov exponent [47,48]. 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 twoparticle method, the principal Lyapunov exponent is defined as where d(0) and d(t) are the 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 [47]. In addition, avoiding saturation of orbits needs renormalization. A global stable system 6 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 punov exponents are regular. Note that each of these values of Lyapunov exponents is given after an integration time of 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) [49,50]. 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 [51]. The modified indicator is of the form An appropriate choice for renormalization within a reasonable amount of time span is important. See Ref. [51] 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 . As a notable question, invariant values of the Lyapunov exponents and the FLIs should require the use of proper time and distances [48]. However, the above computations of the Lyapunov exponents and the FLIs use coordinate distances and time that make the Lyapunov exponents and the FLIs not invariant in general relativity. In spite of this, the dynamical features given by the Lyapunov exponents and the FLIs with coordinate distances and time should be the same as those given by the Lyapunov exponents and the FLIs with proper distances and time. This is because the difference between the distance and the invariant distance is negligible in the present problems and the ratio of the proper time to the coordinate time is a positive finite value [52,53].
Because the FLI is indeed a faster method to identify chaos than the Lyapunov exponent, it is widely used to sketch the global structure of phase space or to provide some insight into the dependence of chaos on a single physical parameter or initial condition [23].

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 an integration time t = 5 × 10 4 . In this way, we plot Fig. 6 in which the dependence of the 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 globally 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 a significant measure. More details of 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. 7 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) are given from a global point of view, but do not mean that an individual orbit must become from strongly chaotic to weakly chaotic (or from chaotic to non-chaotic) when L 1.5so is included in L 1 . This orbit may become more strongly chaotic, or it 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 1PN and L 2PN is an extension to the special case without the PN terms L 1PN and L 2PN in Ref. [38]. As the authors of [38] 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 in 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 .

Conclusions
When a Lagrangian formalism at a certain PN order is transformed into a Hamiltonian formalism, many additional higher-order PN terms usually occur. In this sense, the PN Lagrangian and Hamiltonian approaches at the same order are generally nonequivalent, and these higher-order PN terms truncated are responsible for causing dynamical differences between the PN Lagrangian and Hamiltonian formulations at the same order. When the physical gauge invariant observ- ables are used instead of the coordinates, approximations (or differences) still occur if the truncation of higher-order PN terms exists. The equivalence between the Lagrangian formulation and a Hamiltonian system often requires that the Euler-Lagrangian equations and the Hamiltonian should be up to high enough orders or an infinite order. The integrability or non-integrability of the Lagrangian formulation is that of its equivalent Hamiltonian. For the 2PN 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 through numerical simu-lations. 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 appears in the derived Hamiltonian. The 3PN spin-spin 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 leadingorder spin-orbit coupling can somewhat weaken or suppress the chaos yielded by the leading-order spin-spin coupling in the PN Lagrangian formulation. Numerical results also support this fact. It should be emphasized that the integrability or non-integrability of these systems depends on the number of the first integrals of these systems and is physical, but it does not depend on the choice of coordinates. In other words, the property given by the coordinates and the one given by the physical gauge invariant observables are the same. Although the chaos indicators use the coordinate distances and time, they should give these systems the same dynamical information described by the chaos indicators with the proper distances and time. Here are two reasons. The difference between the distance and the invariant distance is negligible in the present problems, and the ratio of the proper time to the coordinate time is a positive finite value.
The PN Lagrangian or Hamiltonian formulations have been actually used in the construction of theoretical gravitational wave templates and the successful detection of gravitational wave signals is based on the absence of chaos in the considered dynamical systems. Of course, the future detection of gravitational wave signals also requires that chaos be avoided when the PN formulations are planned to be used in the construction of theoretical gravitational wave templates.