Non-truncated strategy to exactly integrate the post-Newtonian Lagrangian circular restricted three-body problem

In this study, we present a novel non-truncated strategy by accompanying the fixed-point iteration with traditional numerical integrators. The proposed non-truncated strategy aims to exactly integrate implicit motion equations that are directly derived from the Lagrangian of the post-Newtonian circular restricted three-body problem. In comparison with the commonly used truncated approach, which cannot exactly but approximately preserve the generalized Jacobian constant (or energy) of the original Lagrangian system, the proposed non-truncated strategy has been determined to preserve this constant well. In fact, the non-truncated strategy and the truncated approach have a difference at second post-Newtonian order. Based on Kolmogorov–Arnold–Moser theory, this difference from the truncation in the equations of motion may lead to destroying the orbital configuration, dynamical behavior of order and chaos, and conservation of the post-Newtonian circular restricted three-body problem. The non-truncated strategy proposed in this study can avoid all these drawbacks and provide highly reliable and accurate numerical solutions for the post-Newtonian Lagrangian dynamics. Finally, numerical results show that the non-truncated strategy can preserve the generalized Jacobian constant in the accuracy of O(10-12)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathscr {O}}(10^{-12})$$\end{document}, whereas the truncated approach at the first post-Newtonian (1PN) order only has an accuracy of O(10-3)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathscr {O}}(10^{-3})$$\end{document}. Moreover, several orbits are observed to be escaping from the bounded region in the 1PN truncated system via the truncated strategy, but these escaping orbits are unobserved via the non-truncated strategy.


Introduction
By considering the relativistic effect, the post-Newtonian (PN) circular restricted three-body problem describes the a e-mail: bxhanm@126.com dynamics of a third body, which moves in a plane of two primaries that drift in a circular orbit around their center of mass. Relativistic correction provides a more accurate description than the Newtonian approach for the bodies in the solar system [1][2][3][4]. Given the so-called chaotic amplification effect on the relativistic problem [2,5], an effective approach to obtaining highly accurate numerical solutions is required once numerical integrators are applied during the study on the dynamics of this type of problems.
Recently, considerable research on the dynamics of the PN circular restricted three-body problem has been conducted. In [5], Huang and Wu studied the relativistic effect of the separation between two primaries on the dynamics by applying scaling transformations to distance and time; these authors used the equations derived by Maindl and Dvorak [1]. To numerically preserve the generalized Jacobian constant, the authors of [6] presented an alternative system for the PN circular restricted three-body problem using the Einstein-Infeld-Hoffmann (EIH) formalism up to the 1PN order. Furthermore, the authors in [6] claimed that c = 10 4 is suitable for the case of a = 1 in the solar system. Recently, the authors of [7] conduct further research on the orbital dynamics, such as equilibrium points and Hill region configurations, for the PN planar circular restricted Sun-Jupiter system by considering the motion equations derived in [6].
The derived motion equations in [5] and [6] are both up to 1PN order, but they do not coincide with each other completely because their derivations are conducted by different approaches. For example, Huang and Wu selected c = 1 for an extreme case for their numerical simulations in [5], whereas Dubeibe et al. [6] selected c = 10 4 for the PN circular restricted three-body problem in the solar system and showed that the generalized Jacobian integral can be preserved with high accuracy. The different selection of c leads to the different behavior of constant conservation and dynamics for several orbits in [5] and [6]. Despite the abovementioned difference, the two system derived in [5] and [6] will be reduced to the same classical Newtonian case once c → ∞. In addition, the results in [5] and [6] support the occurrence of chaos in the PN circular restricted three-body problem.
Typically, the motion equations derived from the PN Lagrangian via the Euler-Lagrange formula are implicit motion equations in terms of accelerations. To conduct quantitative studies with numerical integrators, the explicit expressions of accelerations are required. During this procedure, truncation is a commonly used approach, which is conducted by substituting low-order equations of motion to high-order accelerations or by neglecting higher-order terms. This strategy is extensively used to investigate the dynamics of the PN Lagrangian compact binaries [8][9][10][11][12][13], the equivalence analysis between the PN Lagrangian and the PN Hamiltonian approaches [14][15][16][17][18][19][20], and the dynamics of the PN circular restricted three-body problem [1,5,6]. Among the discussion of the equivalence between the PN Lagrangian and the PN Hamiltonian approaches, the authors of [14][15][16]19] showed the physical equivalence of the two approaches at the same PN order. However, literature [17,21] claimed that the Lagrangian and the Hamiltonian approaches by including spin effects may have different properties of integrability and non-integrability or different dynamical behaviors of order and chaos.
It should be noticed that the truncation for the implicit motion equations may lead to a vast gap between the dynamics of the approximated motion equations and that of the exact motion equations directly derived from the PN Lagrangian because the truncated motion equations no longer correspond to the original PN Lagrangian. This phenomenon is inevitable for highly nonlinear systems where chaos occurs due to the following two facts. First, chaotic trajectory is exponentially sensitive to initial conditions. Second, in the sense of backward error analysis [22], a trajectory in the system of truncated motion equations with exact initial values can be interpreted as a perturbed trajectory in the system of exact motion equations with perturbed initial values. Thus, a small truncation error, even in the magnitude of O(10 −16 ) (see [6]) for the exact motion equations, can still yield the global error in a magnitude of a large order after the exponential amplification for this small truncation error. This phenomenon suggests that we must integrate the exact implicit motion equations rather than the truncated motion equations once the investigation on the dynamics of the PN Lagrangian systems via numerical simulations is conducted.
On account of the abovementioned arguments, in this study, we aim to establish an effective approach to integrating the exact implicit motion equations that are directly derived from the PN Lagrangian via the Euler-Lagrange formula, rather than the truncated motion equations. To this end, we propose a non-truncated strategy that uses the traditional numerical integrators accompanying with the fixed-point iteration to integrate the exact implicit motion equations, and these equations exactly correspond to the PN Lagrangian of the problem. Then, with this new approach, we further study on the dynamics of the PN circular restricted three-body problem. Numerical solutions obtained by this new approach can preserve the Jacobian constant well for all possible values of c and a. Furthermore, in comparison with the truncated approach, we reveal that the truncation for the motion equations will lead to abrupt changes in the orbital configuration and dynamical behavior. The results obtained from the truncated system essentially deviate from those of the exact system that corresponds to the original PN Lagrangian.
The remainder of this paper is organized as follows. In Sect. 2, we first present the 1PN Lagrangian formulation of the PN circular restricted three-body problem and then derive the exact implicit motion equations and the 1PN truncated motion equations. In Sect. 3, we introduce and analyze the new non-truncated strategy that can exactly integrate the implicit motion equations. In addition, the practical implementation issue of the non-truncated strategy is discussed. Section 4 is mainly concerned with the comparison between the new non-truncated strategy and the commonly used truncated approach by conducting numerical simulations on the PN circular restricted three-body problem. Conclusions are provided in the last section.

Formulation of post-Newtonian Lagrangian
In a planar circular restricted three-body problem, the two primaries m 1 and m 2 , separated by distance a, move in a circular orbit around their center of mass with an angular speed ω relative to an inertial frame, and they are at positions r i (i = 1, 2) with velocities v i (i = 1, 2). Then, according to the PN gravitational theory of Einstein, Infeld, and Hoffmann [23], the 1PN Lagrangian of the third particle with mass m, position r = (X, Y ), and velocity u = (Ẋ ,Ẏ ) can be expressed by (see [24]) where n i = (r − r i )/|r − r i | is the unit vector, c is the velocity of light, and G denotes the gravitational constant.
By introducing the transformations X = x cos ωt − y sin ωt, we can transfer the Lagrangian (1) of the third body in the inertial frame to the rotating frame as follows where L = (L +mc 2 )/m (see [1,5]). In this formulation, L 0 denotes the classical Newtonian Lagrangian of the circular restricted three-body problem, L 1 represents a 1PN contribution to the circular orbital frequency, and L 2 is another 1PN contribution from the primaries to the third body. The detailed expressions of L 0 , L 1 , and L 2 can be found in [1,5] and expressed as where Given that the Lagrangian (3) is approximated at the 1PN order, ω 2 calculated by Contopoulos [25] is also at the 1PN order and can be expressed as ω 2 = ω 2 0 (1+2ω 1 /c 2 ).

Exact and truncated motion equations
Because the dissipative terms are excluded from the PN Lagrangian (3), this system is confirmed to be conservative. With the generalized momentum transformation the Lagrangian (3) corresponds to a conserved energy Besides, the Lagrangian (3) also corresponds to a Hamiltonian Although E(ẋ,ẏ, x, y) and H ( p x , p y , x, y) seem to be the same, they have different expressions, because E is a function of the velocities and coordinates, but H is another function of the momenta and coordinates. In addition, E is usually an exact energy of the Lagrangian due to the absence of truncation. However, H will be approximately given in general when the velocities are expressed in terms of the momenta from nonlinear functions of the velocities and momenta, and the truncation has to be performed [17]. Thus, H may be an approximate energy of the Lagrangian. That means E and H are not exactly equal. In an analogy with the Newtonian case, the generalized Jacobian constant in the PN case can be defined as whose detailed expression is indicated by (A.1). By virtue of the Euler-Lagrangian formula d dt and according to the post-Newtonian Lagrangian (3), we can derive the implicit motion equations in terms of the accelerationsẍ = d 2 x/dt 2 andÿ = d 2 y/dt 2 as follows: where the functions F 1 , F 2 , G 1 , and G 2 are detailed in Appendix A. Considering the linearity ofẍ andÿ in (A.2) and (A.3), we can expressẍ andÿ by x, y,ẋ, andẏ theoretically as by solving the linear Eq. (12). However, a careful calculation on the computer shows that the detailed expressions of W 1 and W 2 are unmanageable to explicitly express in the present paper (In particular, the two expressions are over thousands of lines), thereby leading to the obstacle in the theoretical analysis and numerical integration for the PN circular restricted three-body problem. A similar situation also occurs in the PN Lagrangian dynamics of compact binaries, in which the detailed expressions of accelerations cannot be obtained due to the high nonlinearity of the PN Lagrangian formulation.
On noting that the Eq. (12) are much simpler than (13), the functions F 1 , F 2 , G 1 , and G 2 can be easily calculated in a computer for certain values of x, y,ẋ,ẏ,ẍ, andÿ. A common strategy to studying the dynamics of the PN Lagrangian system is using the truncated motion equations, which are only approximate to the exact implicit motion equations at a certain PN order. For example, W 1 and W 2 can be expressed in the series of c −1 as where This strategy is extensively used in the PN Lagrangian of compact binaries [8][9][10][11][12][13], the equivalence analysis between the PN Hamiltonian and the PN Lagrangian approaches at the same order [14][15][16][17][18][19][20], and the dynamics of the PN circular restricted three-body problem [1,5,6]. Following this idea, the approximately truncated motion equations at the 1PN order for (13) are derived in [5] and expressed as and where P and Q are written as

Theoretical difference of dynamics
Given that (16)- (17) are only approximations to (13) at the 1PN level, the energy (8) or the generalized Jacobian constant (10) is also an approximated integral for the truncated system (16)- (17), and this approximation is noted at the 1PN level. Despite the convenient explicit formulations of (16)-(17), their truncations for the exact motion Eq. (13) of the PN Lagrangian may lead to certain serious misleading knowledge on the dynamics. Suppose that X(t) = x(t), y(t),ẋ(t),ẏ(t) is an exact trajectory in the phase space, which obeys a system of the exact motion Eq. (13) with initial values X 0 = (x 0 , y 0 ,ẋ 0 ,ẏ 0 ) at the initial time t 0 . This trajectory X(t) corresponds to an initial energy E 0 = E(ẋ 0 ,ẏ 0 , x 0 , y 0 ) in accordance with (8). Then, we denote the trajectory obtained from integrating the system of the k-PN truncated motion equations with the same initial values X 0 by X(t) = ( x(t), y(t),˙ x(t),˙ y(t)). Since X(t) and X(t) evolve in different systems, a most remarkable point is the non-conservation of the energy (or generalized Jacobian constant) for X(t), i.e., the numerical energy , y(t)) of the trajectory X(t) cannot hold the constant E 0 any longer. Furthermore, from the theory of backward error analysis, we yield that the error of E(t) − E 0 will be at least in the magnitude of O(c −2(k+1) ).
Another drawback of the truncation approach lies in its possible destruction of the orbital configuration or even the dynamical behavior of order and chaos for certain trajectories when the system under consideration is a highly complex problem where chaos may occur. To exhibit this point clearly, we regard the trajectory X(t) evolving in the k-PN truncated system as a perturbed trajectory of the exact one X(t). This perturbation is considered as a perturbation to the initial values X 0 also at the k-PN level. That is, we can consider X(t) to evolve in the exact non-truncated system but with the perturbed initial values X 0 + O(c −2(k+1) ). Then, the difference between X(t) and X(t) could be analyzed in detail via perturbation theory. We present three cases for the difference below. First, if X(t) is located in the region of regular orbits, and O(c −2(k+1) ) is sufficiently small, then it follows from Kolmogorov-Arnold-Moser (KAM) theory (Sec. X of [22]) that X(t) is also regular, and the difference between X(t) and X(t) is roughly in the magnitude of O(c −2(k+1) ). Second, a subtle and complicated case is that X(t) remains in the regular region, while the initial perturbation X 0 + O(c −2(k+1) ) is beyond the regular region that makes X(t) a chaotic orbit. Then, X(t) possesses a distinct dynamical behavior over X(t). The final case is that both X(t) and X(t) are in a chaotic region, in which their difference will be hardly determined and possible in a magnitude of order that is much larger than O(c −2(k+1) ) given the exponential sensitivity of a chaotic orbit to initial conditions. As previously discussed, for the first case in which the exact orbit X(t) and its perturbed counterpart X(t) are regular, the difference error in the magnitude of O(c −2(k+1) ) can become insignificant via two approaches. One approach is fixing c in a large magnitude that O(c −2(k+1) ) = O(10 −16 ) holds for a fixed k. Another approach is enlarging k under a certain c to also satisfy O(c −2(k+1) ) = O(10 −16 ). Thus, the truncated and exact systems are thought of as equivalent with each other in the sense of numerical computation. However, the two approaches can increase the computation cost or significantly confine their potential applications for the fixed c and k.
Furthermore, it will disappoint us in the last two cases in which at least one orbit, i.e., X(t) or X(t), is chaotic. Without loss of generality, we suppose that the exact orbit X(t) is chaotic with a positive Lyapunov characteristic exponent λ [26][27][28], and | X 0 − X 0 | = O(c −2(k+1) ) = C × 10 −n with C, n > 0. Given that the chaotic orbit X(t) is exponentially sensitive to the initial conditions X 0 , a rough estima-tion yields | X(t) − X(t)| = Ce λ(t−t 0 ) × 10 −n . Thus, even an initial perturbation in the magnitude of O(10 −16 ) can be amplified to a vast perturbation to X(t) provided that t − t 0 is sufficiently large. That is, the perturbed orbit X(t) may significantly differ from the exact orbit X(t) in the orbital configuration and dynamical behavior. However, no effective approach to decreasing the difference between X(t) and X(t) is available.
An important point is that the system of motion Eq. (13) exactly corresponds to the Lagrangian (3). Thus, it is a conservative system in which the orbits are all bounded in a limited region. However, the truncated system (16)-(17) only approximately corresponds to (3) at the 1PN level, and its conservation as a dynamical system is indefinite, because we are unsure if such a first integral corresponding to the system (16)- (17) exists. Thus, a bounded orbit in (13) may escape from the limited region when it is evolving in the truncated system (16)-(17) considering the dissipation of energy.
Finally, in the case of non-relativistic limit c → ∞ (thus, 1/c 2 → 0), the exact motion Eq. (13) and the truncated motion Eqs. (16)-(17) will be reduced to the Newtonian case with a Lagrangian L 0 , which corresponds to the Jacobian constant [24]

A class of new implicit numerical integration schemes
In this section, we focus on the exact integration of the fully exact implicit motion Eq. (12) for the PN circular restricted three-body problem. For simplicity, we commence with the n-dimensional autonomous system of the implicit first-order ordinary differential equations where the initial conditions y(t 0 ) = y 0 ∈ R n . Although the implicit motion Eq. (12) under consideration are secondorder, we can transform these motion equations into the form of (21) by adding variables v x =ẋ and v y =ẏ.
In this study, we take the explicit s-stage Runge-Kutta of order r as an example to demonstrate the exact integration of the implicit system (21). If (21) can be equivalently written in an abstract form aṡ then applying the s-stage Runge-Kutta method to (22) yields the following numerical scheme: where a i j , b i , and c i are coefficients of the Runge-Kutta method, y n ≈ y(t n ) and y n+1 ≈ y(t n+1 ) are the numerical solutions, h is the time stepsize, K i ≈ẏ(t n +c i h) are internal stages, and t n = t 0 + nh.
Considering that the explicit form (22) cannot be obtained once the function F( y,ẏ) in (21) demonstrates a high nonlinearity with respect toẏ, we can only obtain an equivalent implicit forṁ rather than the explicit form (22). In the numerical scheme (23), by replacing the function F( y) with F ( y,ẏ) we obtain a class of new numerical schemes:

Analysis of the new scheme
The new numerical scheme (25) is specially designed for the implicit differential equation (21), which is the aim of this paper. The main novelty of the new scheme (25) is that it can solve the exact implicit motion equations of PN Lagrangian systems without any truncation, whereas other existing integrators, such as Runge-Kutta methods (23), symplectic integrators [22,[29][30][31], manifold correction methods [32][33][34][35] and extended space methods [36,37], require the truncation for the exact implicit motion equations to obtain a system of explicit differential equations when they are applied to numerically solve such a problem. As we have analyzed in Sect. 2.2, this truncation for the exact implicit motion equations may result in an essential changing for the dynamics of the system. That is, the new scheme can solve the PN Lagrangian systems more accurately than the classical integrators that are designed for general explicit motion equations.
Given that implicit formulations (21) and (24) are equivalent to the explicit formulation (22), and the new scheme (25) is obtained by replacing the function (22) with (24) in the Runge-Kutta method (23), we can yield that if (23) is of order r for (22), then the new scheme (25) is also of order r for (24). Therefore, we can conveniently borrow the coefficients of Runge-Kutta methods to construct high-order schemes in the new form of (25). An example is the classical fourth-order Runge-Kutta method, whose coefficients can be expressed by the following Butcher tableau: . (26) When these coefficients are used in the new scheme (25), the integrator is certainly fourth-order for (21).
Although the new scheme (25) is designed for the implicit differential equation (21), it can still be applied to explicit differential equations. A thorough argument indicates that the new scheme (25) will be reduced to the Runge-Kutta method (23), provided that this method is used to solve a system of explicit differential equations. In this sense, it can be said that the new scheme (25) is the extension of the Runge-Kutta method (23) from explicit differential equations to general implicit differential equations.
One point must be emphasized that, although (25) is adapted from the explicit Runge-Kutta method (23), this method is essentially an implicit integrator, and the implicit characteristic of (25) is clearly due to (21). Considering this point, we can conclude that only implicit numerical integrators rather than explicit ones, exist for the implicit differential equation (21). Therefore, the implicit numerical integrators, such as the diagonally or fully implicit Runge-Kutta are preferred for (21), because these integrators possess certain superiorities over the explicit integrators in numerical properties (e.g., the numerical stability and global truncation error).
Due to the equivalence between the PN Lagrangian (3) and the PN Hamiltonian (9), structure-preserving (symplectic, symmetric, or energy-preserving) properties [29][30][31][36][37][38][39] of the new scheme (25) can be remarkable. However, since the structure preservation is mainly for Hamiltonian systems, whether a new scheme (25) adapted from a structurepreserving Runge-Kutta method (23) still possesses such structure-preserving property remains an important but difficult problem. This difficulty mainly occurs during the transformation (7). For example, if ( Q n+1 ,Q n+1 ) are numerical solutions solved by the scheme (25) after n + 1 steps, then we must check whether the equation on the second exterior differential forms ∧ holds or not to verify the symplecticity of the scheme (25) by noting However, because of the high complexity of Lagrangian L , this procedure is difficult to accomplish in a general approach. Thus, the study on structure-preserving (symplectic, symmetric, or energy-preserving) properties of the new scheme (25) for implicit differential equations is not the subject of this study, but it will be our further consideration in the future work. Finally, we mention again that integrators, such as the extended space method [36,37], that are specially designed for Hamiltonian systems cannot be directly applied to the implicit differential equation (21).

Implementation issue
In this subsection, we focus on the implementation issue of (25). Considering the implicity of (25), iterations are required to obtain K i . As recommended by Hairer et al.
(see Sec. VIII.6 in [22]), the fixed-point iteration is preferred for the solutions of the implicit numerical scheme (25). To increase the convergence of the fixed-point iteration, we must select the iteration function F ( y,ẏ) appropriately. A possible choice for F ( y,ẏ) is presented below: where 0 < γ 1. From the theory of iteration convergence, the iteration function F ( y,ẏ) in the form of (28) will significantly increase the convergence of the fixed-point iterations for (25). Now, we turn to the PN circular restricted three-body problem. From the detailed expressions of (A.2) and (A.3) in Appendix A, we yielḋ which is a system of implicit first-order differential equations on variables (x, y, v x , v y ). Notably, (29) obeys the form of (28) once c 1. Thus, we can expect a fast convergence for the fixed-point iteration provided that the numerical scheme (25) is applied to the system (29).
We emphasize that the new numerical scheme (25) will be implemented without any truncation for the fully exact implicit motion Eq. (12) and hence is a non-truncated strategy to exactly integrate the PN Lagrangian (3). This strategy can be applied to other problems expressed by Lagrangian in a similar way, such as the PN Lagrangian compact binaries and the PN Lagrangian n-body problems.

Numerical comparisons
In this section, we conduct numerical simulations for the PN circular restricted three-body problem with the new numerical scheme (25) proposed in the previous section. First, we check the effectiveness of the proposed integrators. Second, with this new approach, we further investigate on the dynamics of the PN circular restricted three-body problem. In particular, the dependence of the dynamics on the separation a between two primaries and the velocity of light c is discussed. Third, the dynamics between the exact system (12) and the 1PN truncated system (16)-(17) are compared. To obtain highly accurate numerical solutions, the coefficients of the eighth-and ninth-order Runge-Kutta-Fehlberg algorithm of variable stepsize [RKF8 (9)] are used to accomplish the classical scheme (23) and the proposed implicit scheme (25), respectively with the 1PN truncated system (16)- (17) and the exact implicit system (12). The tolerance of the local truncation error is set to 10 −14 during the variable stepsize procedure.
During the numerical simulation, we highlight that space and time are measured in the unit of a and 1/ω 0 , where ω 0 = G M/a 3 , correspondingly. Hence, the position (x, y) and the velocity (ẋ,ẏ) are measured in the unit of a and aω 0 , respectively. This setting is consistent with that in [5], where the authors implemented scaling transformations for these variables; in the present study, scaling transformations are not performed. In addition, we emphasize once again that geometrized units G = M = 1 are used.
We first use the RKF8(9) to integrate the circular restricted three-body problem in the Newtonian case, i.e., L 0 . By setting y(0) = 0,ẋ(0) = 0 and mass parameter μ 2 = 10 −3 , we consider the following four orbits [5]:ẏ(0) = 0.55, y(0) = 0.56,ẏ(0) = 0.60, andẏ(0) = 0.68, where x(0) is suitably fixed to satisfy C J = 3.07. In this case, the separation between the two primaries is a = 1. The four orbits are denoted by Orbits 1, 2, 3, and 4. We plot the Poincaré surface sections of the four orbits in Fig. 1. In this figure, Orbits 1 and 4 are regular, whereas Orbits 2 and 3 are chaotic. Furthermore, the Jacobian constant C J = 3.07 is preserved favorably in the magnitude of O(10 −12 ) for all the four orbits. All these results in the Newtonian are certainly consistent with those in [5].

Conservation of the Jacobian constant
The generalized Jacobian constant J in (10) should be preserved during the evolution of the exact system (12) that corresponds to the PN Lagrangian (3) but is only approximately preserved during the evolution of the 1PN truncated system (16)- (17). Therefore, we select the conservation of J to indicate the effectiveness of the proposed non-truncated strategy in Sect. 3. In this subsection, Orbit 1 is selected as an example. For selecting the distance a and the velocity of light c, we emphasize that if a 1 or c 1, then the difference between the non-truncated Eq. (12) and the 1PN truncated Eqs. (16)-(17) will be small due to the same limit as c → ∞ or a → ∞ [5]. Then, the 1PN truncated and the non-truncated strategies will show an excellent conservation for the Jacobian constant provided that c 1 or a 1. This phenomenon has been displayed in [6], where the authors showed a favorable Jacobian constant conservation by integrating the 1PN truncated motion equations and setting c = 10 4 1. However, for the case of a = O(1) and c = O(1), the Jacobian constant cannot be preserved well in the 1PN truncated system (16)- (17), because the truncated equations seriously deviate from the exact implicit motion equations. This point has also been confirmed in [6]. By contrast, considering that the fully implicit non-truncated Eq. (12) exactly corresponds to the generalized Jacobian constant J , a favorable conservation of the Jacobian constant for the proposed non-truncated strategy can be expected, even for the case of a = O(1) and c = O(1).
In the light of the above discussions, we are first concerned with the following two cases: Case 1 with a = 30 and c = 1 and Case 2 with a = 1 and c = 10; these two cases are nearly in the magnitude of O (1). The results are depicted in Fig. 2a, b. In both cases, the numerical solutions obtained from integrating the 1PN truncated Eqs. (16)- (17) only preserve the Jacobian constant in the accuracy of O(10 −3 ). For-tunately, the non-truncated strategy of integrating the exact non-truncated Eq. (12) can provide highly accurate numerical solutions that preserve the Jacobian constant in the accuracy of O (10 −12 ), which is in the same magnitude as the Newtonian case.
Moreover, we investigate the dependence of Jacobian constant errors on the two parameters c and a. The corresponding plots are exhibited in Fig. 2c, d. These plots indicate that the truncated approach preserves the Jacobian constant favorably once c or a enlarges, whereas the non-truncated approach preserves the Jacobian constant well for all possible values of c and a. This point confirms our theoretical analysis mentioned above . The results demonstrated in Fig. 2 clearly show that the proposed numerical scheme that integrates the non-truncated motion equations can provide more accurate numerical solutions than the commonly used approach that integrates the truncated motion equations.
Finally, we use the classical fourth-order Runge-Kutta method (RK4) to solve the 1PN truncated system (16)-(17) and the exact system (12) with the stepsize h = 0.01. Numerical solutions obtained from the exact system (12) solved by RKF8(9) are viewed as the reference solutions. The numerical results are plotted in Fig. 3. From this figure, it can be seen that the accuracy of the new scheme (25) is higher than the classical scheme (23) in the aspect of energy conservation and orbital configuration, provided that the two schemes share the same coefficients.

Dependence of the dynamics on c
We emphasize that the parameter c in the PN Lagrangian L of Eq. (3) denotes the readjusted velocity of light in our setting G = M = 1, rather than the real velocity of light. As is known, the PN approximations are essentially expansions in a small parameter, i.e., the PN parameter [6,40] where all physical quantities take their real values. Consequently, we should adjust the parameter c in accordance with our setting G = M = 1 to maintain the same value of ε and thus make the considered system physical and valid. For example, in the case of Sun-Earth coupling, the real PN parameter ε is expressed as which is coincident with the results in [6,41]. For other systems, the scaled parameter c will take different values. In particular, for the two black holes of the first observation event of gravitational wave GW150914 that has a total mass M 70M [42], the readjusted c changes during the orbital inspiral and merger, and subsequent final black hole ringdown. Let R s = 2G M/c 2 210 km denote the Schwarzschild radius of the binary black holes. Then, we should set c = 10 provided that a = 50R s and take c = 100 for the case of a = 5000R s . That is, during the evolution of the two black holes in GW150914, the scaled parameter c will be a variable changing from 1 to ∞. This point is also the reason that we investigate the dynamics of the PN circular restricted three-body problem by varying c in this section. Therefore, we set a = 1 and vary the value of c from 1 to 10 4 to compare the dynamics between the exact non-truncated system (12) and the 1PN truncated system (16)- (17).
We plot the orbital configuration of Orbit 3 in Fig. 4, where the orbits evolve in the three different systems, namely, the Newtonian system L 0 , the 1PN truncated system (16)- (17), and the exact non-truncated system (12) corresponding to the PN Lagrangian (3). In Fig. 4, for the case of c = 10, the difference is much more remarkable in orbits between the three systems (Fig. 4a) than in the case of c = 100 (Fig. 4b). This point is highly consistent with our theoretical analysis in Sect. 2.2 that a large value of c constantly denotes a small difference between the truncated motion Eqs. (16)-(17) and the exact motion Eq. (12). Thus, the 1PN truncation for the exact motion equations of the PN Lagrangian (3) may result in a notable deviation in orbital configuration from the underlying exact orbit once c is in a small magnitude (i.e., a strong PN effect).
To explore a detailed difference between orbits in the two systems (12) and (16)-(17), we use the regular Orbit 1 and the chaotic Orbit 3 as examples by varying c. The dependence of orbital errors on c is displayed in Fig. 5, where ΔX denotes the difference between solutions obtained from (12) and (16)-(17). In Fig. 5, for the regular orbit, the difference will decrease with the increase in c (Fig. 5a). However, for the chaotic orbit, the difference does not decrease with the increase in c (Fig. 5b); that is, the difference is in a large magnitude when even c takes a large value. This point has been theoretically explained in Sect. 2.2, where we empha- size that the difference between orbits in different systems can be estimated by O(c −4 ) only for the case in which the two orbits are regular, whereas the difference cannot be definitely estimated and is typically in a magnitude much larger than O(c −4 ), provided that at least one of the two orbits is chaotic. Figure 6 plots the dependence of the maximum radius R = x 2 + y 2 on c during the evolution of Orbit 1 in the 1PN truncated system (16)-(17) and the exact system (12).
In this figure, orbits evolving in the exact system (12) are all bound in the region R < 20. However, orbits evolving in the 1PN truncated system (16)-(17) may be in a large region because R ≈ 330 for c = 15.85. This phenomenon indicates the orbit expansion due to the loss of energy once the truncation is conducted for (16)- (17), which has been theoretically stressed in Sect. 2.2.
To gain further insight into the dynamical difference between the two systems (12) and (16) In this figure, the orbits in the two systems possess a large divergence for a small c (c = 10). Once c 1, the two plots (Fig. 7c, d) are significantly similar to the Newtonian case as presented in Fig. 1.
Furthermore, we use the fast Lyapunov indicator (FLI) [28,43,44], rather than the Lyapunov characteristic exponent [26][27][28]43], as a more sensitive indicator to determine chaos with reasonable computation costs. The FLIs successfully applied to detect chaos from order in systems of spinning compact binaries [45][46][47]. During the computation of FLIs, the two-particle method [28,43,48] is applied, and the initial separation is set to d(0) = 10 −8 as recommended in [27]. The integration time is set the same as in [5], i.e., T = 3000. Details on the computational techniques for FLIs can be referred to [43]. Here, we select a regular orbit (Orbit 1) and a chaotic orbit (Orbit 3) for comparison and set the threshold value of FLI to 7.5 to distinguish between regularity and chaoticity. The dependence of FLIs on c in the two systems is plotted in Fig. 8. In this figure, the dependence is nearly the same for the 1PN truncated and the exact systems whenever the orbit is regular or chaotic (in the Newtonian case). However, for certain orbits, their dynamics will essentially change in different systems. That is, the same initial values yield regular orbit in one system while they yield chaotic orbit in another. Two examples are illustrated in Fig. 9. In Fig. 9a, for the case of c = 44.67, Orbit 1 will be regular in the system of the exact motion Eq. (12), but it is chaotic in the system of the 1PN truncated motion Eqs. (16)- (17). By contrast, for the case of c = 50.12 (Fig. 9b), Orbit 3 is chaotic in (12), but is regular in (16)- (17), thereby indicating that the orbit evolving in the 1PN truncated system (16)-(17) may fully differ from that evolving in the exact system (12).
Based on the abovementioned results in this section, the 1PN truncated and the exact non-truncated systems have the same limit, i.e., the Newtonian system with the Lagrangian L 0 as c → ∞; however, the two systems can yield fully different results for several orbits in the aspects of the energy (or the Jacobian constant) conservation, the orbital configuration, and the regularity or chaoticity, especially in the case of small c. This result indicates that the 1PN truncation of (16)-(17) will introduce certain essential change in the dynamics for the Lagrangian (3), while the exact system (12) along with its non-truncated integration strategy, is exactly corresponding to the Lagrangian (3).

Dependence of the dynamics on a
This subsection mainly discusses the dependence of the dynamics on another parameter in the three-body problem, i.e., the separation a between two primaries. A similar work has been conducted in [5], where the authors only integrate the truncated motion equations. In this study, with the non-truncated strategy, we can integrate the exact motion equations. Thus, in this subsection we make the comparison between the 1PN truncated motion Eqs. (16)-(17) and the exact non-truncated motion Eq. (12). Throughout this section, we set c = 1. Figure 10 displays two different cases a = 30 and a = 100 for Orbit 1, which is evolving in the truncated and the exact non-truncated systems. It is noted that the evolution of this orbit is totally different in the two cases, i.e., a = 30 and a = 100. Further result on the dependence of the orbital error on a between the two systems is presented in Fig. 11, from which we can yield that chaos enables the orbital difference to be remarkable. Similar to Fig. 6, the maximum radius during the evolution of orbits for varying values of a are illustrated in Fig. 12. The result of the special case a = 80.12 for the 1PN truncated motion equations is not displayed in this figure because the maximum radius, in this case, is approximately R = 10 23.98 , which is much larger than the values in other cases and out of the scope of common bounded orbits in a finite time T =  (12) and the 1PN truncated system (16)- (17) system (3). However, the proposed non-truncated strategy can maintain all the orbits in a bounded region due to R < 20, that is, the non-truncated strategy preserves the conservation of the Lagrangian system (3). The Poincaré surface sections of the four orbits (Orbits 1, 2, 3, and 4) evolving in the two systems with a = 30, 100, 1000, and 10000 respectively, are demonstrated in Fig. 13. Two points can be determined from this figure. One point is that for a large a, the Poincaré surface sections of the two systems will be similar to the Newtonian case, as displayed in Fig. 1. It is due to the fact that a large a indicates a weak PN effect on the classical circular restricted three-body problem, and then the Newtonian part L 0 dominates the dynamical behavior. Another point is that, for the smallest case of a = 30, the Poincaré surface sections indicate distinct dynamics for Orbits 3 and 4 in the two systems. That is, the two orbits are chaotic in the exact system but are regular in the 1PN truncated system.
In addition to the Poincaré surface sections, we use the FLIs to accurately determine the chaos of Orbit 1 and 3 in Fig. 14. It immediately follows from Fig. 14 that for a 10 3.5 , the two orbits behave in the two systems similarly to that in the Newtonian system. This point supports once again that the two systems (12) and (16)-(17) tend to the same Newtonian case as a → ∞. The overall difference between the exact system (12) and the 1PN truncated system (16)- (17) is that the FLIs of Orbit 1 will be larger in the exact system than in the 1PN truncated system once a 10 2.3 . A detailed example that corresponds to Fig. 14 is presented in Fig. 15, in which the same initial values finally yield orbits of opposite dynamics. This point confirms that the exact system (12) and the 1PN truncated system (16)- (17) can generate totally different dynamics, although the later is an approximation to the former at the 1PN level.

Conclusions
In this paper, we presented a non-truncated strategy to exactly integrate the PN circular restricted three-body problem. The proposed non-truncated strategy combines the classical numerical methods (e.g., Runge-Kutta methods) with the fixed-point iteration. In comparison with the commonly used truncated strategy, the proposed approach can exactly integrate the implicit motion equations derived from the PN Lagrangian. This property yields the superiorities of the nontruncated strategy over the commonly used truncated strategy in several aspects. The first remarkable advantage of the new approach is its favorable conservation of the generalized Jacobian constant (or energy) for all possible choice of a and c. However, the truncated strategy can yield this favorable conservation only for the case of c 1 or a 1 (see [5,6]), thereby indicating that the PN contribution is in such a small magnitude that the PN Lagrangian is approximate to the Newtonian case. In addition, by integrating the exact implicit system via the non-truncated strategy, the dynamics of the PN Lagrangian presents qualitatively different behaviors in comparison with that of the 1PN truncated systems once a or c is in a small magnitude. This point confirms that the truncated strategy may destroy the dynamics of the original PN Lagrangian, although the truncated strategy is effective and convenient in the case of c 1 or a 1. In particular, the truncated strategy will also destroy the conservation of PN Lagrangian system, because the orbits are observed to be escaping from the bounded region in the numerical simulations of the 1PN truncated system. This destruction of the boundedness of orbits in conservative system has been observed in [5], where the authors attributed the phenomenon to the PN contributions. However, we consider that this destruction of boundedness is due to the energy dissipation of the 1PN truncation for the exact motion equa-tions. Moreover, the proposed non-truncated strategy can avoid this drawback and preserve the conservation of the PN Lagrangian well. This point also supports our numerical simulations using the non-truncated strategy to integrate the exact PN Lagrangian, where all orbits under consideration are observed to be bounded in a limited region.
In summary, the truncated strategy is common and convenient to use in the theoretical and numerical analyses for the PN circular restricted three-body problem. However, this approach may lead to several drawbacks, such as the nonconservation of the generalized Jacobian constant, the nonconservation of the PN Lagrangian (thus, the orbit escapes from a bounded region), and the destruction of orbital configurations and dynamical behaviors of order and chaos, especially in the case of strong PN contributions. Fortunately, the proposed non-truncated strategy can avoid these drawbacks by integrating the fully exact implicit motion equations that are directly derived from the PN Lagrangian via the Euler-Lagrange formula. This non-truncated strategy can also be extended to the PN Lagrangian systems of compact binaries in its numerical simulation and its equivalence with the PN Hamiltonian at the same PN level.