Contrasting the Implicit Method in Incoherent Lagrangian and the Correction Map Method in Hamiltonian

The equations of motion for a Lagrangian mainly refer to the acceleration equations, which can be obtained by the Euler--Lagrange equations. In the post-Newtonian Lagrangian form of general relativity, the Lagrangian systems can only maintain a certain post-Newtonian order and are incoherent Lagrangians since the higher-order terms are omitted. This truncation can cause some changes in the constant of motion. However, in celestial mechanics, Hamiltonians are more commonly used than Lagrangians. The conversion from Lagrangian to Hamiltonian can be achieved through the Legendre transformation. The coordinate momentum separable Hamiltonian can be computed by the symplectic algorithm, whereas the inseparable Hamiltonian can be used to compute the evolution of motion by the phase-space expansion method. Our recent work involves the design of a multi-factor correction map for the phase-space expansion method, known as the correction map method. In this paper, we compare the performance of the implicit algorithm in post-Newtonian Lagrangians and the correction map method in post-Newtonian Hamiltonians. Specifically, we investigate the extent to which both methods can uphold invariance of the motion's constants, such as energy conservation and angular momentum preservation. Ultimately, the results of numerical simulations demonstrate the superior performance of the correction map method, particularly with respect to angular momentum conservation.


Introduction
Compact binary systems, composed of neutron stars or black holes, etc., are of immense interest to experimental and theoretical researchers as sources of gravitational waves for broadband laser interferometers. The temporal progression of binary systems encompassing compact objects can be elucidated through the implementation of Einstein's equations of general relativity. Explicit symplectic integrators are supposed to be the ideal candidate with several benefits for the numerical simulations in these systems. They are designed to preserve the symplectic structure, which guarantees the precision and stability of numerical solutions over long time intervals. However, Einstein's equations of general relativity describe the motion of strong gravitational systems for which exact solutions are very difficult to obtain, but there have been some efforts to address this issue. For example, Xin Wu [X. Wu(2022)] developed the explicit symplectic integrators for Hamiltonian systems in curved spacetimes, particularly for black hole spacetimes. The papers [Zhou(2022), Zhou(2023)] discuss the motion of charged particles around the Schwarzschild black hole with an external magnetic field. Therefore, the orbital motion described in this paper is a two-body problem, involving the motion of a charged particle around a Schwarzschild black hole. Ying wang [Y. Wang(2021a), Y. Wang(2021b)] constructs the explicit symplectic integrators in general relativity, specifically for the Hamiltonian of Schwarzschild spacetime geometry. The integrators are useful for the long-term integration of N-body Hamiltonian systems and modeling the chaotic motion of charged particles around a black hole with an external magnetic field. Xin Wu ] discusses the construction of explicit symplectic integrators for Kerr black holes in general relativity. The authors introduce a time transformation function to the Hamiltonian of Kerr geometry to obtain a time-transformed Hamiltonian consisting of five splitting parts whose analytical solutions are explicit functions of the new coordinate time. Wei Sun ] proposes an explicit symplectic integrator for the Kerr spacetime geometry to simulate the nonintegrable dynamics of charged particles moving around the Kerr black hole in an external magnetic field. The algorithm shows good numerical performance and is used to study the dynamics of order and chaos of charged particles.
Despite the efficacy of this method, another useful and well-developed approach is the form of the post-Newtonian (PN) Lagrangian or Hamiltonian [G. Pan(2021)] approximation [Blanchet & Iyer(2003), Tanay(2021), G. Pan (2021)]. Arun et al. [Arun(2008)] investigated inspiralling compact binaries in quasi-elliptical orbits and provided a comprehensive analysis of the third post-Newtonian energy flux. Their study focused on understanding the energy loss due to gravitational radiation and its implications for compact binary systems. Tessmer and Schäfer [Tessme(2014)] studied the eccentric motion of spinning compact binaries. They examined the dynamics of these systems with non-circular orbits, considering the effects of spin and exploring the consequences of eccentricity on the gravitational wave signals emitted during inspiral. Hinder et al. [Hinder(2018)] developed an eccentric binary black hole waveform model by combining numerical relativity simulations with post-Newtonian theory. Their work aimed to accurately describe the complete inspiral-merger-ringdown phase of eccentric binary black hole systems, providing insights into the gravitational waveforms emitted during these events. Chattaraj et al. [Chattaraj(2022)] conducted high-accuracy comparisons between post-Newtonian theory and numerical relativity simulations, specifically focusing on eccentric binary black holes. They investigated the influence of higher modes on the waveforms and developed a model that incorporates eccentricity and accurately describes the inspiral, merger, and ringdown phases. Chowdhury and Khlopov [Chowdhury(2022)] studied an eccentric binary black hole system within the framework of post-Newtonian theory. Their research aimed to understand the behavior of binary black holes with non-circular orbits, providing insights into the dynamics and gravitational wave emissions of eccentric binary systems. These approximations provide high-precision theoretical templates of gravitational waveforms, although their higher-order terms are truncated, which affects their equivalence [X. Wu(2015a), X. Wu,(2015b), L. Huang(2016].
The choice of approximation and the selection of the algorithm becomes crucial to ensuring an accurate and effective description of the trajectory evolution of compact binaries systems and matching corresponding gravitational waveforms.
The PN Lagrangian equations of motion are derived from the Euler-Lagrangian equations of a PN Lagrangian formulation, denoted as L(r, v). By calculating the partial derivative of L with respect to velocity, we obtain the generalized momentum p = ∂L/∂v. Similarly, the acceleration equations a = f (r, v, a) = ∂L/∂r can be derived and we can obtain a coherent Lagrangian [Li(2021a), Li(2019b), Li(2021b)]. By limiting the inclusion of accelerations up to a certain PN order in the Lagrangian, the accelerations a in the function f will be modified to a * ; it only has lower-order terms, i.e., a = f (r, v, a * ).The acceleration equations become incoherent, due to the higher-order PN terms disappearing, leading to a loss of some values of the constants of motion during subsequent evolution. The same problem occurs with the post-Newtonian Hamiltonian form. The error of the constant of motion can be used as an indicator to test the performance of different algorithms in both approximate forms.
Various algorithms are available for the calculation of post-Newtonian Lagrangian quantities. For example, in optimizing the fifth-order Runge-Kutta method as a high-precision integrator, Zhong [Zhong(2010)] employed corrections to all integrals within the conservative 3PN order Hamiltonian. Tsang [Tsang(2015)] introduced an implicit symplectic integrator that accounts for 2.5PN gravitational radiation reaction terms in the Newtonian two-body problem. This approach effectively captures the effects of radiation reactions. Lubich [Lubich(2010)] devised an explicit and implicit mixed symplectic integration technique that facilitates the splitting of orbital and spin contributions. By employing this approach, the dynamics of both orbital and spin variables can be accurately simulated. Zhong [Zhong et al(2010)] proposed fourth-order canonical explicit and implicit mixed symplectic methods. These methods offer improved accuracy and stability in the computation of post-Newtonian quantities. Seyrich [Seyrich(2013)] developed Gauss Runge-Kutta implicit canonical symplectic schemes that preserve the structural properties of the system. These schemes ensure long-term numerical stability and accuracy. These algorithms, with their distinct methodologies, contribute to advancing the computation of post-Newtonian Lagrangian quantities, addressing specific aspects such as precision, radiation reaction, spin contributions, stability, and structural preservation.
Regarding the post-Newtonian Hamiltonian, the phase-space expansion method [Pihajoki(2015), Li(2017), Li(2019a)] is a usable algorithm. The Hamiltonian lacks separability and does not possess a coordinate momentum or multiple integrable splitting components. Pihajoki [Pihajoki(2015)] extended the phase space variables by copying the coordinates and momenta. We achieved a Hamiltonian splitting form so that the explicit leapfrog algorithms become available. The permutation map of momentum was designed to suppress the interaction of the original and extended variables. Liu [Liu et al.(2016)] devised a sequential mapping of coordinate and momentum permutations and constructed fourth-order phase-space expansion explicit method compositions of two triple products of the usual second-order leapfrog. These algorithms suffer a clear failure when calculating the chaotic orbits of celestial systems. The interactions between the original variables and the extended one become increasingly strong and show considerably different values, whereas they are supposed to be equivalent. Midpoint and correction maps [Luo et al.(2017), Luo et al.(2021)] have been proposed to ensure the equivalence of the original variables and the copy one. Recently, we proposed a multi-factor correction map that yields a higher accuracy of the phase-space expansion method without significant computational resource increases [Luo et al.(2022)]. This paper aims to design a multi-factor correction map for post-Newton Hamiltonian and examine its performance.
This article is divided into several sections. In Section 2, we revisit the Lagrangian and Hamiltonian equations of motion for compact binary systems within the post-Newtonian (PN) approximation. Section 3 presents the introduction of the phase-space expansion method and the development of a correction map for the post-Newtonian Hamiltonian. In Section 4, we conduct a comparative analysis of the accuracy of numerical solutions obtained using the implicit midpoint method in the computation of the post-Newtonian Lagrangian and the correction map method for the post-Newtonian Hamiltonian. Finally, in Section 5, we conclude

PN Lagrangian and Hamiltonian in Compact Binary
Let us consider a compact binary system governed by a PN Lagrangian L(r, v) up to the m-th order, where r and v represent the position and velocity vectors, respectively. The Euler-Lagrangian equation is given by Equation 1, where the generalized momentum p is defined by Equation 2. The expression for p is given by Equation 2 and represents a nonlinear algebraic equation of v. dp dt = ∂L ∂r . (1) Here p = ∂L ∂v ; (2) We note that Equations 2 and 3 are differential equations, and that (r, p) are treated as the integration variables, whereas v is not. However, we can substitute Equation 2 into Equation 3 to obtain the corresponding acceleration equation, given by Equation 4. Here, aN , a1P N , a2P N , . . . , amP N correspond to the Newtonian term, the 1st, 2nd post-Newtonian-order term to the m-th post-Newtonian-order contributions for the accelerations.
When considering only the m-th PN order term in Equation 4, all terms higher than the m-th PN order are truncated. Consequently, Equation 4 does not align with the PN Lagrangian L, and Equations 3 and 4 are treated as incoherent PN equations of motion in the Lagrangian L. However, when utilizing Equations 3 and 4, the variables (r, v) can be used as a set of integration variables instead of the variables (r, p). Nevertheless, this approach does not fully maintain constants of motion, such as the energy integral expressed by In this paper, L is the dimensionless post-Newtonian (PN) Lagrangian formulation for compact binaries. The evolution of binaries can be given by the expression: LN and L1P N denote the non-relativistic and 1PN contributions to the Lagrangian, respectively. For simplicity, higher-order terms are not considered. The non-relativistic part is expressed as: whereas the 1PN part is given by [Blanchet & Iyer(2003)]: Here η = µ/M is the dimensionless mass parameter. The reduced mass, µ, is defined as M1M2/M = β(1 + β) −2 , and β = M1/M2 is the mass ratio, where M1 and M2 represent the masses of the two bodies constituting the binary system and the total mass is denoted as M = M1 + M2. Additionally, c is the speed of light and G represents the constant of gravity given in natural units with c = G = 1. c is retained in some of the latter equations, and it can be ignored in the actual calculation.
The equations for the evolution of the system can be derived from the Lagrangian formulation. According to Equations 1 and 2, the equation for the evolution of the momentum can be written as: where r is the separation vector between the two masses. The non-relativistic and 1PN contributions to this equation are, respectively, expressed in the first and second terms. The expression for the generalized momentum p till 1pN is given in terms of the velocity v as: The value of momentum can be obtained from Equation 10 once the velocity is known and vice versa. However, velocity needs to be solved iteratively since it cannot be obtained directly from the Lagrangian. The first post-Newtonian relative acceleration equation is used to determine the value of velocity, given by: The sub-terms are Here, aN and a1P N describe the non-relativistic and 1PN contributions to the acceleration. Using Equations 6 and 10, the energy integral in Equation 5 can be expressed as.
In summary, the Lagrangian formulation of the evolution of binaries provides a mathematical framework for studying their motion. The momentum and velocity of the system can be determined from the equations derived from the Lagrangian, which include non-relativistic and relativistic contributions to the acceleration. The dimensionless PN Lagrangian offers valuable insights for the dynamics system, enabling the study of gravitational wave emission caused by binary systems.
With Equations 3 and 11, we can obtain the numerical solution (r · v) by using the fourth-order implicit midpoint method (IM4).
PN Hamiltonian form H can be derived through Lagrangian L using the Legendre transformation, Then we obtain the 1PN Hamiltonian, In order to compare the effect of higher-order PN terms on the error in the constants of motion, we introduce the 2PN post-Newton Hamiltonian, The expressions for the sub-terms in Hamiltonians 16 and 17 are, respectively, given by Due to the disappearance of higher-order terms, H and H * are approximately equal to E, and not strictly equivalent. The integrators used in the Hamiltonian H and H * will be described in the next section.

Phase-Space Expansion Method with a Multi-Factors Correction Map
Since neither the Hamiltonian H nor H * can be separated into multiple integrable parts, the symplectic leapfrog method cannot be applied directly to these Hamiltonians unless they are suitably modified to a splitting form. An effective approach to solving this problem is the phase-space expansion method. Pihajoki [Pihajoki(2015)] introduced a new pair of canonical and conjugate variables ( r, p) from the original variables (r, p). This doubles the phase-space variables, (r, r, p, p) and constructs a new Hamiltonian H(r, r, p, p) using two identical Hamiltonians H1 and H2: H(r, r, p, p) = H1(r, p) + H2( r, p).
where both H1 and H2 should be equal to the original Hamiltonian H. The new Hamiltonian H already exhibits two integrable components. A conventional second-order leapfrog algorithm can be employed for its integration: where h represents the time step, H1 and H2 are Hamiltonian operators. The code corresponding to integration 22 from nth to (n + 1)th step is It can be seen that the calculation of the numerical solution (r, p) of H2 requires the numerical solution ( r, p) of H1, and vice versa, so there is an energy exchange between H1 and H2, and even if the initial conditions are the same, H1 and H2 will become unequal in the later evolution unless the errors are constant equal to 0. To submit the accuracy, we construct a fourth-order algorithm using Yoshida's triplet product The time coefficients λ1, λ2, and λ3 are identical to those presented in [Yoshida(1990)] and are set to λ1 = λ3 = 1/(2 − 2 1/3 ) and λ2 = 1 − 2λ1. Algorithm A4 is utilized to obtain a set of numerical solutions (r, r, p, p). It is crucial to note that the original variables (r, p) and their counterparts ( r, p) are intended to be identical at each integration step, but in reality, they exhibit discrepancies. The interaction between the solutions (r, p) of H1 and ( r, p) of H2 leads to their divergence over time. To ensure the equivalence of the original variables and their copy one, Pihajoki [Pihajoki (2015)] proposed a momentum permutation map, whereas Liu [Liu et al.(2016)] proposed coordinate and momentum permutation maps. These applications were successful in several examples, but not in chaotic orbits. However, our previous work [Luo et al.(2021), Luo et al.(2022)] proposed the manifold corrections map, which effectively overcame the challenge. Unlike that in the original paper [Luo et al.(2021), Luo et al.(2022)], the correction map for the 1PN Hamiltonian is Here, the momentum scaling factor γ and the coordinate scaling factor α are incorporated into M1P N and can be obtained by solving the following formulations: Equation 26 provides α = 2(p 2 + p 2 ) (p+ p) 2 , and Newton's method is used to obtain γ from Equation 27. Then, the fourth-order phase-space expansion method with multi-factor correction map for 1PN Hamiltonian H is established as Similarly, for the 2PN Hamiltonian H * , the aforementioned steps are applicable. The correction map M2P N for the Hamiltonian H * follows the same structure as M1P N .
However, the solution for γ is replaced by the following equation = V (r, r) + H1P N (r, r, p, p) + H2P N (r, r, p, p) 2 .
For convenience, such algorithms are referred to as the correction map method. The correction map method for 2PN Hamiltonian H * is set up as The A4 algorithm is treated as an explicit symplectic method serving to the new Hamiltonian H, ensuring effective preservation of the energy of H, i.e., ∆ H = ∆H1 + ∆H2 ≈ 0. The error evolution of H1 and H2 shows a clear time-axis symmetry, as demonstrated in Figure 1. Specifically, if H1 calculates more energy than the initial energy, H2 will calculate less, and vice versa. Taking advantage of this symmetry, we designed a manifold correction mapping approach to optimize the performance of the A4 algorithm. The M1P N and M2P N corrections imposed on the solutions of A4 serve three primary purposes. The A4 algorithm serves multiple purposes in relation to the new Hamiltonian H. Firstly, it ensures that H1 is equal to H2 to prevent energy discrepancies that could impede the availability of numerical solutions. Secondly, it maintains the constancy of H after the correction, effectively suppressing the growth of energy errors. Thirdly, it reduces the energy deviation of each subterm of H from half of the corresponding subterm of H through the correction process. As the A4 algorithm is an explicit symplectic method serving the new Hamiltonian H, it accurately calculates the total energy as well as the energy of each individual subterm in H; thus, the algorithm altering these energies is not desirable, as it may weaken the algorithm's stability and precision.
In Section 4, we will set initial values and perform numerical simulations of post-Newtonian Lagrangian and post-Newtonian Hamiltonian to compare the differences between the algorithms in terms of maintaining the constants of motion.

Numerical Simulation
This section showcases the outcomes of our numerical simulations, where we compare the post-Newtonian Lagrangian and post-Newtonian Hamiltonian algorithms in maintaining the constants of motion. To this end, we set initial values for a specific orbit, named orbit 1, with initial conditions (β; r, v) = ( 5 4 ; 10, 0, 0, 0, 0.52, 0). The initial value of the momentum p in the post-Newtonian Hamiltonian is obtained from Equation (10). We use the fourth-order implicit midpoint method (IM4) to calculate the 1PN Lagrangian, whereas the algorithms CM1P N and CM2P N are used to calculate the Hamiltonian H and H * , respectively. We take a fixed step size of h = 1 and plot the energy errors in Figure 2a,b. We observe that the CM1P N algorithm designed for the Hamiltonian H has significantly better accuracy in terms of energy error compared to IM4. However, the accuracy of CM1P N drops considerably in the H * error behavior, as expected due to the vanishing of the 2PN term, whereas the error in CM2P N is at an order of 10 −8 . Finally, Figure 2c shows that CM2P N performs the best in terms of accuracy and long-term stability compared to CM1P N and IM4.
Aside from ensuring energy conservation, we also track the preservation of orbital angular momentum L = r × p =[1 + 1 c 2 ( 1−3η 2 v 2 + 3+η r )]r × v, and examine its error as another performance metric for the algorithms. Figure 3 depicts the angular momentum errors ∆ L = L − L0 with L = | L|, where L0 represents the initial value. We deduce that the performance of the IM4 algorithm in terms of angular momentum error is similar to its performance in energy error, with the worst accuracy but good long-term stability. Conversely, both CM1P N and CM2P N exhibit exceptional performance in terms of angular momentum error, with very little difference between them and significantly superior to IM4. However, there is a noticeable error growth in CM1P N and CM2P N . Our simulations show that the algorithm CM1P N in the 1PN Hamiltonian has a significant advantage over the IM4 algorithm in maintaining the conservation of orbital angular momentum and a small accuracy advantage in maintaining energy integrals. The algorithm CM2P N in the 2PN Hamiltonian also has a significant advantage in maintaining angular momentum, while being comparable to CM1P N , with some improvement in the accuracy of the energy.
To validate the aforementioned conclusion, additional numerical simulations will be conducted in a different orbit, referred to as orbit 2. The initial conditions for orbit 2 are set as (β; r, v) = ( 5 4 ; 10, 0, 0, 0, 0.52, 0). In Figure 4, three categories of energy error in orbit 2 will be depicted as follows: (a) Energy error analysis of orbit 2 for IM4, CM1P N , and CM2P N will be presented in Figure 4a. (b) Figure 4b will display the energy error analysis of orbit 2 specifically for IM4. (c) The energy error analysis of orbit 2, focusing on IM4, will be illustrated in Figure 4c.
It is evident from the figures that the performance of IM4, CM1P N , and CM2P N in orbit 2 closely resembles that of orbit 1. IM4 continues to exhibit the highest error, CM1P N demonstrates a widening gap with IM4, and CM2P N remains the most advanced in terms of accuracy.
Turning to the angular momentum errors depicted in Figure 5, it is observed that there is no significant improvement for IM4, which still exhibits considerable deviation compared to the first-order post-Newtonian approximation, CM1P N . Furthermore, the inclusion of the second-order post-Newtonian term in CM2P N does not contribute significantly to reducing the angular momentum error.
Summarizing the findings from the numerical simulations conducted for both orbit 1 and orbit 2, we can conclude that CM1P N and CM2P N perform better for the calculation of the post-Newtonian approximation to the Hamiltonian. They exhibit a slight advantage in terms of energy error while demonstrating a notably superior accuracy in the calculation of angular momentum.

Summary
The exact equations of motion for a post-Newtonian Lagrangian formalism are the Euler-Lagrange equations, which consist of a coherent Lagrangian without any truncated terms. However, when the post-Newtonian Lagrangian form of general relativity maintains only a certain post-Newtonian order, it is referred to as the incoherent Lagrangian, with higher-order terms of the acceleration truncated. Incoherent Lagrangian can be numerically simulated using the Runge-Kutta method and implicit algorithms. Therefore, in the incoherent Lagrangian, motion constants such as energy integrals are only approximately conserved. The retention of the Hamiltonian in a certain post-Newtonian (PN) order leads to a high-order truncation problem. In addition, if the Hamiltonian is separable, symplectic algorithms can be used, which provide excellent performance. For the case where the Hamiltonian is inseparable, symplectic-like algorithms such as the phase-space expansion method with correction map, namely correction map method, can be used. The phase-space expansion method with correction map is referred to as the correction map method, which utilizes the symmetry of energy errors in H1 and H2 in the new Hamiltonian H to improve the accuracy and stability of the algorithm.
A comparison was made between the performance of the implicit midpoint method in the incoherent ) and its initial value (H(0)). (b) The energy error of H * , denoted as ∆H * , is determined as the absolute difference between the value of the Hamiltonian H * at time t (H * (t)) and its initial value (H * (0)). (c) The energy error of E, denoted as ∆E * , is computed as the absolute difference between the value of E at time t (E(t)) and its initial value (E(0)). The algorithm IM4 is represented by a solid black line, whereas C1PN and C2PN are indicated by dashed red and blue lines, respectively. The performance of each algorithm in orbit 2 is similar to that of orbit 1. Lagrangian and the correction map method in the PN Hamiltonian. Under the 1PN Newtonian approximation, the correction map method performed better in terms of energy error, exhibiting higher accuracy and comparable stability. On the other hand, with regards to angular momentum error, the correction map method was significantly higher, reaching an order of 10 −11 , whereas the implicit midpoint method was only 10 −1 . Similarly, under the 2PN Newtonian Hamiltonian, the manifold correction mapping method further improved the accuracy of energy error, but there was no noticeable impact on the angular momentum error.
In conclusion, we compared the implicit midpoint methods for solving the equations of motion in post-Newtonian Lagrangians and the correction map method for PN Hamiltonians and investigate the extent to which both methods can uphold invariance of the motion's constants, such as energy conservation and angular momentum preservation. Ultimately, the results of numerical simulations demonstrate the superior performance of the correction map method, particularly with respect to angular momentum conservation. Compared to incoherent Lagrangian, we recommend using the manifold correction map method for the Hamiltonian of compact binaries as a numerical tool.