Dissipated Correction Map Method with Trapezoidal Rule for the Simulations of Gravitational Waves from Spinning Compact Binary

The correction map method means extended phase-space algorithm with correction map. In our research, we have developed a correction map method, specifically the dissipated correction map method with trapezoidal rule, for numerical simulations of gravitational waves from spinning compact binary systems. This new correction map method, denoted as $CM3$, has shown remarkable performance in various simulation results, such as phase space distance, dissipated energy error, and gravitational waveform, closely resembling the high-order precision implicit Gaussian algorithm. When compared to the previously used midpoint map which denoted as $C_2$, the $CM3$ consistently exhibits a closer alignment with the highly accurate Gaussian algorithm in waveform evolution and orbital trajectory analysis. Through detailed comparisons and analyses, it is evident that $CM3$ outperforms other algorithms, including $CM2$ and $C_2$ mentioned in this paper, in terms of accuracy and precision in simulating spinning compact binary systems. The incorporation of the trapezoidal rule and the optimization with a scale factor $\gamma$ have significantly enhanced the performance of $CM3$, making it a promising method for future numerical simulations in astrophysics. With the groundbreaking detection of gravitational waves by the LIGO/VIRGO collaboration, interest in this research domain has soared. Our work contributes valuable insights for the application of matched filtering techniques in the analysis of gravitational wave signals, enhancing the precision and reliability of these detection.


Introduction
In our recent work, our group has developed correction map for the extended phase-space method for the compact spinning binary systems.Extended phase-space algorithms are proposed by [Pihajoki:2015], which replicate the origin variables of position and momentum to generate corresponding copy variables, thereby expanding the phase space representation.All of these processes is particularly aimed at enabling the application of second-order explicit leapfrog algorithms to inherently inseparable Hamiltonian systems, exemplified by the complex dynamics of spinning compact binary systems modeled through post-Newtonian approximations as detailed in [Blanchet:2002mb, Tanay2021, Zotos 2019, Li 2020].By employing such extended phase-space techniques, researchers gain access to numerical solutions not only for the original position and momentum but also for the duplicated variables.In an ideal scenario, the original and replicated variables should maintain strict equality throughout the simulation.However, practical computations reveal that temporal discrepancies tend to accumulate between these paired variables due to their inherent interactions within the system.To counteract this unwanted drift, [Pihajoki:2015] proposes the implementation of a momentum permutation map, which systematically exchanges the values of the original and copied momenta at each time-step, reducing the difference between original variables and copy ones.Therefore, the design of the mapping matrix has a great impact on both the accuracy and stability of the extended phase-space method.Building upon the foundational work laid out in [Pihajoki:2015], [Liu:2016MNRAS] further refined these developments by devising the coordinate and momenta sequent permutation maps for a fourth-order extended phase-space explicit algorithm.This advanced algorithm was ingeniously crafted using a combination of two Yoshida's triple products derived from the second-order leapfrog algorithm, as initially presented by [Yoshida 1990PRD].The purpose was to enhance the accuracy in terms of energy error behavior during simulations.
However, despite these improvements, the algorithm exhibited significant shortcomings when applied to numerical simulations of chaotic orbits.Over time, the disparities between the original and extended variables grew due to their mutual interactions, which posed a substantial challenge.While these differences were negligible for non-chaotic or regular orbits, they could be exacerbated in numerically sensitive chaotic systems, leading to a detrimental feedback loop.To address this critical issue, [Luo:2017ApJ] proposed the midpoint map.This method effectively ensures that the original and extended variables remain strictly identical throughout the simulation process.Moreover, it necessitates the use of only a single Yoshida's triple product to construct a fourth-order algorithm, thereby achieving a twofold increase in computational efficiency without compromising the accuracy of results.Moreover, [pan2021extended] successfully employed the midpoint map in the context of coherent post-Newtonian Euler-Lagrange equations, achieving commendable computational outcomes.Further underscoring the efficacy of the midpoint map, [hu2019novel] demonstrated its superior performance relative to various alternative algorithms in their recent study.In order to further improve the accuracy of this algorithm, we design the correction map.The introduction of manifold correction technology into this framework significantly enhances the precision and stability of the algorithms used to simulate such celestial dynamical behaviors.The manifold correction ensures a more accurate preservation of the intrinsic geometric structure of physical systems during long-term evolution, thus reducing cumulative errors that might arise from conventional numerical integration methods and maximizes the energy conservation of the conservative system [Wu 2007APR, Ma:2008ApJ, Wang 2018AAS].And the correction map mitigates the divergence issues in chaotic orbits.
However, it is important to acknowledge that in reality compact spinning binary stars do indeed experience gravitational dissipation.When gravitational dissipation is taken into account, these systems cease to be conservative and instead continuously emit gravitational waves as they lose energy over time-a prediction rooted in Einstein's General Theory of Relativity.This theoretical prediction was spectacularly confirmed by direct observation in 2016 with the LIGO experiment ([LIGOScientific:2016aoc]), which detected gravitational waves emitted from the coalescence of compact binaries, thereby validating a cornerstone of relativity theory and inaugurating the new era of gravitational wave astronomy.It also arouses the enthusiasm of researchers to study the dynamic evolution of compact celestial bodies.[wu2022explicit] presents such integrators tailored to black hole spacetimes, providing a robust framework for studying the dynamics of particles in these highly curved regions.[zhou2022note] delve into the specific scenario of charged particles orbiting a Schwarzschild black hole under the influence of an external magnetic field.These studies elucidate the intricate interplay between gravity, charge, and the magnetic field, thereby enriching our understanding of the two-body problem where a charged particle interacts with the black hole.Ying Wang's work in the papers[wang2021construction] further extends the scope of explicit symplectic integrators to general relativity, focusing on the Hamiltonian associated with Schwarzschild spacetime geometry.The resulting integrators demonstrate their efficacy in accurately simulating the long-term behavior of N-body Hamiltonian systems, as well as in unraveling the chaotic motion of charged particles in the presence of both a black hole and an external magnetic field.[wu2021construction] addresses the more complex case of Kerr black holes, constructing explicit symplectic integrators within the context of general relativity.To achieve this, a time transformation is ingeniously applied to the Hamiltonian of Kerr geometry, yielding a transformed Hamiltonian composed of five separable components, each with analytically tractable solutions expressed in terms of a new coordinate time.This innovative approach enhances the efficiency and accuracy of numerical simulations for Kerr black hole systems.Wei Sun's contribution in paper[sun2021applying] introduces an explicit symplectic integrator tailored to the Kerr spacetime geometry, aimed at accurately capturing the non-integrable dynamics of charged particles orbiting a Kerr black hole under the influence of an external magnetic field.
The orbital motion of these compact objects is more dissipative than that of other small objects and the gravitational waves are easier to detect and identify.When the gravitational dissipation term is cosidered, in this thesis we readjust the scheme of the correction map to better serve the numerical simulation of gravitational waveforms.As a comparison, we will introduce an implicit algorithm, there are many kinds of implicit algorithms, for example, [Suzuki 1990]'s fourth-order composition in this context is realized as a composite of five individual second-order integrators.Expanding upon this work, [Zhong:2010PRD] incorporated conjugate canonical spin variables (as presented by [Wu 2010APS]) to develop fourth-order canonical explicit-implicit mixed-symplectic methods.In these methods, separable Hamiltonians are handled using the second-order explicit leapfrog algorithm, while nonseparable Hamiltonians are computed using the secondorder implicit midpoint method.Subsequent elaborations and applications of this construction were further explored in [Wu 2011] and [Mei 2013EPJC].[Seyrich:2013PRD] contributed to the field by proposing Gauss Runge-Kutta implicit canonical symplectic schemes, which inherently preserve the underlying structure of the system being modeled.In addition, [Tsang 2015APJL] advanced the development of implicit slimplectic integrators, tailored for the integration of general nonconservative systems, such as a Newtonian two-body problem subject to 2.5PN gravitational radiation reaction terms.These integrators demonstrate the continued evolution and refinement of numerical techniques for accurately simulating complex astrophysical phenomena.
We include the 4-stage Gauss Runge-Kutta implicit canonical symplectic algorithm([Seyrich:2013PRD]) as a numerical simulation algorithm tool, whose high accuracy numerical solutions can be treated as true solutions.Implicit algorithms are alternative solutions for computing binary star gravitational waves because they do not have the restriction of the use of separable or inseparable Hamiltonian.Nevertheless, implicit algorithms have multiple iterations that consume computational resources heavily, as well as convergence problems of the iterations leading to unavailability of numerical solutions.
The manuscript is structured as follows.In section 2, we present the Hamiltonian formulation for spinning compact binary stars, considering the 2.5PN dissipative term, and present the corresponding gravitational wave equations for the two polarization states.Then the dissipated correction map algorithm with the trapezoidal rule will be introduced, detailing its implementation and underlying principles.In section 3, Initial values for various parameters are specified, and the numerical solutions for spinning compact binary stars are computed using the dissipated correction map algorithm.The accuracy of this approach is assessed via phase space distance, dissipated energy error, and gravitational waveform comparison, utilizing the numerical solutions from the 4-stage implicit Gaussian algorithm as the reference truth.Finally, conclusions are drawn based on the results and analyses presented throughout the paper in the section 4.
2 Spinning Compact Binary and Gravitational Wave 2.1 The Hamiltonian of spinning compact binary Before introducing the algorithm, we wish to provide a concise description of the Hamiltonian for a spinaligned compact binary system composed of two objects such as black hole or neutron star, along with the associated formula for gravitational dissipation term.We assume that the individual masses of the two bodies are denoted by m1 and m2, with their total mass given by M = m1 + m2, the reduced mass defined as µ = m 1 m 2 M , the mass ratio β = m 1 m 2 , and the symmetric mass ratio η = m 1 m 2 M 2 = β (1+β) 2 .The dimensionless coordinate r is expressed in units of M , the momentum p is measured in units of the reduced mass, and t is denoted as time.
We define the vector r to represent the relative separation between the two bodies, with its magnitude representing the distance.The unit vector pointing along r is then given by n = r |r| .The provided equations describe the Hamiltonian H of a binary system, composed of two interacting bodies, taking into account various post-Newtonian (PN) corrections and spin-orbit and spin-spin interactions.The Hamiltonian is expressed as a sum of several sub-Hamiltonians: (1) Each of these components represents a distinct aspect of the system's dynamics and is respectively written as: HSOSS = H1.5SO+ H2SS denote Spin-orbit and spin-spin interactions.This term captures the influence of the bodies' intrinsic spins on their orbital motion.It is composed of two parts: where S is given by: and S0 assumes the form: The positional vector r and its associated momentum p obey the Hamiltonian canonical equations, which govern their temporal evolution as follows: (10) These equations embody the fundamental dynamics of the non-spin components of the system.Simultaneously, the spin vectors Sj experience time evolution according to: wherein the cross product × encapsulates the intrinsic rotational character of the spin dynamics.Collectively, these expressions define the comprehensive time-dependent behavior of the system as prescribed by the Hamiltonian H.
Observing equation 11, it becomes evident that neither of the spin variables exhibits the canonical conjugacy property.To rectify this, [Wu 2010APS] capitalized on the conserved nature of the individual spin magnitudes and introduced a novel set of generalized coordinates θj and the corresponding generalized momenta ξj.This transformation allows for the unit spin vectors to be recast in the following form: The Hamiltonian, as expressed by equation 1, now manifests a structure with ten distinct, canonically conjugate phase space variables, constituted by the tuple (r, θ1, θ2; p, ξ1, ξ2).This reconfiguration ensures that each coordinate variable is paired with its proper momentum counterpart, thereby restoring the fundamental symplectic structure inherent in Hamiltonian mechanics.Consequently, under this newfound parametrization, the Hamiltonian assumes the appellation of Before considering gravitational dissipation The constants of motion mentioned in this systems include the total energy E = H, total angular moment vector J = S1 + S2 + r × p and the spin lengths S 2 j = S 2 j .These constants of motion play crucial roles in constraining the behavior of the binary system and serve as essential diagnostics when numerically integrating the equations of motion derived from the Hamiltonian (1).The sophisticated symplectic integration schemes mentioned above, such as those developed by [Lubich:2010mj], [Zhong:2010PRD], [Seyrich:2013PRD], and [Tsang 2015APJL], are designed to efficiently and accurately simulate the dynamics of such systems while preserving these fundamental conservation laws.

Gravitational radiation reaction
Before discussing the gravitational dissipation of binary stars, we write down the second-order post-Newtonian approximation Lagrangian corresponding to the Hamiltonian 1 where the terms LSOSS consisted by L1.5SO and L2SS are provided in [Hartl:2004xr], the Langrangian 14 and Hamiltonian 1 can be transformed into each other by Legendre transformation The Hamiltonian 1 that governs a binary system is conserved, unless dissipative effects from gravitational radiation at the 2.5PN order are taken into account, which render the system non-conservative.Conventionally, such dissipative terms are not directly integrated into the Lagrangian expression of the equation 14; instead, they should be appended to the equations of motion directly.With doubling the degrees of freedom, r → (r, r), [Tsang 2015APJL] introduce the dissipative terms in Lagrangian form: Here r+ = (r + r)/2 and r− = r − r.Within the adopted unit system, the description of the radiative reaction force, as presented in the works of [Galley.PhysRevD.79.124027] and [Galley.galley2012radiation], is given by Analogous to the conventional Lagrangian L, which yields the conservative Euler-Lagrange equations of motion, the nonconservative Lagrangian Λ gives rise to a set of nonconservative Euler-Lagrange equations governing the system's dynamics.
The abbreviation "PL" denotes the physical limit, a condition in which the two variables r and r are required to be identical (i.e., r = r and ṙ = ṙ, or equivalently r+ = r, r− = 0, and ṙ+ = ṙ) after the necessary derivatives have been computed in equation 24.It is crucial to note that enforcing the PL in equation 22 would render it devoid of practical utility.When the PL is taken into account, the extended coordinate r and its time derivative ṙ cease to feature in the equations.
Although in the PL, r and r are equated, it is important to emphasize that r is distinct from the extended coordinate r which will be introduced in following section and used in the construction of extended phasespace algorithms.Similarly, while r and r are set equal at every integration step, they do not share the same conceptual basis as r and r in the PL.
Employing the Legendre transformation, [Galley.galley2013classical]derived a New Hamiltonian from the newly formulated non-conservative Lagrangian, It is worth mentioning that the momenta p and p are the standard conjugate momenta associated with the conservative scenario discussed earlier.Furthermore, the term LRR appearing in equation 25 is derived by substituting p and p for ṙ and ṙ, respectively, in equation 25.The Hamiltonian's canonical equations governing the dynamics of this nonconservative system Γ can be expressed as: Indeed, the dissipative terms are conspicuously absent in equation 26, given that LRR in equation 31 is explicitly independent of the momentum component p−.Analogously to the conservative scenario, the dissipative effects manifest solely in the acceleration equation 29.
Upon the emission of gravitational radiation, the total time derivative of the Hamiltonian H fails to vanish, instead assuming the form: Knowing the Hamiltonian canonical equations under the effect of gravitational dissipation, we can obtain numerical solutions for the evolution of binary orbits over time, and these numerical solutions will become the variables used in the calculation of gravitational waves, of course, the actual observed waveform depends on the direction of the wave source and the observer.Assuming an observer located in the xoz plane, let p denote the direction along the intersection line of the orbital plane with the horizon.Define q = N × p, where N represents the orientation of the observer.The two polarization states of a gravitational wave are given by The indices i, j = 1, 2, 3 correspond to the x, y, and z components, respectively, with repeated indices denoting Einstein's summation convention.The 2PN approximation of the waveform hij ([will1996gravitational]) expands as Here, D represents the distance between the observer and the wave source.The superscripts, such as P 1.5 , indicate the effective post-Newtonian (PN) order, while the subscripts specify the nature of each term.The sub-term of h ij are given by: where Onece obtained numerical solutions for the time-evolution of q and p of a spinning compact binary system through numerical methods, we can construct the corresponding gravitational waveforms.In the subsequent chapter, we will devise an optimized correction map in extended phase-space algorithm for this system.

Dissipative correction map in extended phase-space method
Since the Hamiltonian H can't 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 extended phase-space method.[Pihajoki:2015] introduced a new pair of conjugate and canonical variables ( r, p) from the original variables (r, p).This doubles the phase-space variables,(r, p)→(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). (43) When it comes to the Hamiltonian 13, the formulation 43 should be rewrited as: H(r, θ j , r, θ j ; p, ξ j , p, ξ j ) = H1(r, θ j , p, ξ j ) + H2( r, θ j , p, ξ j ).( 44) Observing that both (r, θ j , p, ξ j ) and ( r, θ j , p, ξ j ) constitute two pairs of conjugate canonical variables, it is immediately apparent that the newly formed Hamiltonian H comprises two distinct integrable components.
Given this property, it is reasonable to anticipate that the second-order leapfrog algorithm would be well suited to numerically whole the Hamiltonian H.The implementation of such splitting methods proceeds as follows: When H1 and H2 represent the operators that facilitate the analytical solution of the individual Hamiltonians H1 and H2, respectively, and h denotes the chosen time-step, a standard leapfrog algorithm can be expressed as: It is noteworthy that, given identical initial conditions, the solution pairs (r, p) and ( r, p) are anticipated to coincide at each time step.Nevertheless, their trajectories rapidly deviate in subsequent time steps due to the intricate interdependence between the solution (r, p) derived from H1 and the solution ( r, p) associated with H2.This phenomenon manifests itself as a compensatory relationship, where any increase in the value of H1 over half of the entire Hamiltonian H, is accompanied by a commensurate decrease in H2, and vice versa.Motivated by the inherent symmetry between H1 and H2, [Luo:2017ApJ] devise a midpoint map and ensure the equality between H1 and H2, 2 , 0, 0, 0, 0, 0, 0 1 2 , 1 2 , 0, 0, 0, 0, 0, 0 0, 0, 1 2 , 1 2 , 0, 0, 0, 0 0, 0, 1 2 , 1 2 , 0, 0, 0, 0 0, 0, 0, 0, 1 2 , 1 2 , 0, 0 0, 0, 0, 0, 1 2 , 1 2 , 0, 0 0, 0, 0, 0, 0, 0, 1 2 , 1 2 0, 0, 0, 0, 0, 0, The purpose of map M1 is to take the midpoint values between original variables and their corresponding duplicate variables and reassign these midpoints to both the original and duplicate variables, for example, r = r = (r+ r)/2.This operation effectively aligns initially unequal pairs of original and replica variables, ensuring they become identical and thus preventing further divergence during subsequent evolutionary processes.Subsequently, the leapfrog algorithm combined with the midpoint map can be formulated as: From the nth time step to the subsequent (n + 1)th step, the corresponding numerical solutions can be represented as When gravitational dissipation is taken into account, the constants of motion associated with the Hamiltonian are no longer conserved.Consequently, in evaluating the absolute energy error for Algorithm C2, we cannot employ the conventional approach E − E0, where E is the instantaneous energy and E0 is the initial energy.Instead, we must resort to Equation 32, which accounts for the energy dissipation due to gravitational radiation.To proceed, let us first reexpress Equation 32 in a suitable form for our purposes.

Numerical simulations
Our primary focus lies in assessing the performance of the algorithms outlined in Section 2 when applied to controlling numerical errors in post-Newtonian (PN) systems of spinning compact binaries, as modeled by the Hamiltonian formulation given in Equation 1.This ten-dimensional canonical spin Hamiltonian possesses four integrals of motion: the total energy and the three components of the total angular momentum vector.However, the lack of a fifth integral renders the system nonintegrable, potentially giving rise to chaotic behavior in certain spin Hamiltonians, as demonstrated in [Zhong:2010PRD], [Mei 2013EPJC, Mei:2013uqa], [Luo:2020] and [Luo 2022].
For comparative purposes, an 4-stage implicit Gaussian algorithm with eighth order accuracy, denoted Gauss4, will also be employed to solve Hamiltonian 1, serving as a reference solution.
In Figure 1, we plot the absolute energy errors △EC2 for each algorithm with fixed time step h = 0.5, revealing that C2 exhibits the poorest error behavior, with CM 2 demonstrating marginally better Figure 2: Phase space distance D ps between Gauss4 and other algorithms as functions of time steps.The distances in the ascending order are CM 3 (green dot), CM 2 (red dash), C 2 (black).The D ps of CM 2 starts out the smallest distance but keeps growing and eventually becomes the second closest, while CM 3 rises to the smallest distance.C 2 is the longest distance and ranking last.
performance.Both algorithms display a slight energy offset, which is not uncommon in dissipative systems.On the other hand, CM 3 exhibits the highest precision, approaching the limit of double-precision arithmetic on our computing platform, and has excellent stability.Given CM 3's inherent characteristic of correction map M 3 to minimize △EC2, this superior accuracy is anticipated; however, it renders the comparison of energy errors alone insufficient for impartially evaluating the overall performance of algorithms.
To address this, in Figure 2, we present the phase space distances between the numerical solutions produced by each algorithm and the reference solutions obtained using 4-stage implicit Gaussian method.Here Dps = [(r, θ j ; p, ξ j )gauss − (r, θ j ; p, ξ j )] 2 , the solutions of Gauss4 are denoted as (r, θ j ; p, ξ j )gauss.Initially, CM 2 displays the shortest phase space distance, but it rapidly deteriorates and becomes inferior to CM 3, which subsequently maintains the closest proximity to the true solution.Throughout this period, C2 consistently remains the farthest from the reference.This plot, while ultimately reaching similar conclusions as Figure 1 regarding algorithmic performance, reveals a little differences in the temporal evolution, particularly highlighting CM 2 as initially outperforming the others.
Beyond phase space distance, we further assess algorithmic performance by subtracting the dissipation energy calculated by each algorithm from that computed using the highly accurate 4-stage implicit Gaussian Since the latter provides the most reliable estimate of dissipation energy, this approach offers an objective basis for comparison.In accordance with this strategy, we generate Figure 3.It can be seen that figure 3 demonstrates striking similarities in the evolutionary trends with those observed in Figure 2, with CM 2 initially outperforming its counterparts, followed by a decline in performance and eventual supersession by CM 3, while C2 consistently underperforms throughout.n=0 △H(n) Gauss calculated by the higher-order algorithm Gauss4.C 2 is a solid black line, CM 2 is represented by a red dash line, and CM 3 is a green dotted line.It can be seen that the evolution law of each algorithm is highly similar to the phase space distance in Figure 2.
Figure 4 presents gravitational waveform plots generated by various algorithms with setting the direction p = (1, 0, 0) and the orientation of the observer N = (0, sin(π/4), cos(π/4)).Specifically, Figure 4a depicts the waveform for hx, where it is immediately apparent that the waveforms produced by C2, CM 2, CM 3, and Gauss4 are virtually indistinguishable to the naked eye due to their near-complete overlap.However, upon closer inspection through local magnification, discernible differences become evident.Figure 4b serves this purpose, providing an amplified view of the hx waveform evolution.Here, it is clear that CM 3 yields the waveform closest in resemblance to that produced by the high-order Gauss4 algorithm, followed by CM 2, with C2 exhibiting the greatest deviation from Gauss4.
Figure 4c then displays the hy component of the gravitational waves, generated by each algorithm.Similarly, at a cursory glance, the waveforms appear largely similar.To reveal the nuances, Figure 4d offers a magnified look at the local details of the hy waveforms, again revealing that CM 3 maintains the closest alignment with the highly accurate Gauss4 waveform, while the other algorithms follow in decreasing proximity.
Lastly, Figure 5 presents a graphical rendering of the binary star's orbital trajectory in configuration space as calculated by Gauss4.This visual representation serves to provide readers with a more intuitive understanding of the underlying physical dynamics, complementing the waveform analysis and offering a comprehensive perspective on the system's behavior as modeled by the highest precision algorithm under consideration.
To delve deeper into the comparative performance of these algorithms, we conduct numerical simulations for a distinct orbit, designated Orbit 2, characterized by the following initial conditions:(β; r, p) = (1; 23, 0, 0, 0, 0.24, 0),χ1 = χ2 = 1, Ŝ1 = (ρ1 cos π 4 , ρ1 sin π 4 , −0.983734), Ŝ2 = (ρ2 cos π 4 , ρ2 sin π 4 , −0.983734), ρ1 = ρ2 = 1 − (−0.983734) 2 .and setting fix time-step h = 0.6 While there is a degree of consistency in the overall performance of each algorithm between Orbit 1 and Orbit 2, subtle differences emerge, except the patterns of energy errors.From figure 6 we see that the energy error △EC2 evolution law drawn in orbit 2 and orbit 1 has the same conclusion, and CM 3 is still the most accurate and stable long-term evolution, CM 2 ranks second with a gap of several orders of magnitude, and C2 ranks last slightly below CM 2. The change of the phase space distance diagram in Figure 7 compared for orbit 1 is that CM 3 maintains the minimum distance from the exact solution.CM 2 does not have an advantage like orbit 1 at the beginning and then lower than CM 3, but slightly worse than CM 3 in the long run.Of course, both CM 2 and CM 3 outperform C2 in approximating the exact phase space behavior, the comparison of their performances in orbit 2 reveals that CM 3 sustains its role as the superior integrator.
Consistent with the observations from the phase space distance analysis, the energy error diagram in Figure 8, which employs the Gauss4-calculated dissipative energy as the reference truth value, further confirms the performance hierarchy of the investigated algorithms.Once again, CM 3 emerges as the most accurate, followed closely by CM 2, while C2 occupies the third position.The evolutionary trends displayed in this energy error plot strikingly resemble those encountered in Figure 7, highlighting the strong correlation between the phase space distance and energy error metrics in characterizing the efficacy of these numerical integration schemes, whether it's orbit 1 or orbit 2. The congruity between the phase space distance and energy error profiles underscores the fact that both measures are effectively capturing the same fundamental aspect of algorithmic performance: the ability to accurately track the true dynamics of the spinning compact binary system.This strong correlation implies that the conclusions drawn from analyzing any one indicator in isolation are likely to echo those of the other, enhancing the reliability of the overall assessment.(A) and (C) show the overall gravitational waveforms for the 'hx' and 'h+' polarizations, respectively, providing a global perspective on the radiation emitted by the spinning compact binary during orbit 2. Across these comprehensive views, all algorithms yield waveforms that are essentially indistinguishable from one another, indicating a high degree of agreement in their overall representation of the gravitational wave signal.
(B) and (D) investigate the local amplification details of the 'hx' and 'h+' polarizations, where subtle differences in the waveforms become apparent.In these magnified sections, CM 3 is revealed to exhibit the closest resemblance to the benchmark Gauss4 calculation, and slightly closer than CM 2, with C2 demonstrating the greatest deviation.This finding aligns with the conclusions drawn from the analysis of orbit 1.
Lastly, Figure 10 illustrates the orbital evolution of the binary system in configuration space as modeled by Gauss4 throughout orbit 2. This visualization offers a complementary perspective on the underlying dynamical processes driving the observed gravitational wave patterns, providing a more holistic understanding of the system's behavior.
Despite these distinctions, the general trends observed across the two orbits provide valuable insights into the strengths and weaknesses of the employed numerical methods, ultimately informing their suitability for accurately simulating a wide range of spinning compact binary configurations.

Summary
The present work has focused on the development and evaluation of the dissipated correction map method with the trapezoidal rule for numerical simulations of gravitational waves emitted by spinning compact binary systems.Our objective was to advance the frontiers of simulating inherently intricate and complex celestial phenomena, particularly concentrating on enhancing the precision and stability of extended phase-space algorithms.The proposed dissipated correction map builds upon the foundation of extended phase-space techniques, tackling key challenges encountered in modeling the complex dynamics of these systems using post-Newtonian approximations.Here, we highlight the substantial enhancements offered by the dissipated correction map over the midpoint map and their implications for future research in gravitational wave astronomy.
In the numerical simulations and performance assessment, we rigorously test the effectiveness of the CM 3, we carried out extensive numerical simulations using a ten-dimensional canonical spin Hamiltonian, known for its tendency to display chaotic dynamics due to the absence of a fifth integral.This choice of this Hamiltonian allowed us to appraise the correction map method's performance under conditions that are particularly taxing for numerical integrators.The simulations were performed for Orbit 1 and 2, and the outcomes were contrasted against those obtained using the midpoint map which was the already established algorithms.Our analyses revealed several compelling advantages of the CM 3 over the C2 and alternative approaches.Energy error evaluations demonstrated that the CM 3 consistently achieved the lowest energy errors, signifying superior conservation attributes and a heightened level of accuracy in tracing the system's energy evolution.Furthermore, the dissipated energy comparison showed that the CM 3 closely paralleled the results obtained from the high-precision Gaussian algorithm, further supporting its accuracy in representing the dynamics of the binary system.Temporal stability evaluations, quantified via phase space distance analysis, unmistakably showed that the CM 3 outperformed the CM 2 and C2 over time.This means that the numerical solution of CM 3 is closer to the Gauss4 algorithm as the truth value.
Visual scrutiny of the simulated gravitational waveforms provided compelling testament to the CM 3 superiority.Gravitational waveforms computed using the CM 3 closely mimicked those generated by the high-precision Gaussian algorithm, particularly in regions of local magnification, indicating that the CM 3 captures the fine-grained details of the wave signal with better performance.Therefore, compared with the previous C2 algorithm, CM 2 or CM 3, especially CM 3, is recommended to simulate the gravitational waveform of spinning compact binaries .

Figure 1 :
Figure1: Energy error of △EC2 calculated by extended phase-space method with all maps.the CM 3 algorithm (green dot) consistently displays the highest levels of accuracy and long-term stability among the tested methods.Conversely, the energy calculations produced by C 2 (black) exhibit the most pronounced bias.The CM 2 scheme (red dashes) offers a marginal improvement in accuracy compared to C 2 .

Figure 4 :
Figure 4: The gravitational waveform h x and h + in orbit 1. (A).The global evolution diagram of h x drawn by C 2 , CM 2, CM 3 and Gauss4, they are almost the same.(B).h x local magnification diagram, it can be seen that the result of CM 3 is closest to Gauss4, followed by CM 2, and finally C 2 .(C).The global evolution graph of h + , with almost no difference among all algorithms.(D).h + locally enlarged figure shows that the rankings of evolution closest to Gauss4 are CM 3, CM 2, and C 2 .

Figure 5 :
Figure 5: The projection of orbit 1 onto x − y − z space plot by Gauss4.

Figure 6 :
Figure6: The extended phase-space method, incorporating all maps, was utilized to compute the energy error of △EC2.In particular, the CM 3(green dot) algorithm consistently demonstrated the highest degrees of precision and long-term stability among the methods tested.On the other hand, the energy calculations generated by C 2 (black) exhibited the most significant deviation.The CM 2(red dash) scheme offered a slight enhancement in accuracy compared to C 2 .

Figure 7 :
Figure 7: Phase space distance D ps between Gauss4 and other algorithms as functions of time steps.The distances in the ascending order are CM 3 (green dot), CM 2 (red dask), C 2 (black).The D ps of CM 3 keeps the smallest distance, while CM 2 is farther away.C 2 is the longest distance and ranking last.

Figure 9 :
Figure 9: The gravitational waveform h x and h + in orbit 2. (A).The global evolution diagram of h x drawn by C 2 , CM 2, CM 3 and Gauss4, they are almost the same.(B).h x local magnification diagram, It can be seen that the waveforms of CM 3 and CM 2 are between that of Gauss4 and C 2 .And CM 3 is slightly closer to Gauss4 than CM 2. (C).The global evolution graph of h + , with almost no difference among all algorithms.(D).h + locally enlarged figure shows that the rankings of evolution closest to Gauss4 are similar to the ranking in h x .

Figure 9
Figure9presents four distinct diagrams, labeled (A), (B), (C), and (D), each offering specific insights into the gravitational wave signatures generated by the studied algorithms with p = (1, 0, 0) and N = (0, sin(π/4), cos(π/4)):(A) and (C) show the overall gravitational waveforms for the 'hx' and 'h+' polarizations, respectively, providing a global perspective on the radiation emitted by the spinning compact binary during orbit 2. Across these comprehensive views, all algorithms yield waveforms that are essentially indistinguishable from one another, indicating a high degree of agreement in their overall representation of the gravitational wave signal.(B)and (D) investigate the local amplification details of the 'hx' and 'h+' polarizations, where subtle differences in the waveforms become apparent.In these magnified sections, CM 3 is revealed to exhibit the closest resemblance to the benchmark Gauss4 calculation, and slightly closer than CM 2, with C2 demonstrating the greatest deviation.This finding aligns with the conclusions drawn from the analysis of orbit 1.Lastly, Figure10illustrates the orbital evolution of the binary system in configuration space as modeled by Gauss4 throughout orbit 2. This visualization offers a complementary perspective on the underlying dynamical processes driving the observed gravitational wave patterns, providing a more holistic understanding of the system's behavior.Despite these distinctions, the general trends observed across the two orbits provide valuable insights into the strengths and weaknesses of the employed numerical methods, ultimately informing their suitability for accurately simulating a wide range of spinning compact binary configurations.

Figure 10 :
Figure 10: The projection of orbit 2 onto x − y − z space plot by Gauss4.