General relativistic effect on Sitnikov Three-Body Problem: Restricted Case

1. Introduction

In 1911, MacMillan [1] demonstrated the existence of a periodic solution to the Newtonian restricted three-body problem under the following configuration: the two central primaries (a binary system) have equal mass m and revolve around their barycenter in circular orbits, while the third body is a test particle that moves along a trajectory perpendicular to the orbital plane of the binary and passes through their common barycenter, under the influence of Newtonian gravitational attraction. In 1961, within the same configuration, Sitnikov [2] showed that a periodic solution also exists when the two primaries move along elliptical orbits, see Figure 1. This type of Newtonian restricted three-body problem is often referred to as the “Sitnikov problem.”

Although the Sitnikov three-body problem is one of the simplest dynamical models, it is known to exhibit a wide range of chaotic behaviors. As a result, the Sitnikov problem has been extensively investigated both analytically and numerically; see, for example, [3,4,5,6,7,8,9,10,11,12,13,14] and references therein.

In many previous studies, the three-body problem, including the Sitnikov problem, has been investigated under Newtonian gravity. However, in recent years, with the successful detection of gravitational waves and increasing interest in black hole formation, general relativistic effect in the three-body problem have attracted growing attention.

The general relativistic effect on Euler’s collinear solution was investigated in [15], and the papers e.g., [16,17] studied the contribution of general relativity to Lagrange’s triangle solution, and the stability of the Lagrange points $L_4$ and $L_5$ was examined through linear analysis by [18]. Further, It has been shown that the figure-eight orbit of the three-body system, discovered by Moore [19] under Newtonian gravity (see also [20,21,22,23]) exists under general relativistic gravity, see e.g., [24,25]. And, some papers such as [26,27,28] considered gravitational wave emission in the three-body problem forming Lagrange’s triangular solution.

On the other hand, The discovery of PSR J0337+1715 has drawn increasing attention to test of general relativity in hierarchical triple systems [29]. The effect of general relativity in hierarchical three-body systems has been actively discussed, for example, [30,31,32,33,34,35] and references therein.

The effect of general relativity on the Sitnikov problem has been considered by [36,37]. Kovács et al. [36] explore the Sitnikov three-body problem within the first post-Newtonian approximation to assess how relativistic corrections influence its chaotic phase-space structure. Their numerical analysis shows that increasing the dimensionless gravitational radius $\lambda$ leads to bifurcations in phase-space; altering chaotic regions and escape dynamics, and results in systematic changes to the chaotic saddle, including shifts in the positions of KAM islands and escape corridors. This work demonstrates that even weak relativistic effects can significantly impact chaotic scattering in the classical Sitnikov model.

While, Bernal et al. [37] investigate a relativistic Sitnikov problem using a first post-Newtonian approximation. They explore how the dimensionless gravitational radius $\lambda$ affects asymptotic behavior and chaotic scattering in the system. They identify (i) a metamorphosis of the KAM islands and associated escape regions as $\lambda$ increases, and (ii) two distinct inflection points in the unpredictability of the final state at $\lambda \simeq 0.02$ and $\lambda \simeq 0.028$. These inflection features are quantified using basin entropy analysis, providing deeper insight into relativistic chaotic scattering in the Sitnikov framework.

In this paper, we numerically investigate the Sitnikov three-body problem within the framework of the first post-Newtonian approximation and present practical and physically concrete results, rather than focusing on chaotic behavior and phase-space perspectives mainly discussed by [36,37].

This paper is organized as follows. In Section 2, we derive the general relativistic equations of motion used in this study from the Einstein-Infeld-Hoffmann equations. In Section 3, we present the results of numerical integration of these equations of motion. Section 4 is devoted to conclusions.

2. Equation of Motion

First, we derive the equations of motion used in this paper from the Einstein–Infeld–Hoffmann (EIH) equations, which describe the motion of N point masses within the first post-Newtonian approximation. The Einstein–Infeld–Hoffmann equations are given by: see e.g., [38,39],

$$\begin{array}{cc}\hfill {\displaystyle \frac{d{\mathit{v}}_k}{dt}}& =-\sum _{i\ne k}{\displaystyle \frac{G{m}_i}{r_{ki}^3}}{\mathit{r}}_{ki}+{\displaystyle \frac{1}{c^2}}\left\{\sum _{i\ne k}{\displaystyle \frac{G{m}_i}{r_{ki}^3}}{\mathit{r}}_{ki}\left[\sum _{j\ne i}{\displaystyle \frac{G{m}_j}{r_{ij}}}+4\sum _{j\ne k}{\displaystyle \frac{G{m}_j}{r_{kj}}}\right.\right.\hfill \\ \hfill & -\left.{\displaystyle \frac{1}{2}}\sum _{j\ne i}{\displaystyle \frac{G{m}_j}{r_{ij}^3}}({\mathit{r}}_{ki}·{\mathit{r}}_{ij})+{\displaystyle \frac{3}{2}}{\displaystyle \frac{{({\mathit{r}}_{ki}·{\mathit{v}}_i)}^2}{r_{ki}^2}}-2{\mathit{v}}_i·{\mathit{v}}_i-{\mathit{v}}_k·{\mathit{v}}_k+4{\mathit{v}}_k·{\mathit{v}}_i\right]\hfill \\ \hfill & +\left.\sum _{i\ne k}{\displaystyle \frac{G{m}_i}{r_{ki}^3}}{\mathit{v}}_{ki}[(4{\mathit{v}}_k-3{\mathit{v}}_i)·{\mathit{r}}_{ki}]-{\displaystyle \frac{7}{2}}\sum _{i\ne k}{\displaystyle \frac{G{m}_i}{r_{ki}}}\sum _{j\ne i}{\displaystyle \frac{G{m}_j}{r_{ij}^3}}{\mathit{r}}_{ij}\right\}+\mathcal{O}\left(c^{-4}\right),\hfill \end{array}$$

(1)

where G is the Newtonian gravitational constant, c is the speed of light in vacuum, $m_i$ is the mass of body i, ${\mathit{r}}_i$ is the position vector of body i in the barycentric coordinates system, ${\mathit{r}}_{ij}={\mathit{r}}_i-{\mathit{r}}_j$, $r_{ij}=\left|{\mathit{r}}_{ij}\right|$, ${\mathit{v}}_i=d{\mathit{r}}_i/dt$ is the coordinate velocity of body i (t is the coordinate time), and ${\mathit{v}}_{ij}={\mathit{v}}_i-{\mathit{v}}_j$. From this equation, we derive the equations of motion for the central binary system and the third body.

2.1. Motion of primary bodies

In this paper, we assume that the third body has an infinitesimally small mass, setting $m_3=0$. Then, the motion of the central binary, consisting of two massive primary bodies (1 and 2), can be described by the relativistic two-body problem. We introduce the relative coordinate,

$$\begin{array}{c}\hfill \mathit{R}={\mathit{r}}_1-{\mathit{r}}_2,\mathit{V}={\displaystyle \frac{d\mathit{R}}{dt}}={\mathit{v}}_1-{\mathit{v}}_2,\end{array}$$

(2)

and from Eq. (1), The equation of motion of the binary becomes,

$$\begin{array}{cc}\hfill {\displaystyle \frac{d\mathit{V}}{dt}}& =-{\displaystyle \frac{GM}{R^3}}\mathit{R}\hfill \\ \hfill & +{\displaystyle \frac{GM}{c^2{R}^3}}\{\left[{\displaystyle \frac{GM}{R}}(4+2\nu )+{\displaystyle \frac{3}{2}}\nu {\displaystyle \frac{{(\mathit{R}·\mathit{V})}^2}{R^2}}-(1+3\nu )(\mathit{V}·\mathit{V})\right]\mathit{R}\hfill \\ \hfill & +(4-2\nu )(\mathit{R}·\mathit{V})\mathit{V}\}\hfill \end{array}$$

(3)

in which $R=\left|\mathit{R}\right|$, $M=m_1+m_2$, and $\nu =m_1{m}_2/M^2$. In the Sitnikov problem, the central binary stars have equal masses, and we set $m_1=m_2=m$, so that $M=2m$ and $\nu =1/4$. Thus, Eq. (3) can be rewritten as

$$\begin{array}{c}\hfill {\displaystyle \frac{d\mathit{V}}{dt}}=-{\displaystyle \frac{2Gm}{R^3}}\mathit{R}+{\displaystyle \frac{2Gm}{c^2{R}^3}}\left\{\left[9{\displaystyle \frac{Gm}{R}}+{\displaystyle \frac{3}{8}}{\displaystyle \frac{{(\mathit{R}·\mathit{V})}^2}{R^2}}-{\displaystyle \frac{7}{4}}(\mathit{V}·\mathit{V})\right]\mathit{R}+{\displaystyle \frac{7}{2}}(\mathit{R}·\mathit{V})\mathit{V}\right\}.\end{array}$$

(4)

Here, we transform the coordinate time t, $\mathit{R}$ and $\mathit{V}$ into a dimensionless variables T, $\overline{\mathit{R}}$ and $\overline{\mathit{V}}$ as foloows,

$$\begin{array}{c}\hfill t={\displaystyle \frac{1}{n}}T,\mathit{R}=a\overline{\mathit{R}},R=a\overline{R},\mathit{V}=na\overline{\mathit{V}},\overline{\mathit{V}}={\displaystyle \frac{d\overline{\mathit{R}}}{dT}},n=\sqrt{{\displaystyle \frac{2Gm}{a^3}}},\end{array}$$

(5)

where a is the semi-major axis and n is the mean motion of central binary. Further let us introduce dimensionless gravitational radius $\lambda$ as foloows,

$$\begin{array}{c}\hfill \lambda \equiv {\displaystyle \frac{2Gm}{a{c}^2}},\end{array}$$

(6)

Then, $\lambda =0$ corresponds to the Newtonian case. The value of $\lambda$ increases as the mass of the central bodies increases or as the two primary bodies come closer together, and represents the strength of relativistic effect.

Using Eqs. (5) and (6), Eq. (4) is rewritten as,

$$\begin{array}{c}\hfill {\displaystyle \frac{d\overline{\mathit{V}}}{dT}}=-{\displaystyle \frac{\overline{\mathit{R}}}{{\overline{R}}^3}}+{\displaystyle \frac{\lambda }{{\overline{R}}^3}}\left\{\left[{\displaystyle \frac{9}{2}}{\displaystyle \frac{1}{\overline{R}}}+{\displaystyle \frac{3}{8}}{\displaystyle \frac{{(\overline{\mathit{R}}·\overline{\mathit{V}})}^2}{{\overline{R}}^2}}-{\displaystyle \frac{7}{4}}(\overline{\mathit{V}}·\overline{\mathit{V}})\right]\overline{\mathit{R}}+{\displaystyle \frac{7}{2}}(\overline{\mathit{R}}·\overline{\mathit{V}})\overline{\mathit{V}}\right\}.\end{array}$$

(7)

Eq. (7) is more suitable for numerical integration than Eq. (4) since the variables $\overline{\mathit{R}},\overline{\mathit{V}}$ are of order $\mathcal{O}\left(1\right)$.

2.2. Motion of third body

Next, let us write down the equation of motion of the third body (test particle). In the Sitnikov problem, the test particle moves along a trajectory perpendicular to the orbital plane of the central binary and passes through the binary’s center of mass. We therefore assume that the motion of the primary bodies is confined to the $\overline{X}\overline{Y}$-plane. And because $m_1=m_{2}2=m$, we set the positions and velocities of the primaries as follows:

$$\begin{array}{c}\hfill {\mathit{r}}_1=-{\mathit{r}}_2={\displaystyle \frac{1}{2}}\mathit{R},\mathit{R}=(X,Y,0),{\mathit{v}}_1=-{\mathit{v}}_2={\displaystyle \frac{1}{2}}\mathit{V},\mathit{V}=V_X,V_Y,0).\end{array}$$

(8)

Since, in this paper, we assume the third body is a test particle with $m_3=0$, we set the position and velocity of the third body as follows, respectively:

$$\begin{array}{c}\hfill {\mathit{r}}_3=(0,0,z),{\mathit{v}}_3=(0,0,v_z),\end{array}$$

(9)

From Eq. (1), the equation of motion of third body becomes,

$$\begin{array}{cc}\hfill {\displaystyle \frac{d{v}_z}{dt}}& =-{\displaystyle \frac{2Gm}{r^3}}z\hfill \end{array}$$

$$\begin{array}{cc}\hfill & +{\displaystyle \frac{Gm}{c^2{r}^3}}\left\{\left[{\displaystyle \frac{5}{2}}{\displaystyle \frac{Gm}{R}}+16{\displaystyle \frac{Gm}{r}}+{\displaystyle \frac{3}{16}}{\displaystyle \frac{{(\mathit{R}·\mathit{V})}^2}{r^2}}-\mathit{V}·\mathit{V}+6v_z^2\right]z+{\displaystyle \frac{3}{2}}(\mathit{R}·\mathit{V})v_z\right\},\hfill \end{array}$$

(10)

$$\begin{array}{cc}\hfill {\displaystyle \frac{d{v}_x}{dt}}& ={\displaystyle \frac{d{v}_y}{dt}}=0,\hfill \end{array}$$

(11)

where

$$\begin{array}{c}\hfill r=|{\mathit{r}}_{31}|=|{\mathit{r}}_{32}|=\sqrt{z^2+{\left({\displaystyle \frac{R}{2}}\right)}^2}.\end{array}$$

(12)

As (5), we use following transformations,

$$\begin{array}{c}\hfill z=a\overline{z},r=a\overline{r},\overline{r}=\sqrt{{\left({\displaystyle \frac{\overline{R}}{2}}\right)}^2+{\overline{z}}^2},v_z=na{\overline{v}}_z,{\overline{v}}_z={\displaystyle \frac{d\overline{z}}{dT}},\end{array}$$

(13)

and using Eq. (6), Eq. (10) is replaced by,

$$\begin{array}{cc}\hfill {\displaystyle \frac{d{\overline{v}}_z}{dT}}& =-{\displaystyle \frac{\overline{z}}{{\overline{r}}^3}}\hfill \end{array}$$

$$\begin{array}{cc}\hfill & +{\displaystyle \frac{\lambda }{{\overline{r}}^3}}\left\{\left[{\displaystyle \frac{5}{8}}{\displaystyle \frac{1}{\overline{R}}}+{\displaystyle \frac{4}{\overline{r}}}+{\displaystyle \frac{3}{32}}{\displaystyle \frac{{(\overline{\mathit{R}}·\overline{\mathit{V}})}^2}{{\overline{r}}^2}}-{\displaystyle \frac{1}{2}}(\overline{\mathit{V}}·\overline{\mathit{V}})+3{\overline{v}}_z^2\right]\overline{z}+{\displaystyle \frac{3}{4}}(\overline{\mathit{R}}·\overline{\mathit{V}}){\overline{v}}_z\right\},\hfill \end{array}$$

(14)

$$\begin{array}{cc}\hfill {\displaystyle \frac{d{\overline{v}}_x}{dT}}& ={\displaystyle \frac{d{\overline{v}}_y}{dT}}=0.\hfill \end{array}$$

(15)

The configuration of the dynamical system considered in this paper is shown in Figure 1.

3. Numerical experiments

Let us present the results of the numerical integration of Eqs. (7) and (14). The setup of our numerical integration is summarized in Table 1.

The initial condition of the central binary is chosen as follows:

$$\begin{array}{c}\hfill \overline{X}=1-e,\overline{Y}=0,\dot{\overline{X}}=0,\dot{\overline{Y}}={\displaystyle \frac{\sqrt{1-e^2}}{(1-e)}},\end{array}$$

(16)

where e is the orbital eccentricity of the central binary, with $0\le e\le 0.9$, varied in intervals of $\Delta e=0.1$. The initial values of the third body, ${\overline{z}}_0$ and ${\dot{\overline{z}}}_o$, are taken as $0\le {\overline{z}}_o\le 6$ with intervals of $\Delta {\overline{z}}_0=0.05$, and ${\dot{\overline{z}}}_0=0$. The range of the dimensionless gravitational radius $\lambda$ is $0\le \lambda \le 0.035$.

Numerical integration is performed using Gragg’s extrapolation method based on the Aitken–Neville algorithm, and the calculation time per orbit is set to a maximum of 1000 Keplerian periods for the central binary.

When the energy of the test particle, E, becomes positive, the particle is considered to have escaped from the system;

$$\begin{array}{c}\hfill E={\displaystyle \frac{1}{2}}{\overline{v}}_z^2-{\displaystyle \frac{1}{\sqrt{{\displaystyle {\overline{z}}^2+{\left(\frac{\overline{R}}{2}\right)}^2}}}}\le 0.\end{array}$$

(17)

Figure 2 shows the parameter space of $({\overline{z}}_0,\lambda )$, where blue indicates the parameter region where the test particle was ejected from the system, while red depicts the region where a collision with the central binary star occurred. The collision of the central binary is caused by large values of $\lambda$ as well as by the orbital eccentricity e.

Through numerical experiments, it has been found that collisions between the primaries occur when the periastron distance of the binary becomes smaller than several tens of times the Schwarzschild radius.

$$\begin{array}{c}\hfill a(1-e)\lesssim \delta {\displaystyle \frac{2Gm}{c^2}}=\delta a\lambda ,\delta =\mathcal{O}\left(10\right).\end{array}$$

(18)

Then, throughout Figure 2, Figure 3, Figure 4, Figure 5, Figure 6 and Figure 7, we present the results excluding the red-colored parameter region shown in Figure 2.

Figure 3 illustrates the escape time of the test particle from the system. In the case of $e=0.0$ (i.e., the orbit of the primaries is circular), the escape time of the test particle is very long, and it can remain in the system for an extended period. As the dimensionless gravitational radius $\lambda$ and the eccentricity of the binary orbit e increase, the tendency for the test particle to escape from the system becomes more pronounced.

As a general trend, when the orbital eccentricity of the central binary is small, it takes longer for the test particle to escape from the system. However, as the eccentricity increases, the particle escapes in a shorter time. An increase in the orbital eccentricity of the central binary leads to larger oscillations in the relative distance between the two bodies, and it is thought that the periodic variations in the relative distance affect the escape of the test particles.

Figure 4 and Figure 5 show the escape velocity and the $\overline{z}/(\overline{R}/2)$ ratio of the test particle at the moment of escape from the system, respectively. Due to the orbital eccentricity of the binary, the distance between the two primaries is not constant but fluctuates; therefore, the ratio $\overline{z}/(\overline{R}/2)$ is used, instead of $\overline{z}$, to represent the spatial relationship between the test particle and the binary. A small value of this ratio indicates that the test particle escaped near the orbital plane of the binary, and in the case of a large value of the ratio, the test particle leaves the system at a location farther away from the binary.

From Figure 4 and Figure 5, it can be observed that when the eccentricity of the binary orbit is small, the escape velocity of the test particle tends to be relatively low, and the ratio $\overline{z}/(\overline{R}/2)$ tends to be large — that is, the test particle is more likely to escape from a location far from the binary. On the other hand, as the eccentricity of the central binary increases, the test particle tends to escape at a higher velocity from a position closer to the orbital plane of the binary.

Figure 6 shows the acceleration acting on the test particle at the moment it escapes from the system. It can be observed that, as the eccentricity of the binary orbit increases, regions with higher acceleration (indicated in yellow or cyan) appear. However, the overall trend is that the acceleration acting on the test particle is nearly zero when it moves far away from the system.

Figure 7 shows the survival rate of test particles that remain in the system up to 1000 orbital periods of the central binary. Given the eccentricity e of the binary orbit and the dimensionless gravitational radius $\lambda$, the number of initial values of $z_0$ is $N=120$, since $z_0$ varies from 0 to 6 in increments of 0.05. Therefore, Figure 7 represents the ratio $n\left(P\right)/N$ as a function of orbital period of centrap binary, P where $n\left(P\right)$ mewns the number of particles remaining in the system at period P.

For the range $0.0\le e\le 0.4$, the number of surviving test particles decreases with increasing values of both the dimensionless gravitational radius $\lambda$ and the orbital eccentricity of the binary e. However, $e\ge 0.5$, there exists a region where the number of s urviving test particles increases despite the increase in $\lambda$; for example, see the dark-colored area around $\lambda =0.03$ in the upper left of Figure 7, where $e=0.5$.

Looking at Figure 2, the region where the number of surviving particles increases with increasing $\lambda$ is the region where there are many test particles remaining in the system. This region exists near the lower limit of the region where the central binary collision occurs.

Throughout the numerical experiments, it can be seen that general relativistic effect does not bind test particles to the system, but rather acts to eject them from the system. In addition, we can see that the increase in the eccentricity of the central binary orbit is also contributing to this tendency.

4. Conclusions

In this paper, we investigated the effect of general relativity on the Sitnikov problem. To study the general relativistic contributions, we first derived the equations of motion for both the binary and the test particle based on the first post-Newtonian Einstein–Infeld–Hoffmann equation, and integrated these equations numerically. We examine the behavior of the test particle (third body) as a function of the orbital eccentricity of the central binary e, the dimensionless gravitational radius $\lambda$, which characterizes the strength of general relativistic effect, and the initial position of the test particle $z_0$.

Our numerical calculations revealed the following: (a) Due to general relativistic effect and increasing eccentricity of the central binary, the test particle is likely to be ejected from the system. (b) When the binary’s eccentricity is low, it takes longer for the test particle to be ejected from the system; however, as the eccentricity increases, the particle tends to be ejected more quickly. (c) The escape velocity of the test particle tends to increase with the binary’s eccentricity. (d) When the binary’s eccentricity is small, the test particle is ejected farther from the orbital plane of binary; as the eccentricity increases, the particle is likely to be ejected closer to the plane. And (e) although the acceleration acting on the test particle can be seen to increase as the eccentricity and dimensionless gravitational radius become larger, overall, the acceleration acting when the test particle escapes from the system is generally zero.

As mentioned at the end of Section 3, in the Sitnikov three-body problem, the increase both in general relativistic effects and in the eccentricity of the central binary orbit contribute to the ejection of the test particle from the system rather than binding it to the system. Then, as one possible application of the Sitnikov problem, it may be conceivable to explain the mechanism of astrophysical jets, especially Relativistic jets.

The Blandford–Payne mechanism and the Blandford–Znajek mechanism have been proposed as leading theories for the formation of relativistic jets; however, a complete understanding has not yet been achieved. To resolve this issue, it is essential to determine which of the two mechanisms predominantly operates, or whether both do, and to provide detailed observational verification of the black hole spin and magnetic field structure. In addition, a coherent explanation of the continuous process by which the jet is formed, accelerated, and collimated (i.e., from formation to acceleration and stabilization) is required. If jets are emitted from a (nearly) equal-mass binary black hole system, then, as a consequence of the Sitnikov problem, particle motion in the z-direction (perpendicular to the binary orbital plane) may be enhanced. This could potentially contribute to a better understanding of the jet acceleration and collimation processes. This is because, as demonstrated by our numerical calculations, in the Sitnikov problem, relativistic effect and an increase in the orbital eccentricity of binary tend to expel test particles from the system. When ejected, the test particles either possess high velocities or experience nearly zero acceleration, indicating that they are not gravitationally bound.

In this paper, we have considered the third body as a test particle, but if we regard the third body as a massive particle, the central binary star will be affected by the gravity of this third body. Then, the behavior of the system will become more complex. Gravitational waves from such three-body system have not been studied so far, therefore it would be worthwhile to investigate them in the future.