Chaos in a generalized Euler's three-body problem

The Kepler problem in celestial mechanics—two bodies interact with each other by the inverse-square force law—has provided us with many insights into the motion of astrophysical bodies. This problem was generalized by Euler to a certain three-body problem—motion of a particle subjected to the inverse-square force sourced by two fixed masses in space—which is known as Euler's three-body problem [1].

These problems have a common feature of the Liouville integrability. The Kepler problem is integrable, i.e. there exist three constants of motion, energy, angular momentum around an axis, and the squared total angular momentum, commutable with each other by the Poisson brackets. As a result, the equations of motion reduce to decoupled ordinary differential equations for each variables. In addition, there is a vector-type constant of motion in this system, the so-called Laplace–Runge–Lenz vector (see, e.g. reference [2]). Since all six of those constants are Poisson commutable, the Kepler problem is superintegrable and is solved in an algebraic way. On the other hand, the Euler problem admits two constants of motion, energy and angular momentum. Additionally, there is a nontrivial constant of a second-order polynomial in momentum, the so-called Whittaker constant [3, 4]. These three constants are Poisson commutable with each other, and therefore, the system is integrable. The nontrivial constants in both cases are relevant to hidden symmetry developed in general relativity—the Killing tensors (see, e.g. reference [5]). Though the constants are composed by reducible Killing tensors (i.e. the solution to the Killing tensor equations in the Euclidean space is given as a linear combination of the symmetric tensor products of the Killing vectors.) because of the maximal symmetry in the Euclidean space, they still play an important role in terms of the integrability of particle dynamics.

Such nontrivial constants were discovered by exploring the separation of variables of the corresponding Hamilton–Jacobi equation. This insight was applied to the geodesic systems in the Kerr spacetime, and it was found that the separation of variables of the Hamilton–Jacobi equation occurs and there exists the Carter constant [6], which is relevant to the irreducible and nontrivial Killing tensor [7]. This discovery leads us to study particle motion in the Kerr geometry further in analytical ways (e.g. reference [8]).

It is a natural but nontrivial question to ask whether integrability is preserved if the system is generalized somehow. In generalizing an integrable Newtonian particle system to a general relativistic particle system, we often encounter examples where integrability breaks down. For example, the dynamics of a massive particle in the four-dimensional (4D) Majumdar–Papapetrou dihole spacetime [9, 10]—a generalization of the Euler problem to general relativity—is nonintegrable because it exhibits chaotic behavior [11–17]. In other relativistic Euler's three-body problems, chaotic orbits also appear in general [18–20].

Recently, with the progress in the study of higher-dimensional black holes [23], the relationship between the hidden symmetry and the Liouville integrability has been gradually revealed through the dynamics of particles in higher dimensions [24]. One of the generalizations of the Euler problem to higher-dimensional space was also discussed in the context of Newtonian theory, in which the spatial dimension was parameterized while the potential still being inversely proportional to distance [25]. In this case, even if the spatial dimension is arbitrary, there exists the same type of the Whittaker constant; in other words, the integrability is preserved. However, in generalization of Newtonian gravity to higher-dimensional space, the power-law of gravitational force should be modified accordingly. ¹ Therefore, in this paper, we consider the Euler problem with the inverse-square potential, which can be seen as a natural generalization of the three-body problem to higher-dimensional Newtonian theory. The purpose of this paper is to clarify whether the integrable structure of the original Euler problem can be preserved even in such a generalized Euler's three-body problem.

It should be noted that the timelike geodesics in the higher-dimensional Majumdar–Papapetrou dihole spacetime show chaotic behavior [27]. Thus, we may also ask whether the Newtonian limit of this relativistic system recovers integrability. As in the case of the black ring, ² the integrable structure may be recovered in the Newtonian limit. One of our purposes is to fill in the missing pieces and clarifies the boundary of integrability.

To approach this problem, we focus on stable circular massive particle orbits in the 5D Majumdar–Papapetrou dihole spacetime, which are caused by the many-body effect of the sources [29, 30]. Although it should be noted that relativistic corrections in higher dimensions can affect particle dynamics at infinity, ³ this fact suggests that stable circular orbits also exist in Newtonian gravity. If such orbits exist in our system, stable bound orbits inevitably appear near these orbits (discussed in detail below). If the system is nonintegrable, then the chaos of the stable bound orbits can be determined by the Poincaré map method.

This paper is organized as follows. In section 2, we consider the existence of stationary particle orbits in the Newtonian gravitational potential generated by two fixed centers in ${\mathbb{E}}^{4}$. After formulating a method for finding stationary orbits and determining their stability, we show the sequences of stable/unstable stationary orbits. In section 3, we discuss the appearance of chaos for stable bound orbits by using the Poincaré map method. Section 4 is devoted to a summary and discussions.

We focus on the Newtonian potential generated by two point mass sources fixed at different points in ${\mathbb{E}}^{4}$. Let r be a position vector. Let r = ∓ a denote the positions of the sources with mass M_±, respectively, where a is a nonzero constant vector. Then the mass density distribution is

$$\begin{equation}\sigma (\boldsymbol{r})={M}_{+}\delta (\boldsymbol{r}+\boldsymbol{a})+{M}_{-}\delta (\boldsymbol{r}-\boldsymbol{a}),\end{equation} \tag{ 1 }$$

where δ denotes the delta function. Solving the Newtonian field equation with the source term σ( r ), ⁴ we obtain a Newtonian gravitational potential

$$\begin{equation}{\Phi}(\boldsymbol{r})=-\frac{G{M}_{+}}{2\enspace \vert \enspace \boldsymbol{r}+\boldsymbol{a}\enspace {\vert }^{2}}-\frac{G{M}_{-}}{2\enspace \vert \enspace \boldsymbol{r}-\boldsymbol{a}\enspace {\vert }^{2}},\end{equation} \tag{ 3 }$$

where G is the gravitational constant. Hereafter, the masses are assumed to be equal, M_± = M. Introduce the cylindrical coordinates (ρ, θ, ϕ, z) in which the Euclidean metric takes the form

$$\begin{equation}\mathrm{d}{\ell }^{2}=\mathrm{d}{\rho }^{2}+{\rho }^{2}(\mathrm{d}{\theta }^{2}+{\mathrm{sin}}^{2}\enspace \theta \enspace \mathrm{d}{\phi }^{2})+\mathrm{d}{z}^{2}.\end{equation} \tag{ 4 }$$

Without loss of generality, we may put the point sources at z = ±a on the z axis, where a = | a |. Then the potential Φ( r ) in these coordinates is given by

$$\begin{equation}{\Phi}(\boldsymbol{r})=-\frac{GM}{2}\left(\frac{1}{{r}_{+}^{2}}+\frac{1}{{r}_{-}^{2}}\right),\end{equation} \tag{ 5 }$$

where

$$\begin{equation}{r}_{{\pm}}=\sqrt{{(z{\pm}a)}^{2}+{\rho }^{2}}.\end{equation} \tag{ 6 }$$

Let us consider freely falling particle motion in Φ( r ). Let m be particle mass and let p be a canonical momentum conjugate with coordinates of the particle. The Hamiltonian of this mechanical system is given by

$$\begin{equation}H=\frac{\vert \enspace \boldsymbol{p}\enspace {\vert }^{2}}{2m}+m{\Phi}\end{equation} \tag{ 7 }$$

$$\begin{equation}=\frac{1}{2m}\left({p}_{z}^{2}+{p}_{\rho }^{2}+\frac{{Q}^{2}}{{\rho }^{2}}\right)-\frac{\alpha \enspace m}{2}\left(\frac{1}{{r}_{+}^{2}}+\frac{1}{{r}_{-}^{2}}\right),\end{equation} \tag{ 8 }$$

where α = GM and Q² is defined by

$$\begin{equation}{Q}^{2}={p}_{\theta }^{2}+\frac{{p}_{\phi }^{2}}{{\mathrm{sin}}^{2}\enspace \theta },\end{equation} \tag{ 9 }$$

which is a constant of motion associated with the S² rotational symmetry of Φ. We use units in which m = 1 in what follows. The Hamiltonian H is equivalent to particle energy and takes a constant value E. The energy conservation H = E leads to the energy equation

$$\begin{equation}\frac{1}{2}({\dot {z}}^{2}+{\dot {\rho }}^{2})+V=E,\end{equation} \tag{ 10 }$$

$$\begin{equation}V(\rho ,z;{Q}^{2})=\frac{{Q}^{2}}{2{\rho }^{2}}-\frac{\alpha }{2}\left(\frac{1}{{r}_{+}^{2}}+\frac{1}{{r}_{-}^{2}}\right),\end{equation} \tag{ 11 }$$

where the dots denote the derivatives with respect to time, and the Hamilton equations, ${p}_{z}=\dot {z}$ and ${p}_{\rho }=\dot {\rho }$, have been used. We call V the effective potential in what follows.

We focus on stationary orbits where the ρ and z coordinates of a particle remain constant. Note that all such orbits are circular because of the S² rotational symmetry of Φ. To move on the circular orbits, a particle must stay at a stationary point of the effective potential V, where the following conditions must hold:

$$\begin{equation}{V}_{z}=\alpha \left(\frac{z+a}{{r}_{+}^{4}}+\frac{z-a}{{r}_{-}^{4}}\right)=0,\end{equation} \tag{ 12 }$$

$$\begin{equation}{V}_{\rho }=-\frac{{Q}^{2}}{{\rho }^{3}}+\alpha \rho \left(\frac{1}{{r}_{+}^{4}}+\frac{1}{{r}_{-}^{4}}\right)=0,\end{equation} \tag{ 13 }$$

$$\begin{equation}V=E,\end{equation} \tag{ 14 }$$

where V_i = ∂_i V (i = z, ρ). The condition (12) leads to the three real roots

$$\begin{equation}z=0,\end{equation} \tag{ 15 }$$

$$\begin{equation}z={\pm}{z}_{0}(\rho ){:=}{\pm}\sqrt{\frac{3{a}^{2}-{\rho }^{2}}{1+2a/\sqrt{{\rho }^{2}+{a}^{2}}}},\end{equation} \tag{ 16 }$$

and the others are imaginary roots. Note that z₀ is well-defined only in the range $0{< }\rho {\leqslant}\sqrt{3}a$. Furthermore, we obtain the following values ${Q}_{0}^{2}$ and E₀ by solving equations (13) and (14) for Q² and E:

$$\begin{equation}{Q}^{2}={Q}_{0}^{2}{:=}\alpha {\rho }^{4}\left(\frac{1}{{r}_{+}^{4}}+\frac{1}{{r}_{-}^{4}}\right),\end{equation} \tag{ 17 }$$

$$\begin{equation}E={E}_{0}{:=}V(\rho ,z;{Q}_{0}^{2})=\frac{\alpha {\rho }^{2}}{2}\left(\frac{1}{{r}_{+}^{4}}+\frac{1}{{r}_{-}^{4}}\right)-\frac{\alpha }{2}\left(\frac{1}{{r}_{+}^{2}}+\frac{1}{{r}_{-}^{2}}\right).\end{equation} \tag{ 18 }$$

Since ${Q}_{0}^{2}$ is always non-negative, we can find a stationary orbit at any point on the sequences

$$\begin{equation}\left.{\gamma }_{0}{:=}\left\{\enspace (\rho ,z){\backslash}(0,{\pm}a)\enspace \right\vert \enspace z=0\quad \mathrm{o}\mathrm{r}\quad z={\pm}{z}_{0}\enspace \right\}.\end{equation} \tag{ 19 }$$

Let us divide γ₀ into two parts, depending on whether the stationary orbit at each point on γ₀ is stable or unstable. If V at an extremum point on γ₀ is locally minimized, then the circular orbit is stable. If not, that is, if it is locally maximized or has a saddle point, then the circular orbit is unstable. Let (V_ij ) be the Hessian matrix of V, where V_ij (ρ, z; Q²) := ∂_j ∂_i V (i, j = ρ, z). Let h and k be the determinant and the trace of (V_ij ), respectively, i.e. h(ρ, z; Q²) :=  det(V_ij ) and k(ρ, z; Q²) := tr(V_ij ). Evaluating h and k at ${Q}^{2}={Q}_{0}^{2}$, we obtain

$$\begin{align}{h}_{0}(\rho ,z)& {:=}h(\rho ,z;{Q}_{0}^{2})=\frac{64\enspace {\alpha }^{2}{a}^{2}{\rho }^{2}}{{r}_{+}^{6}{r}_{-}^{6}}-12{\alpha }^{2}\left(\frac{1}{{r}_{+}^{4}}+\frac{1}{{r}_{-}^{4}}\right)\left[\frac{{(z+a)}^{2}}{{r}_{+}^{6}}+\frac{{(z-a)}^{2}}{{r}_{-}^{6}}\right],\end{align} \tag{ 20 }$$

$$\begin{equation}{k}_{0}(\rho ,z){:=}k(\rho ,z;{Q}_{0}^{2})=\alpha \left(\frac{1}{{r}_{+}^{4}}+\frac{1}{{r}_{-}^{4}}\right),\end{equation} \tag{ 21 }$$

where ${V}_{ij,0}(\rho ,z){:=}{V}_{ij}(\rho ,z;{Q}_{0}^{2})$ are given by

$$\begin{equation}{V}_{\rho \rho ,0}=4\enspace \alpha \left(\frac{1}{{r}_{+}^{4}}+\frac{1}{{r}_{-}^{4}}\right)-4\enspace \alpha {\rho }^{2}\left(\frac{1}{{r}_{+}^{6}}+\frac{1}{{r}_{-}^{6}}\right),\end{equation} \tag{ 22 }$$

$$\begin{equation}{V}_{zz,0}=\alpha \left(\frac{1}{{r}_{+}^{4}}+\frac{1}{{r}_{-}^{4}}\right)-4\enspace \alpha \left[\enspace \frac{{(z+a)}^{2}}{{r}_{+}^{6}}+\frac{{(z-a)}^{2}}{{r}_{-}^{6}}\enspace \right],\end{equation} \tag{ 23 }$$

$$\begin{equation}{V}_{\rho z,0}={V}_{z\rho ,0}=-4\enspace \alpha \rho \left(\frac{z+a}{{r}_{+}^{6}}+\frac{z-a}{{r}_{-}^{6}}\right).\end{equation} \tag{ 24 }$$

Note that k₀ > 0 everywhere. Hence, the circular orbits on γ₀ are stable if h₀ > 0, unstable if h₀ < 0, and marginally stable if h₀ = 0.

We consider the stability of the circular orbits on ${z}^{2}={z}_{0}^{2}$. The function h₀ on this branch satisfies

$$\begin{equation}{h}_{0}(\rho ,{\pm}{z}_{0})={\alpha }^{2}\frac{2\enspace {a}^{3}+(2\enspace {a}^{2}-{\rho }^{2})\sqrt{{\rho }^{2}+{a}^{2}}}{16\enspace {a}^{4}{\rho }^{4}{({\rho }^{2}+{a}^{2})}^{3/2}}{\geqslant}0,\end{equation} \tag{ 25 }$$

where the equality holds for $\rho =\sqrt{3}a$. This result indicates that all the circular orbits are stable on this branch. Next we consider the stability of the circular orbits on z = 0, where h₀ reduces to

$$\begin{equation}{h}_{0}(\rho ,0)=16\enspace {\alpha }^{2}{a}^{2}\frac{{\rho }^{2}-3\enspace {a}^{2}}{{({\rho }^{2}+{a}^{2})}^{6}}.\end{equation} \tag{ 26 }$$

This implies that the circular orbits are stable (i.e. h₀ > 0) for $\rho { >}\sqrt{3}a$, unstable (i.e. h₀ < 0) for $\rho {< }\sqrt{3}a$, and marginally stable (i.e. h₀ = 0) for $\rho =\sqrt{3}a$. Note that this stability behavior is caused by switching of the stability in the z direction at $\rho =\sqrt{3}a$,

$$\begin{equation}{V}_{zz,0}(\rho ,0)=2\enspace \alpha \frac{{\rho }^{2}-3\enspace {a}^{2}}{{({\rho }^{2}+{a}^{2})}^{3}}.\end{equation} \tag{ 27 }$$

The sequence of the circular orbits and their stability are shown in figure 1. Blue solid curves show the sequence of stable circular orbits, and blue dashed segment shows the sequence of unstable circular orbits.

Zoom In Zoom Out Reset image size

It is worth noting that while there is no stable circular orbit for a single point mass source in ${\mathbb{E}}^{4}$, there are stable circular orbits for any nonzero value of a (up to the asymptotic region). To see the effect of a on the existence of stable circular orbits, we consider V on z = 0,

$$\begin{equation}V(\rho ,0)=-\frac{\alpha }{{\rho }^{2}+{a}^{2}}+\frac{{Q}^{2}}{2{\rho }^{2}}.\end{equation} \tag{ 28 }$$

If Q² ⩾ 2α, the potential V(ρ, 0) decreases monotonically as ρ increases, which implies that no local minimum exists. If Q² < 2α, however, it always has a local minimum at $\rho =a{[Q/(\sqrt{2\alpha }-Q)]}^{1/2}$, which is caused by V → ∞ as ρ → 0 and V↗0 as ρ → ∞. In the limit $Q\to \sqrt{2\alpha }$, the local minimum point goes to infinity. These results imply that the existence of nonzero a due to the two-body effect is essential for the potential (28) to form a local minimum point. Therefore, the stable circular orbit can be interpreted as a result of the many-body effect.

We use stable stationary orbits to find stable bound orbits—a particle moves in a spatially bounded region without reaching infinity or the sources even if small perturbations are applied. As discussed in the previous section, a particle on a stable stationary orbit must stay at a local minimum point of the effective potential V. If the energy level E₀ at the local minimum point increases slightly (i.e. some positive energy ΔE > 0 is injected), then it will start to move away from the local minimum. If the energy displacement |ΔE| is small enough, the effective potential contour at E₀ + ΔE will be closed. Thus, the particle with energy E₀ + ΔE oscillates in the vicinity of the local minimum and is bounded inside the contour, i.e. the particle moves on a stable bound orbit. Therefore, the stationary stable orbits inevitably induce the existence of stable bound orbits in its vicinity.

Figure 2 shows some typical effective potential contours in the upper panels. We use units in which α = 1 and a = 1 in what follows. Black solid curves denote the contours of V with ${Q}^{2}={Q}_{0}^{2}({\rho }_{0},0)$, where (a) ρ₀ = 2.0, (b) ρ₀ = 1.9, and (c) ρ₀ = 1.8. The position (ρ, z) = (ρ₀, 0) corresponds to the local minimum point of V, at which the values of E₀ and ${Q}_{0}^{2}$ are evaluated as (a) $({E}_{0},{Q}_{0}^{2})=(-0.04,1.28)$, (b) $({E}_{0},{Q}_{0}^{2})\simeq (-0.047\enspace 05,1.226)$, and (c) $({E}_{0},{Q}_{0}^{2})\simeq (-0.055\enspace 62,1.167)$. Red solid curves show the contour level with V = −0.025. Blue solid curves denote stable bound orbits with E = −0.025 and ${Q}^{2}={Q}_{0}^{2}({\rho }_{0},0)$, which are confined in the red closed curves. In the case (a), the stable bound orbit shows ordered patterns inside the red closed curve. These are formed by the superposition of oscillations in two directions, a Lissajous curve. In the case (b)—the local minimum point is a bit closer to the sources than the case (a)—the stable bound orbit shows a slightly disturbed from a Lissajous curve. This can be seen as evidence that the chaotic nature of the orbit begins to emerge. In the case (c)—the local minimum point is the closest among these three cases—the stable bound orbit denotes no longer a Lissajous curve. The completely irregular pattern suggests the emergence of chaotic nature of the orbits.

Zoom In Zoom Out Reset image size

Using the Poincaré map method, ⁵ we evaluate more quantitatively the signs of chaos in stable bound orbits. The lower panels of figure 2 show the Poincaré maps of stable bound orbits for 50 different initial conditions (denoted by different colors) with the same energy and the angular momentum as in the upper case. The Poincaré section is defined at z = 0, where ρ and p_ρ are recorded when a particle passes through it with $\dot {z}{ >}0$. In the case (a), we find closed dotted loops in the (ρ, p_ρ ) plane for each stable bound orbit. This implies that the orbits in the phase space lie on a torus, which means that the chaotic nature of the orbits is not manifest. In the case (b), the Poincaré sections form closed dotted curves in some cases while most are scattered to fill bounded regions in the (ρ, p_ρ ) plane. The latter behavior is a result of the chaotic nature of stable bound orbits. In the case (c), the point set that forms a closed dotted curve no longer appears, and all orbits show chaos. Therefore, we can conclude that our generalized Euler's three-body problem shows chaos in general.

We have considered the generalized Euler's three-body problem, the dynamics of a freely falling particle in the Newtonian gravitational potential generated by point masses fixed, respectively, at two different points in ${\mathbb{E}}^{4}$. We presented the conditions under which particles remain in stationary orbits in terms of the effective potential and proposed a systematic procedure for solving them. We used it to identify the sequences of stationary circular orbits in the generalized Euler problem. Furthermore, the linear stability of each circular orbit was clarified. As a result, we have found a family of stable circular orbits extending from two point masses to infinity, whose existence is independent of the separation between the sources. It is worth noting that stable stationary orbits do not exist in the single point source case, i.e. the Kepler problem in ${\mathbb{E}}^{4}$ [35]. Therefore, we can conclude that the existence of stable circular orbits in this system is caused by the many-body effect of the sources.

We compare our results with the sequences of stable circular orbits in the 5D Majumdar–Papapetrou dihole spacetime with equal mass M_* and separation a_* [29]. In the vicinity of each pair of sources, we can see the difference between Newtonian gravity, where there exist sequences of stable circular orbits up to an arbitrary neighborhood of the point masses, and general relativity, where a pair of the innermost stable circular orbits appear due to the relativistic effects. On the other hand, in the asymptotic region, stable circular orbits exist when the dihole separation is large (${a}_{{\ast}}/\sqrt{{M}_{{\ast}}}{\geqslant}\sqrt{3}$ in geometrized units), which is consistent with the result in our Euler problem. However, our results do not reproduce the disappearance of stable circular orbits in the asymptotic region for small separation case (${a}_{{\ast}}/\sqrt{{M}_{{\ast}}}{< }\sqrt{3}$). This difference is due to a higher-order multipole, which contains a relativistic correction in this case, affecting the existence of stable circular orbits in the asymptotic region of higher-dimensional spacetime. Such relativistic effects on particle dynamics at infinity are unique to higher-dimensional spacetimes.

In the latter part of this paper, we used the stable stationary orbits to find stable bound orbits. Calculating the Poincaré maps for the stable bound orbits, we revealed the appearance of chaotic behavior. Furthermore, we found that the chaos increases as the region where the stable bound orbits lie approaches the sources. These results imply that our generalized Euler problem is nonintegrable. In other words, there is no additional constant of motion such as the Whittaker constant. Therefore, we conclude that the chaos observed in the timelike geodesics of the 5D Majumdar–Papapetrou dihole spacetimes does not disappear even in the Newtonian limit.

The fact that our system has no integrable structure is preserved even if we parameterize the spatial dimension without changing the inverse-square potential. On the other hand, if the exponent of the inverse power-law of the potential is larger than 2, we need further study. In such cases, however, the Poincaré map method may not be applicable because no stable bound orbits are expected to exist. Euler's three-body problem in a higher-dimensional Euclidean space with compactified extra-dimensional space is also an interesting problem for the future.

This work was supported by Grant-in-Aid for Early-Career Scientists from the Japan Society for the Promotion of Science (JSPS KAKENHI Grant No. JP19K14715).

All data that support the findings of this study are included within the article (and any supplementary files).