The orbital evolution and gravitational waves of OJ 287 in the 4th post-Newtonian order

The supermassive black hole (SMBH) binary system OJ 287 is one of the few most frequently observed astrophysical objects, which presents us with an excellent test subject for General Relativity. Owning to the fact that all of its initial parameters are given to high accuracy by [1, 2], we can analyze its entire orbital evolution. OJ 287 is located at redshift z = 0.306, and the masses of the primary and secondary SMBHs are m₁ = 1.84 × 10¹⁰ M_⊙ and m₂ = 1.46 × 10⁸ M_⊙ which makes the components of this system one of the most massive objects in the observable universe. For centuries, astronomers have observed and examined its quasi-periodic light variations, yet it has only recently been recognized as a blazar. The optical light curve of OJ 287 of the last $\sim 120$ years has been showing a double periodicity of 60 years and 12 years [3, 4]. In 1996 a detailed model was introduced by Letho and Valtonen [5, 6]. In the model, the proposed 12 years orbital period is related to the major outbursts of the primary SMBH. The mechanism that forms the pair of outbursts occurring, 1 to 2 years separated, was assumed to be caused by impacts on the accretion disk of the primary SMBH by the secondary. During the transition of the secondary SMBH, matter from the accretion disk is taken away, and bremsstrahlung bursts are formed as the hot plasma becomes optically thin. This outburst structure allows us to reproduce the orbits of the binary with high accuracy. A more detailed analysis of the outburst structure can be found in [1, 2]. From the luminosity of the source, we can determine both the masses and the relative velocity. It turned out that this system has a very high velocity of v/c ≃ (0.06–0.26), and therefore it is an ideal candidate to test strong gravitational effects.

From the timing of outbursts, Valtonen and his associates [1, 2, 7] could accurately estimate the spin of the primary SMBH and show its significance in the orbital evolution. Pihajoki [8] showed that the secondary SMBH could also produce bremsstrahlung outbursts, which implies that the secondary SMBH can also spin up. We found that in the analysis of gravitational waves (GWs), the contribution of the secondary spin is not significant enough to be considered in this article. Here we used the highest spin–spin (SS) [9, 10], spin–orbit (SO) [11, 12] terms and post-Newtonian (PN) [13–18] up to the newly released 4th PN order [19].

OJ 287 is composed of a pair of close SMBHs, and such a binary system emits intense gravitational radiation. The large amplitude of its radiation makes OJ 287 an excellent target for low-frequency gravitational-wave detectors [20]. Earlier, Zhang and his colleagues analyzed the gravitational waveform of OJ 287. They made their analysis up to 3.5PN order and their results for three periods are in harmony with the findings of Valtonen. Furthermore, they investigated the detectability of GWs emitted by the system [21]. From the frequency of the binary system, it appears that we can detect them with pulsar timing array (PTA) [22] experiments. In this paper, we extend their results by going until the start of the merger phase, with the use of fourth post-Newtonian (4PN) and spin terms. We will examine the orbital evolution and waveforms after the last stable orbit (LSO) and use OJ 287 as an experimental case for testing 4PN results.

The inspiral of compact astrophysical objects can be described accurately in the weak-field limit, where the gravitational potential and orbital velocity are regarded as small parameters. In the PN regime [18], the motion of the binary system can be characterized as perturbed Keplerian orbits. The equation of motion and the emitted radiation of these systems are analyzed in detail by [9, 10, 15–17, 23–27]. As a result, thoroughly tested and ready to use formulas are available up to 4PN order.

To study the orbital motion, the evolution of spins and angular momenta of the binary system, we used the computational tool called CBwaves [13]. The equation of motion was integrated numerically by a fourth-order Runge–Kutta method. We determined the radiation field by a simultaneous evaluation of the analytic waveform involving all the contributions up to 2PN order. As an output of the code, both time and frequency domain waveforms are available. The PN contributions to the acceleration and the radiation field are listed in the appendices of [13].

Improving previous results on spin effects, we added 3.5PN SO contribution to the equations of motion, see e.g. [9, 14], and spin precession and the total-energy flux emitted in GWs, as it can be found in [11, 12]. Furthermore, in 2018, Blanchet et al [19] released the 4PN term for acceleration in the center-of-mass frame. The differential equation of the motion solved by CBwaves can be written as the sum of the following terms

$$\begin{align}\hfill \mathbf{a}& ={\mathbf{a}}_{\text{PN}}+{\mathbf{a}}_{2\text{PN}}+{\mathbf{a}}_{\text{SO}}+{\mathbf{a}}_{\text{SS}}+{\mathbf{a}}_{\text{RR}}+{\mathbf{a}}_{\text{PNSO}}+{\mathbf{a}}_{3\text{PN}}\hfill \\ \hfill & \quad +{\mathbf{a}}_{\text{RR}1\text{PN}}+{\mathbf{a}}_{3.5\text{PNSO}}+{\mathbf{a}}_{\text{RRSO}}+{\mathbf{a}}_{\text{RRSS}}+{\mathbf{a}}_{4\text{PN}}.\hfill \end{align} \tag{ 1 }$$

The exact form of the different contributions can be found up to 3PN in the appendix of [13]. For the sake of completeness, the newly introduced 3.5PN SO terms, a_3.5PNSO, spin precession, and 4PN corrections are presented here in the appendixes. We used here the radiation-reaction term, derived from the Burke–Thorne radiation-reaction potential. Detailed description of the gauge parameters can be found in [17, 28] and in this paper we collected them in table 1.

We have obtained the next-to-leading order correction to the SO interaction in the center-of-mass frame up to 3.5PN order from [11], and we used is equation (3.8), which can be expressed as

$$\begin{equation}{\mathbf{a}}_{\mathrm{3.5}\mathrm{P}\mathrm{N}\mathrm{S}\mathrm{O}}=\frac{Gm}{{r}^{3}}\left[{}^{(3.5)}\mathbf{b}_{1}+{}^{(3.5)}\mathbf{b}_{2}\frac{Gm}{r}+{}^{(3.5)}\mathbf{b}_{3}\frac{{G}^{2}{m}^{2}}{{r}^{2}}\right],\end{equation} \tag{ 2 }$$

where the exact form of the vector coefficients (^(3.5) b₁,^(3.5) b₂,^(3.5) b₃) can be found in [11]. The spin evolution have been calculated in CBwaves by

$$\begin{equation}{\dot{\mathbf{S}}}_{1}=\frac{G}{{c}^{2}{r}^{3}}\left\{{{\Omega}}_{1}{\mathbf{L}}_{\text{N}}-{\mathbf{S}}_{2}-3({\hat{\mathbf{n}}\cdot S}_{2})\hat{\mathbf{n}}+\frac{{G}^{2}\mu m}{{c}^{5}{r}^{2}}\left[\frac{2}{3}(\mathbf{v}\cdot {\mathbf{S}}_{2})+30\dot{r}(\hat{\mathbf{n}}\cdot {\mathbf{S}}_{2})\right]\mathbf{n}\right\}\times {\mathbf{S}}_{1}.\end{equation} \tag{ 3 }$$

The above equation contains the standard SO and SS [16], the 2PN SO [10], the 3PN SO [11] and the 3.5PN radiative SS [26] contributions. Let us note that in CBwaves, we presume the use of covariant spin supplementary condition, ${S}_{A}^{\mu \nu }=0$, where ${u}_{A}^{\mu }$ is the four-velocity of the center-of-mass world line of body A, with A = 1, 2. Furthermore, we can find the spin precession correction Ω₁ to the SO contribution terms in appendix B. Here we restrict ourselves to the following form

$$\begin{equation}{{\Omega}}_{1}=\left[\frac{1}{{c}^{2}}({\alpha }_{\text{PN}}+{\alpha }_{\text{PN,Kidder}})+\frac{1}{{c}^{4}}({\alpha }_{2\text{PN}}+{\alpha }_{2\text{PN,Blanchet}})+\frac{1}{{c}^{6}}{\alpha }_{3\text{PN}}+\mathcal{O}\left(\frac{1}{{c}^{7}}\right)\right].\end{equation} \tag{ 4 }$$

Simultaneously with the orbital evolution, the radiation field h_ij is determined by the analytic waveform contributions up to 2PN order in harmonic coordinates, valid for general eccentric and spinning sources. This can be decomposed as

$$\begin{align}\hfill {h}_{ij}& =\frac{2G\mu }{{c}^{4}D}\left[{Q}_{ij}+{P}^{0.5}{Q}_{ij}+P{Q}_{ij}+{P}^{1.5}{Q}_{ij}+{P}^{2}{Q}_{ij}+P{Q}_{ij}^{\text{SO}}\right.\hfill \\ \hfill & \quad \left.+{P}^{1.5}{Q}_{ij}^{\text{SO}}+{P}^{2}{Q}_{ij}^{\text{SO}}+P{Q}_{ij}^{SS}+{P}^{1.5}{Q}_{ij}^{\text{tail}}\right],\hfill \end{align} \tag{ 5 }$$

where D is the distance to the source, and μ = m₁ m₂/(m₁ + m₂) is the reduced mass of the binary. Q_ij is the quadrupole (or Newtonian) term, P^0.5 Q_ij , PQ_ij , P^1.5 Q_ij , P² Q_ij [29] are higher-order relativistic corrections, ${P}^{1.5}{Q}_{ij}^{\text{tail}}$ is the 1.5PN tail term, $P{Q}_{ij}^{\text{SO}}$, ${P}^{1.5}{Q}_{ij}^{\text{SO}}$, ${P}^{2}{Q}_{ij}^{\text{SO}}$, and $P{Q}_{ij}^{\text{SS}}$ [16] are the SO and SS corrections. The full expressions of these terms are found, in detail, in the appendix of [13]. Using a fourth-order Runge–Kutta method, we could model the evolution and the emitted GWs of general binary systems with high accuracy.

Looking at a plane wave traveling in direction $\hat{\mathbf{N}}$, which is the unit spatial vector pointing from the center of mass of the source to the observer, the transverse-traceless (TT) part of the radiation field is provided by

$$\begin{equation}{h}_{ij}^{\text{TT}}={{\Lambda}}_{ij,kl}{h}_{kl},\end{equation} \tag{ 6 }$$

where

$$\begin{equation}{{\Lambda}}_{ij,kl}(\hat{\mathbf{N}})={\mathbf{P}}_{ik}{\mathbf{P}}_{jl}-\frac{1}{2}{\mathbf{P}}_{ij}{\mathbf{P}}_{kl}\end{equation} \tag{ 7 }$$

and

$$\begin{equation}{\mathbf{P}}_{ij}(\hat{\mathbf{N}})={\delta }_{ij}-{N}_{i}{N}_{j}.\end{equation} \tag{ 8 }$$

The radiation field can be defined by choosing an orthonormal triad as in [13], where

$$\begin{equation}\hat{\mathbf{N}}=(\mathrm{sin}\,\iota \,\mathrm{cos}\,\phi ,\,\mathrm{sin}\,\iota \,\mathrm{sin}\,\phi ,\,\mathrm{cos}\,\iota ),\end{equation} \tag{ 9 }$$

$$\begin{equation}\hat{\mathbf{p}}=(\mathrm{cos}\,\iota \,\mathrm{cos}\,\phi ,\,\mathrm{cos}\,\iota \,\mathrm{sin}\,\phi ,-\mathrm{sin}\,\iota ),\end{equation} \tag{ 10 }$$

$$\begin{equation}\hat{\mathbf{q}}=(-\mathrm{sin}\,\phi ,\mathrm{cos}\,\phi ,0),\end{equation} \tag{ 11 }$$

and ι and ϕ are the polar angles determining the relative position of the radiation frame with respect to the source frame $(\hat{\mathbf{x}},\hat{\mathbf{y}},\hat{\mathbf{z}})$, see the first figure of [13]. The polarization states of the gravitation wave can be given, with respect to the orthonormal radiation frame, as

$$\begin{equation}{h}_{+}=\frac{1}{2}\left({\hat{p}}_{i}{\hat{p}}_{j}-{\hat{q}}_{i}{\hat{q}}_{j}\right){h}_{ij}^{\text{TT}},\quad {h}_{\times }=\left({\hat{p}}_{i}{\hat{q}}_{j}-{\hat{q}}_{i}{\hat{p}}_{j}\right){h}_{ij}^{\text{TT}}.\end{equation} \tag{ 12 }$$

The strain produced by the binary system at the detector can be given with the combination

$$\begin{equation}h(t)={F}_{+}{h}_{+}(t)+{F}_{\times }{h}_{\times }(t),\end{equation} \tag{ 13 }$$

where F₊ and F_× are the antenna pattern functions, which can be found in detailed form in [13].

In our numerical analysis, we used a program package named CBwaves developed by Csizmadia et al [13] and maintained by our group. This package provides a fast and accurate computational tool to examine the GWs emitted by spinning eccentric binary systems. We achieve this by integrating the center-of-mass equation of motion within the PN framework up to the 4th term, plus the spin precession equations up to the 3rd term developed by Bohé et al [11, 12].

The initial parameters for our analysis are chosen from the estimates of Voltanen [1, 2], and Zhang et al [20]. The values are listed in table 2.

Table 2. Initial parameters.

m⊙	(Meter)	1476.62504
m1	(M⊙)	1.84 × 1010
m2	(M⊙)	1.46 × 108
r0	(M)	$\sim 100$
${S}_{1}/{m}_{1}^{2}$		0.381
${S}_{2}/{m}_{2}^{2}$		0.996
Δϕ	(deg)	39.1
e0		0.658

Where m₁ is the mass of the primary SMBH and m₂ is the mass of the secondary SMBH, r₀ is the initial separation which we set to be approximately 100(m₁ + m₂ = M), Δϕ is the precession rate of the orbit per period, e₀ is the initial eccentricity of the secondary mass' orbit. We have chosen a different initial separation from the value, r₀ ∼ 0.535 pc ≃ 60M, which can be found in [20] for the sake of simplicity.

In [8], we found that both SMBH owns spin, but in our analysis, we assumed that only the primary SMBH has a spin, and it points to the positive $\hat{\mathbf{z}}$ direction, perpendicular to the orbital plane. In earlier computations, we have found that the contribution of the secondary SMBH's rotation is relatively small in the waveform and the orbital evolution of OJ 287. We had run simulations of the binary system in the inspiral phase and let the separation of SMBHs evolve beyond the last stable or into the so-called unstable orbit. We have estimated that the merger of the two SMBH would start slightly later than $\sim 17\,000\enspace \text{years}$. We stopped our simulation before the final unstable circular orbit ends, in other words, before the separation of the two bodies becomes smaller than r = 3M. With this requirement, we also prevented the event horizons from coming to contact with each other, which would be beyond the capabilities of the program package we used. The orbit of the secondary SMBH and the total and radiated energies have been presented in figure 1. As we can see, as long as the two SMBH are separated more than by 6M, the total energy of the system remains seemingly constant, which means the orbit of the secondary SMBH is stable. But by looking at the radiated energy we see it almost stops growing as the orbit of the secondary SMBH converges in. It gives us some confidence that the calculations still give physically relevant results even beyond the LSO.

Zoom In Zoom Out Reset image size

In order to further prove our point, we have analyzed the maximal orbital separation (apoapsis points of the elliptic orbit) and PN parameter ∼ (v/c)² ∼ GM/(rc²) of the 3PN- and 4PN-accurate simulations, as shown in figure 2.

Zoom In Zoom Out Reset image size

As we can see in figure 2, the separation in both 3PN and 4PN approximations, marked by solid green and dashed orange lines, changes similarly. On the other hand, the PN parameters in the 3PN and 4PN approximations are different beyond the LSO. At the 3PN order equation, the PN parameter ∼ (v/c)² has an inflection point, and it rapidly starts to decrease shortly after leaving the LSO. On the other hand, at the 4PN order equation, the PN parameter grows steadily even after the LSO. It shows that the 3PN order equation of motion changes its sign shortly after the LSO, which means the binary system starts to decouple. Consequently, the 4PN contribution provides us with a more accurate description of the orbital evolution beyond the LSO than the computations with the 3PN-accurate equation. According to reference [13], PN formalism is valid if ∼ 0.08–0.1, however in section 9.6 of reference [18], Blanchet argued that the limit of the formalism is when the orbital velocities are of the order of 50% of the speed of light. In our computation, the orbital velocity was $v=\sqrt{{\epsilon}}\sim 0.35c$, in other words, $\sim 35\%$ of the speed of light, which clearly stays well below the limit set by Blanchet.

Following the computations presented above, now we may analyze the unstable orbit in the 4PN regime. We have shown its orbital evolution in 3PN and 4PN accuracy in figure 3. The total and radiated energies of the system are shown in figure 4.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figures 3 and 4 show further proof of our assumption, as we find more differences in the behavior orbit at 3PN and 4PN orders. We notice that the orbit at the 3PN order makes a similar number of cycles before and after the LSO. The orbit at the 4PN order reaches the 3M limit within a few revolutions beyond the LSO. Furthermore, if we look at figure 4, we see that the total and radiated energies of the binary start to rapidly decrease and grow much earlier than in the 4PN term. Following the arguments presented above, we can declare that the analysis of the unstable orbit of eccentric, extreme mass ratio, supermassive binary systems require the 4PN term.

We have shown in figures 5 and 6 the strain of the GWs emitted by the binary, first for the entire orbital evolution, then for after the LSO. In order to obtain the proper amplitudes, one has to know the distance of the source from the observer. We used the distance D found in [20] as follows

$$\begin{equation}{D}_{\text{L}}(z)=\frac{c}{{H}_{0}}(1+z)\underset{0{}}{\overset{z}{\int }}\frac{\mathrm{d}{z}^{\prime }}{\sqrt{{{\Omega}}_{m}(1+{z}^{\prime })+{{\Omega}}_{{\Lambda}}}},\end{equation} \tag{ 14 }$$

where Ω_m = 0.685, Ω_Λ = 1 − Ω_m are cosmological constants in a spatially flat ΛCDM universe, and H₀ = 67.3  (km s⁻¹) Mpc⁻¹ is the Hubble constant. Then the distance of OJ 287 is D = D_L(0.306) = 1.646 89 × 10⁹ pc, from here we can calculate the prefactor amplitude of equation (5), which is $\frac{2G\mu }{{c}^{4}D}=8.6539\times 1{0}^{-15}$.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As it has been demonstrated in figure 7, OJ 287 is a perfect candidate for detection techniques that target astrophysical objects like neutron stars. These methods are known as PTAs [22], and the currently ongoing experiments are the European pulsar timing array (EPTA) [30], the International pulsar timing array (IPTA) [31] and square kilometer array (SKA) [32]. Nevertheless, OJ 287's orbital frequency is insufficiently low to be detected by our current ground-based gravitational-wave interferometers. The detector sensitivity curve of these detectors and the characteristic strain of OJ 287 are presented in figure 7. The characteristic strain is attained by Fourier transformation of the GW amplitude,

$$\begin{equation}\tilde{h}=\mathcal{F}\left\{h(t)\right\}(f)=\underset{-\infty {}}{\overset{\infty }{\int }}h(t){\text{e}}^{-2\pi \text{i}ft}\,\mathrm{d}t,\end{equation} \tag{ 15 }$$

as it is found in reference [33]. Then we calculated the characteristic strain for the first orbit of our simulation and the unstable orbit in 3PN and 4PN order. As we earlier have seen, the contribution of the 4PN term to the equation of motion is initially small (Δa_3PN ∼ 10⁻³² and Δa_4PN ∼ 10⁻³⁵) compared to 3PN order.

Zoom In Zoom Out Reset image size

In this paper, we have studied the orbital evolution of the SMBH binary system OJ 287. We applied the PN [18] approximation up to the fourth-order correction [19] term and used the next-to-next-leading order of spin effects [9–12]. We have examined the evolution of the system beyond the LSO. By explicit computations of CBwaves, we have found that the 4PN order is essential in the evolution of orbits beyond the LSO since it makes a more numerically accurate evolution than the 3PN order. We have proved that by analyzing the total and radiated energies and the PN parameter of the system in detail, see figures 2 and 4. We have shown the waveform for the entire inspiral phase and presented the difference between the waveforms for 4PN and 3PN orders, see figure 6.

As it has been demonstrated in figure 7, OJ 287 is a perfect candidate for detection techniques that target astrophysical objects like neutron stars. These methods are known as PTAs [22] and currently, ongoing experiments are the EPTA [30], the IPTA [31] and SKA [32]. Nevertheless, OJ 287's orbital frequency is insufficiently low to be detected by our current ground-based gravitational-wave interferometers.

This work was supported by COST Action PHAROS (CA16214), the Hungarian NKFIH Grant Nos. 124366 and 124508. Furthermore, we would like to acknowledge the support of Wigner Scientific Computational Laboratory.

The data that support the findings of this study are available upon reasonable request from the authors.

We can find the exact calculations of the 4th PN order in [19], therefore, for the sake of simplicity, we can write the 4th order acceleration as

$$\begin{equation}{\mathbf{a}}_{4\text{PN}}=-\frac{Gm}{{r}^{2}}(A\mathbf{n}+B\mathbf{v}),\end{equation} \tag{ A.1 }$$

where we get the A coefficient as the sum of the following terms

$$\begin{align}\hfill {A}^{(0)}& =\left(\frac{315}{128}\eta -\frac{2205}{128}{\eta }^{2}+\frac{2205}{64}{\eta }^{3}-\frac{2205}{128}{\eta }^{4}\right){\dot{r}}^{8}+\left(-\frac{175}{16}\eta +\frac{595}{8}{\eta }^{2}-\frac{2415}{16}{\eta }^{3}\right.\hfill \\ \hfill & \quad +\left.\frac{735}{8}{\eta }^{4}\right){\dot{r}}^{6}{v}^{2}+\left(\frac{135}{8}\eta -\frac{1875}{16}{\eta }^{2}+\frac{4035}{16}{\eta }^{3}-\frac{1335}{8}{\eta }^{4}\right){\dot{r}}^{4}{v}^{4}+\left(-\frac{21}{2}\eta \right.\hfill \\ \hfill & \quad +\frac{1191}{16}{\eta }^{2}-\left.\frac{327}{2}{\eta }^{3}+99{\eta }^{4}\right){\dot{r}}^{2}{v}^{6}+\left(\frac{21}{8}\eta -\frac{175}{8}{\eta }^{2}+61{\eta }^{3}-54{\eta }^{4}\right){v}^{8}\hfill \end{align} \tag{ A.2 }$$

$$\begin{align}\hfill {A}^{(1)}& =\frac{m}{r}\left[\left(\frac{2973}{40}\eta +407{\eta }^{2}+\frac{181}{2}{\eta }^{3}-86{\eta }^{4}\right){\dot{r}}^{6}+\left(\frac{1497}{32}\eta -\frac{1627}{2}{\eta }^{2}\right.\right.\hfill \\ \hfill & \quad \left.\left.-81{\eta }^{3}+228{\eta }^{4}\right){\dot{r}}^{4}{v}^{2}+\left(-\frac{2583}{16}\eta +\frac{1009}{2}{\eta }^{2}+47{\eta }^{3}-104{\eta }^{4}\right){\dot{r}}^{2}{v}^{4}\right.\hfill \\ \hfill & \quad +\left.\left(\frac{1067}{32}\eta -58{\eta }^{2}-44{\eta }^{3}+58{{\epsilon}}^{4}\right){v}^{6}\right]\hfill \end{align} \tag{ A.3 }$$

$$\begin{align}\hfill {A}^{(2)}& =\frac{{m}^{2}}{{r}^{2}}\left\{\left[\left(\frac{2094\,751}{960}+\frac{45\,255}{1024}{\pi }^{2}\right)\eta +\frac{326\,101}{96}{\eta }^{2}-\frac{4305}{128}{\eta }^{3}-\frac{1959}{32}{\eta }^{4}\right]{\dot{r}}^{4}\right.\hfill \\ \hfill & \quad +\left[-\left(\frac{1636\,681}{1120}+\frac{12\,585}{512}{\pi }^{2}\right)\eta -\left(\frac{255\,461}{112}+\frac{3075}{128}{\pi }^{2}\right){\eta }^{2}-\frac{309}{4}{\eta }^{3}\right.\hfill \\ \hfill & \quad +63{\eta }^{4}\left.\right]{\dot{r}}^{2}{v}^{2}+\left[\left(\frac{1096\,941}{11\,200}+\frac{1155}{1024}{\pi }^{2}\right)\eta +\left(\frac{7263}{70}-\frac{123}{64}{\pi }^{2}\right){\eta }^{2}\right.\hfill \\ \hfill & \quad \left.+\frac{145}{2}{\eta }^{3}-16{\eta }^{4}\right]{v}^{4}\hfill \end{align} \tag{ A.4 }$$

$$\begin{align}\hfill {A}^{(3)}& =\frac{{m}^{3}}{{r}^{3}}\left\{\left[-2+\left(\frac{1297\,943}{8400}-\frac{2969}{16}{\pi }^{2}\right)\eta +\left(\frac{1255\,151}{840}+\frac{7095}{32}{\pi }^{2}\right){\eta }^{2}-17{\eta }^{3}\right.\right.\hfill \\ \hfill & \quad -\left.\left.24{\eta }^{4}\right]{\dot{r}}^{2}\left[\left(\frac{1237\,279}{252\,000}+\frac{3835}{96}{\pi }^{2}\right)\eta -\left(\frac{693\,947}{2520}+\frac{229}{8}{\pi }^{2}\right){\eta }^{2}+\frac{19}{2}{\eta }^{3}\right]{v}^{2}\right\}\hfill \end{align} \tag{ A.5 }$$

$$\begin{equation}{A}^{(4)}=\frac{{m}^{4}}{{r}^{4}}\left[25+\left(\frac{6625\,537}{12\,600}-\frac{4543}{93}{\pi }^{2}\right)\eta +\left(\frac{477\,763}{420}+\frac{3}{4}{\pi }^{2}\right){\eta }^{2}\right]\end{equation} \tag{ A.6 }$$

and we get B coefficient similarly to A as:

$$\begin{align}\hfill {B}^{(0)}& =\left(\frac{105}{16}\eta -\frac{245}{8}{\eta }^{2}+\frac{385}{16}{\eta }^{3}+\frac{35}{8}{\eta }^{4}\right){\dot{r}}^{7}+\left(-\frac{165}{8}\eta +\frac{1665}{16}{\eta }^{2}\right.\hfill \\ \hfill & \quad \left.-\frac{1725}{16}{\eta }^{3}-\frac{105}{4}{\eta }^{4}\right){\dot{r}}^{5}{v}^{2}+\left(\frac{45}{2}\eta -\frac{1869}{16}{\eta }^{2}+129{\eta }^{3}+54{\eta }^{4}\right){\dot{r}}^{3}{v}^{4}\hfill \\ \hfill & \quad +\left(-\frac{157}{16}\eta +54{\eta }^{2}-69{\eta }^{3}-24{\eta }^{4}\right)\dot{r}{v}^{6}\hfill \end{align} \tag{ A.7 }$$

$$\begin{align}\hfill {B}^{(1)}& =\frac{m}{r}\left[\left(-\frac{54\,319}{160}\eta -\frac{901}{8}{\eta }^{2}+60{\eta }^{3}30{\eta }^{4}\right){\dot{r}}^{5}+\left(\frac{25\,943}{48}\eta +\frac{1199}{12}{\eta }^{2}\right.\right.\hfill \\ \hfill & \quad \left.\left.-\frac{349}{2}{\eta }^{3}-98{\eta }^{4}\right){\dot{r}}^{3}{v}^{2}+\left(-\frac{5725}{32}\eta -\frac{389}{8}{\eta }^{2}+118{\eta }^{3}+44{\eta }^{4}\right)\dot{\left(\right.}\left.r\right){v}^{4}\right]\hfill \end{align} \tag{ A.8 }$$

$$\begin{align}\hfill {B}^{(2)}& =\frac{{m}^{2}}{{r}^{2}}\left\{\left[-\left(\frac{9130\,111}{3360}+\frac{4695}{256}{\pi }^{2}\right)\eta -\left(\frac{184\,613}{112}-\frac{1845}{64}{\pi }^{2}\right){\eta }^{2}\right.\right.\hfill \\ \hfill & \quad \left.+\frac{209}{2}{\eta }^{3}+74{\eta }^{4}\right]{\dot{r}}^{3}+\left[\left(\frac{8692\,601}{5600}+\frac{1455}{256}{\pi }^{2}\right)\eta \right.\hfill \\ \hfill & \quad \left.+\left(\frac{58\,557}{70}-\frac{123}{8}{\pi }^{2}\right){\eta }^{2}-70{\eta }^{3}-34{\eta }^{4}\right]\dot{r}{v}^{2}\hfill \end{align} \tag{ A.9 }$$

$$\begin{equation}{B}^{(3)}=\frac{{m}^{3}}{{r}^{3}}\left[2-\left(\frac{619\,267}{525}-\frac{791}{16}{\pi }^{2}\right)\eta -\left(\frac{28\,406}{45}+\frac{2201}{32}{\pi }^{2}\right){\eta }^{2}+66{\eta }^{3}16{\eta }^{4}\right]\dot{r}\end{equation} \tag{ A.10 }$$

where η = μ/(m₁ + m₂) the reduced mass, μ = m₁ m₂/(m₁ + m₂), m = m₁ + m₂ the total mass of the binary system and n = r/r the normal vector. Furthermore v² is the absolute value of the speed of the center-of-mass, which we derive from v² = v ⋅ v, where v = v₁ − v₂ is the relative velocity. At last $\dot{r}=\mathbf{n}\cdot \mathbf{v}$ is the scalar product of the relative velocity and the normal vector.

Now we can turn our attention to the 4th order corrections of the E_4PN = E/μ energy and $\mathcal{L}=\mathcal{L}/{\mathcal{L}}_{\text{N}}$ orbital angular momentum of the system. First let see the terms of the energy corrections as it can be found in the 3rd section of [19]

$$\begin{equation}{\mathcal{E}}^{(0)}=\left(\frac{63}{256}-\frac{1089}{256}\eta +\frac{7065}{256}{\eta }^{2}-\frac{10\,143}{128}{\eta }^{3}+\frac{21\,735}{256}{\eta }^{4}\right){v}^{10}\end{equation} \tag{ A.11 }$$

$$\begin{align}\hfill {\mathcal{E}}^{(1)}& =\frac{m}{r}\left[\left(-\frac{35}{128}\eta +\frac{245}{128}{\eta }^{2}-\frac{245}{64}{\eta }^{3}+\frac{245}{128}{\eta }^{4}\right){\dot{r}}^{8}+\left(\frac{25}{32}\eta -\frac{125}{16}{\eta }^{2}+\frac{185}{8}{\eta }^{3}\right.\right.\hfill \\ \hfill & \quad \left.\left.-\frac{595}{32}{\eta }^{4}\right){\dot{r}}^{6}{v}^{2}+\left(\frac{27}{64}\eta +\frac{243}{32}{\eta }^{2}-\frac{1683}{32}{\eta }^{3}+\frac{4851}{64}{\eta }^{4}\right){\dot{r}}^{4}{v}^{4}+\right.\hfill \\ \hfill & \quad \left.\times \left(-\frac{147}{32}\eta +\frac{369}{32}{\eta }^{2}+\frac{423}{8}{\eta }^{3}-\frac{4655}{32}{\eta }^{4}\right){\dot{r}}^{2}{v}^{6}\right.\hfill \\ \hfill & \quad \left.+\left(\frac{525}{128}-\frac{4011}{128}\eta +\frac{9507}{128}{\eta }^{2}-\frac{357}{64}{\eta }^{3}-\frac{15\,827}{128}{\eta }^{4}\right){v}^{8}\right]\hfill \end{align} \tag{ A.12 }$$

$$\begin{align}\hfill {\mathcal{E}}^{(2)}& =\frac{{m}^{2}}{{r}^{2}}\left[\left(-\frac{4771}{640}\eta -\frac{461}{8}{\eta }^{2}-\frac{17}{2}{\eta }^{3}+\frac{15}{2}{\eta }^{4}\right){\dot{r}}^{6}+\left(\frac{5347}{384}\eta +\frac{19\,465}{96}{\eta }^{2}\right.\right.\hfill \\ \hfill & \quad \left.-\frac{439}{8}{\eta }^{3}-\frac{135}{2}{\eta }^{4}\right){\dot{r}}^{4}{v}^{2}+\left(\frac{15}{16}-\frac{5893}{128}\eta -\frac{12\,995}{64}{\eta }^{2}+\frac{18\,511}{64}{\eta }^{3}\right.\hfill \\ \hfill & \quad \left.\left.+\frac{2846}{16}{\eta }^{4}\right){\dot{r}}^{2}{v}^{4}+\left(\frac{575}{32}-\frac{4489}{128}\eta +\frac{5129}{64}{\eta }^{2}-\frac{8289}{64}{\eta }^{3}+\frac{975}{16}{\eta }^{4}\right){v}^{6}\right]\hfill \end{align} \tag{ A.13 }$$

$$\begin{align}\hfill {\mathcal{E}}^{(3)}& =\frac{{m}^{3}}{{r}^{3}}\left\{\left[-\left(\frac{2599\,207}{6720}-\frac{6465}{1024}{\pi }^{2}\right)\eta -\left(\frac{103\,205}{224}-\frac{615}{128}{\pi }^{2}\right){\eta }^{2}+\frac{69}{32}{\eta }^{3}\right.\right.\hfill \\ \hfill & \quad \left.\left.+\frac{87}{4}{\eta }^{4}\right]{\dot{r}}^{4}+\left[\frac{21}{4}\left(\frac{1086\,923}{1680}+\frac{333}{512}{\pi }^{2}\right)\eta +\left(\frac{206\,013}{560}+\frac{123}{64}{\pi }^{2}\right){\eta }^{2}\right.\right.\hfill \\ \hfill & \quad \left.\left.-\frac{2437}{16}{\eta }^{3}-\frac{141}{2}{\eta }^{4}\right]{\dot{r}}^{2}{v}^{2}+\left[\frac{273}{16}-\left(\frac{22\,649\,399}{100\,800}-\frac{1071}{1024}{\pi }^{2}\right)\eta \right.\right.\hfill \\ \hfill & \quad \left.\left.+\left(\frac{521\,063}{10\,080}-\frac{205}{128}{\pi }^{2}\right){\eta }^{2}+\frac{2373}{32}{\eta }^{3}-\frac{45}{4}{\eta }^{4}\right]{v}^{4}\right\}\hfill \end{align} \tag{ A.14 }$$

$$\begin{align}\hfill {\mathcal{E}}^{(4)}& =\frac{{m}^{4}}{{r}^{4}}\left\{\left[\frac{9}{4}-\left(\frac{1622\,437}{12\,600}-\frac{2645}{96}{\pi }^{2}\right)\eta -\left(\frac{289\,351}{2520}+\frac{1367}{32}{\pi }^{2}\right){\eta }^{2}+\right.\right.\hfill \\ \hfill & \quad \left.\left.\frac{213}{8}{\eta }^{3}+\frac{15}{2}{\eta }^{4}\right]{\dot{r}}^{2}+\left[\frac{15}{16}+\left(\frac{1859\,363}{16\,800}-\frac{149}{32}{\pi }^{2}\right)\eta \right.\right.\hfill \\ \hfill & \quad \left.\left.+\left(\frac{22\,963}{5040}+\frac{311}{32}{\pi }^{2}\right){\eta }^{2}-\frac{29}{8}{\eta }^{3}+\frac{1}{2}{\eta }^{4}\right]\right\}\hfill \end{align} \tag{ A.15 }$$

$$\begin{equation}{\mathcal{E}}^{(5)}=\frac{{m}^{5}}{{r}^{5}}\left[-\frac{3}{8}-\left(\frac{1697\,177}{25\,200}+\frac{105}{32}{\pi }^{2}\right)\eta -\left(\frac{55\,111}{720}+11{\pi }^{2}\right){\eta }^{2}\right]\end{equation} \tag{ A.16 }$$

and the terms of the angular momentum, where ${\mathcal{L}}_{\text{N}}=\mu \mathbf{r}\times \mathbf{v}$ the Newtonian term

$$\begin{equation}{\mathcal{L}}^{(0)}=\left(\frac{35}{128}-\frac{605}{128}\eta +\frac{3925}{128}{\eta }^{2}-\frac{5635}{64}{\eta }^{3}+\frac{12\,075}{128}\right){v}^{8}\end{equation} \tag{ A.17 }$$

$$\begin{align}\hfill {\mathcal{L}}^{(1)}& =\frac{m}{r}\left[\left(-\frac{5}{8}\eta +\frac{15}{8}{\eta }^{2}+\frac{45}{16}{\eta }^{3}-\frac{85}{16}{\eta }^{4}\right){\dot{r}}^{6}+\left(3\eta -\frac{45}{4}{\eta }^{2}-\frac{135}{16}{\eta }^{3}\right.\right.\hfill \\ \hfill & \quad \left.\left.+\frac{693}{16}{\eta }^{4}\right){\dot{r}}^{4}{v}^{2}+\left(-\frac{53}{8}\eta +\frac{423}{16}{\eta }^{2}+\frac{299}{16}{\eta }^{3}-\frac{1995}{16}{\eta }^{4}\right){\dot{r}}^{2}{v}^{4}\right.\hfill \\ \hfill & \quad \left.+\left(\frac{75}{16}-\frac{151}{4}\eta +\frac{1553}{16}{\eta }^{2}-\frac{425}{16}{\eta }^{3}-\frac{2261}{16}{\eta }^{4}\right){v}^{6}\right]\hfill \end{align} \tag{ A.18 }$$

$$\begin{align}\hfill {\mathcal{L}}^{(2)}& =\frac{{m}^{2}}{{r}^{2}}\left[\left(\frac{14\,773}{320}\eta +\frac{3235}{48}{\eta }^{2}-\frac{155}{4}{\eta }^{3}-27{\eta }^{4}\right){\dot{r}}^{4}+\left(\frac{3}{4}-\frac{5551}{60}\eta \right.\right.\hfill \\ \hfill & \quad \left.\left.-\frac{256}{3}{\eta }^{2}+\frac{4459}{16}{\eta }^{3}+\frac{569}{4}{\eta }^{4}\right){\dot{r}}^{2}{v}^{2}+\left(\frac{345}{16}-\frac{65\,491}{960}\eta +\frac{12\,427}{96}{\eta }^{2}\right.\right.\hfill \\ \hfill & \quad \left.\left.-\frac{3845}{32}{\eta }^{3}+\frac{585}{8}{\eta }^{4}\right){v}^{4}\right]\hfill \end{align} \tag{ A.19 }$$

$$\begin{align}\hfill {\mathcal{L}}^{(3)}& =\frac{{m}^{3}}{{r}^{3}}\left\{\left[\frac{7}{2}+\left(\frac{7775\,977}{16\,800}+\frac{426}{256}{\pi }^{2}\right)\eta +\frac{121\,449}{560}{\eta }^{2}-\frac{1025}{8}{\eta }^{3}-47{\eta }^{4}\right]{\dot{r}}^{2}+\left[\frac{91}{4}\right.\right.\hfill \\ \hfill & \quad -\left.\left.\left(\frac{13\,576\,009}{50\,400}-\frac{469}{256}{\pi }^{2}\right)\eta +\left(\frac{276\,433}{5040}-\frac{41}{16}{\pi }^{2}\right){\eta }^{2}+\frac{637}{8}{\eta }^{3}-15{\eta }^{4}\right]{v}^{2}\right\}\hfill \end{align} \tag{ A.20 }$$

$$\begin{equation}{\mathcal{L}}^{(4)}=\frac{{m}^{4}}{{r}^{4}}\left[\frac{15}{8}+\left(\frac{3809\,041}{25\,200}-\frac{85}{8}{\pi }^{2}\right)\eta -\left(\frac{20\,131}{420}-\frac{663}{32}{\pi }^{2}\right){\eta }^{2}-\frac{15}{4}{\eta }^{3}+{\eta }^{4}\right]\end{equation} \tag{ A.21 }$$

As it is found in [10, 11, 16], we can write the 3rd PN order of the Ω₁ expression in the center of mass frame, as follows

$$\begin{equation}{{\Omega}}_{1}=\left[\frac{1}{{c}^{2}}({\alpha }_{\text{PN}}+{\alpha }_{\text{PN,Kidder}})+\frac{1}{{c}^{4}}({\alpha }_{2\text{PN}}+{\alpha }_{2\text{PN,Blanchet}})+\frac{1}{{c}^{6}}{\alpha }_{3\text{PN}}+\mathcal{O}\left(\frac{1}{{c}^{7}}\right)\right].\end{equation} \tag{ B.1 }$$

where the individual corrections can be written as

$$\begin{equation}{\alpha }_{\text{PN}}=\frac{Gm}{{r}^{2}}\left(\frac{3}{4}+\frac{1}{2}\eta -\frac{3}{4}\frac{\delta m}{m}\right),\end{equation} \tag{ B.2 }$$

$$\begin{equation}{\alpha }_{\text{PN,Kidder}}=\frac{G}{{r}^{3}}\frac{4+3{\zeta }_{1}}{2}\end{equation} \tag{ B.3 }$$

$$\begin{align}\hfill {\alpha }_{2\text{PN}}& =\frac{Gm}{{r}^{2}}\left[\left(-\frac{3}{2}\eta +\frac{3}{4}{\eta }^{2}-\frac{3}{2}\eta \frac{\delta m}{m}\right){\dot{r}}^{2}+\left(\frac{1}{16}+\frac{11}{8}\eta -\frac{3}{8}{\eta }^{2}\right.\right.\hfill \\ \hfill & \quad \left.\left.+\frac{\delta m}{m}\left(-\frac{1}{16}+\frac{1}{2}\eta \right)\right){v}^{2}\right]+\frac{{G}^{2}{m}^{2}}{{r}^{3}}\hfill \\ \hfill & \quad \times \left(-\frac{1}{4}-\frac{3}{8}\eta +\frac{1}{2}{\eta }^{2}+\frac{\delta m}{m}\left(\frac{1}{4}-\frac{1}{8}\eta \right)\right),\hfill \end{align} \tag{ B.4 }$$

$$\begin{align}\hfill {\alpha }_{2\text{PN,Blanchet}}& =\frac{G{m}_{1}^{2}}{8{m}^{2}{r}^{3}}\left[\right.(16+22{\zeta }_{1}+10{\zeta }_{1}^{2}+{\zeta }_{1}^{3}){v}^{2}-2\left(4+7{\zeta }_{1}+7{\zeta }_{1}^{2}\right.\hfill \\ \hfill & \quad \left.\left.+2{\zeta }_{1}^{3}\right)\frac{Gm}{r}-6(4+3{\zeta }_{1}){\dot{r}}^{2}\right]\hfill \end{align} \tag{ B.5 }$$

$$\begin{align}\hfill {\alpha }_{3\text{PN}}& =\frac{Gm}{{r}^{2}}\left[\left(\frac{15}{8}\eta -\frac{195}{32}{\eta }^{2}+\frac{15}{16}{\eta }^{3}+\frac{\delta m}{m}\left(\frac{15}{8}\eta -\frac{75}{32}{\eta }^{2}\right)\right){\dot{r}}^{4}\right.\hfill \\ \hfill & \quad \left.+\left(-3\eta +\frac{291}{32}{\eta }^{2}-\frac{45}{16}{\eta }^{3}+\frac{\delta m}{m}\left(-3\eta +\frac{177}{32}{\eta }^{2}\right)\right){\dot{r}}^{2}{v}^{2}\right.\hfill \\ \hfill & \quad \left.+\left(\frac{1}{32}+\frac{19}{16}\eta -\frac{31}{8}{\eta }^{2}+\frac{17}{16}{\eta }^{3}+\frac{\delta m}{m}\left(-\frac{1}{32}+\frac{3}{4}\eta -\frac{11}{8}{\eta }^{2}\right)\right){v}^{4}\right]\hfill \\ \hfill & \quad +\frac{{G}^{2}{m}^{2}}{{r}^{3}}\left[\left(\frac{1}{4}-\frac{525}{32}\eta -\frac{159}{16}{\eta }^{2}+\frac{13}{4}{\eta }^{3}+\frac{\delta m}{m}\left(-\frac{1}{4}-\frac{75}{32}\eta \right.\right.\right.\hfill \\ \hfill & \quad \left.\left.\left.-\frac{87}{16}{\eta }^{2}\right)\right){\dot{r}}^{2}+\left(\frac{3}{16}+\frac{27}{4}\eta +\frac{75}{32}{\eta }^{2}-\frac{9}{8}{\eta }^{3}\right.\right.\hfill \\ \hfill & \quad \left.\left.+\frac{\delta m}{m}\left(-\frac{3}{16}+\frac{9}{8}\eta +\frac{35}{32}{\eta }^{2}\right)\right){v}^{2}\right]+\frac{{G}^{3}{m}^{3}}{{r}^{4}}\hfill \\ \hfill & \quad \times \left(\frac{7}{16}-\frac{9}{4}\eta -\frac{9}{8}{\eta }^{2}+\frac{1}{2}{\eta }^{3}+\frac{\delta m}{m}\left(-\frac{17}{16}-\frac{1}{8}\eta -\frac{1}{8}{\eta }^{2}\right)\right).\hfill \end{align} \tag{ B.6 }$$

Here δm = m₁ − m₂ is the mass difference and ζ₁ = m₁/m₂. We can get the expressions for Ω₂ with a simple δm → −δm and ζ₁ → ζ₂ = m₂/m₁ exchange.