A&A 394, 323-328 (2002)
DOI: 10.1051/0004-6361:20021144

Symmetric doubly asymptotic orbits at collinear equilibrium points in the general three-body problem

E. A. Perdios - V. S. Kalantonis

Department of Engineering Sciences, University of Patras, 26500 Patras, Greece

Received 22 March 2002 / Accepted 23 July 2002

Abstract
This paper concerns the numerical determination of doubly asymptotic orbits at the collinear equilibrium points of the general three-body problem. To compute a doubly asymptotic orbit we use a linear analysis to determine the relevant outgoing eigenvector and a combined iterative method based on the characteristic bisection and Newton's method to solve the resulting system of equations for the mass parameters of the problem. Many such orbits are presented.

Key words: celestial mechanics - methods: numerical

1 Introduction

Asymptotic orbits are an interesting case of motion in dynamical systems, when the moving particle remains in the vicinity of unstable equilibrium points for an infinitely long time, and they can be used as reference orbits. Asymptotic orbits to equilibrium points were first studied by Str $\rm\ddot o$mgen (1935), and were associated with families of periodic orbits as their terminations with infinite period. These were spiralling asymptotics to the equilateral triangle equilibrium points of the restricted three-body problem. Str $\rm\ddot o$mgen's work was extended by Markellos et al. (1995), Perdios (1996, 1997), within the framework of the photogravitational restricted problem (see also Szebehely & Nacozy 1967; Danby 1967, 1984; Markellos 1991, among others, for the classical case). Non-spiralling asymptotic orbits to the collinear equilibrium points were studied by Deprit & Henrard (1965), Perdios & Markellos (1990), Perdios (1993) and Perdios et al. (2001). These orbits correspond to period discontinuities in the evolution of orbits within the families to which they belong. Asymptotic orbits to the collinear points and to the associated Lyapunov orbits around them were also studied in detail, both analytically and numerically, by Llibre et al. (1985) and Koon et al. (2000).

In this paper we compute non-spiralling doubly asymptotic orbits at the collinear equilibrium points of the general three-body problem. Analytical results for the existence of such orbits in the general three-body problem, such as those of Llibre et al. (1985) and Koon et al. (2000) in the restricted problem, are not available, to our knowledge. In the present paper we provide numerical evidence for their existence. We find several such orbits at each of the three collinear equilibrium points for both possible directions of the outgoing eigenvector. The determination of a doubly asymptotic orbit proceeds through a linear analysis to determine the relevant outgoing eigenvector, and a combination of the characteristic bisection and Newton's method to solve numerically the resulting system of equations for the mass parameters of the problem.

In the next section we give the equations of motion of the general three-body problem. In Sect. 3 we describe the procedure employed to determine doubly asymptotic orbits. Finally, in Sect. 4 we present the doubly asymptotic orbits computed.

2 Equations of motion

In a coordinate system ${\rm O}xy$, rotating with angular velocity $\dot\theta$, the equations of motion of three particles with masses m₁, m₂ and m₃, attracting each other according to the Newtonian law of gravitation, are expressed in the form (Hadjidemetriou & Christides 1975):

 

$$\displaystyle\ddot x - 2\dot\theta \dot y = bx+x\dot\theta^2+ \ddot\theta y + \mu ax_2,$$

$$\displaystyle\ddot x_2 = \left( m_3 b^\star + \dot\theta^2 \right)x_2 -\frac{(1-m_3)(1-\mu)^3}{x_2^2} \displaystyle + m_3(1-\mu)ax,$$

$$\displaystyle\ddot\theta = -\frac{2\dot\theta \dot x_2}{x_2} + \frac{m_3(1-\mu)ay}{x_2},$$ (1)

where:

$$\displaystyle a = \displaystyle-\left(\frac{1}{r_{13}^3}-\frac{1}{r_{23}^3}\right),$$

$$\displaystyle b = \displaystyle-\left(\frac{1-\mu}{r_{13}^3}+\fra... ...ce*{0.4cm} b^\star = -\left(\frac{\mu}{r_{13}^3}+\frac{1-\mu}{r_{23}^3}\right),$$

$$\displaystyle r_{13} = \displaystyle\left[\left(x+\frac{\mu}{1-\mu}x_2\right)^2+y^2 \right]^{1/2},$$

$$\displaystyle r_{23} = \displaystyle\left[\left(x-x_2\right)^2+y^2 \right]^{1/2},$$

$$\displaystyle \mu = \displaystyle\frac{m_2}{m_1+m_2}, \hspace*{0.5cm} m_1+m_2+m_3=1,$$ (2)

and we assume that $m_1\geqslant m_2\geqslant m_3$. The variables involved are the position of m₃(x,y), the distance x₂ of m₂ from the origin and the angle $\theta$which describes the orientation of the x-axis containing m₁ and m₂, with respect to an inertial frame. The position of m₁ on the x-axis is given by the equation:

$$x_1=-\frac{\mu}{1-\mu}x_2.$$ (3)

In the following we set:

$$\displaystyle X_1=x, \qquad X_2=y, \qquad X_3=x_2, \qquad X_4=\theta,$$

$$\displaystyle X_5=\dot x, \qquad X_6=\dot y, \qquad X_7= \dot x_2, \qquad X_8=\dot\theta.$$ (4)

3 Determination of doubly asymptotic orbits

To compute a doubly asymptotic orbit we first have to determine, for each pair of values of the mass parameters $(\mu ,m_3)$, the collinear equilibrium points L_i, i=1,2,3 on the rotating X₁-axis. To this end we set X_i(t)=0, i=2,5,6,7and X₈(t)=1 in the system of the equations of motion and the following system arises:

 

$$\displaystyle (1+b_{(0)})X_{1(0)}+\mu a_{(0)}X_{3(0)} = 0,$$

$$\displaystyle (1+m_3b_{(0)}^\star)X_{3(0)}-\frac{(1-m_3)(1-\mu)^3}{X_{3(0)}^2} + m_3(1-\mu) a_{(0)}X_{1(0)}=0,$$ (5)

where the subscript (0) throughout this paper denotes evaluation at the equilibrium configuration. By solving the above system numerically we can determine the exact location X₁₍₀₎ of m₃ and X₃₍₀₎ of m₂ on the X₁-axis for each one of the three collinear equilibrium configurations. We note here that "equilibrium configurations" is a proper term since each such configuration consists of three points on the X₁-axis. However, throughout this paper we use the term "equilibrium point", instead of equilibrium configuration, for simplicity. Due to the condition X₈(t)=1, these equilibrium points correspond to circular orbits of the three bodies in a fixed coordinate system.

For motion in the vicinity of an equilibrium point L we transfer the origin of the coordinate system by means of the translation:

$$X_1=X_{1(0)}+\xi_1,\quad X_5=\xi_5,$$

$$X_2=\xi_2, \qquad\qquad\! X_6=\xi_6,$$

$$X_3=X_{3(0)}+\xi_3,\quad\! X_7=\xi_7,$$

$$X_4=\xi_4, \qquad\qquad\!\! X_8=1+\xi_8,$$ (6)

and expand the equations of motion to first order terms in the variables. Then, the equations of motion are written in the form:

$$\dot\xi=A\xi,$$ (7)

where

$$\xi=\left(\xi_1,\ldots ,\xi_8\right)^\top \hspace*{0.5cm} {\r... ...\left( \begin{array}{cc} O & I \\ P & Q \end{array}\right),$$ (8)

with O, I being the zero and identity matrices respectively and

$$P=\left( \begin{array}{cccc} c_5^1 & 0 & c_5^3 & 0 \\ [1mm] ... ... 0 & 0 & c_7^8 \\ [1mm] 0 & 0 & c_8^7 & 0 \end{array}\right),$$ (9)

where we have abbreviated:

$$\displaystyle c_5^1 = 1+b_{(0)}+\left(\frac{\partial b}{\partial ... ...)_{(0)}X_{1(0)}+ \mu\left(\frac{\partial a}{\partial X_1}\right)_{(0)}X_{3(0)},$$

$$\displaystyle c_5^3 = \mu a_{(0)}+\left(\frac{\partial b}{\partia... ...)_{(0)}X_{1(0)}+\mu \left(\frac{\partial a}{\partial X_3}\right)_{(0)}X_{3(0)},$$

$$\displaystyle c_5^6= 2,\hspace*{0.5cm} c_5^8 =2X_{1(0)},$$

$$\displaystyle c_6^2= 1+b_{(0)}-\frac{m_3(1-\mu)a_{(0)}X_{1(0)}}{X_{3(0)}},$$

$$\displaystyle c_6^5= -2,\hspace*{0.5cm} c_6^7=\frac{2X_{1(0)}}{X_{3(0)}},$$

$$\displaystyle c_7^1= m_3\Bigg[(1-\mu)a_{(0)}+(1-\mu)\left(\frac{\... ...1(0)} +\left(\frac{\partial b^\star}{\partial X_1}\right)_{(0)}X_{3(0)} \Bigg],$$

$$\displaystyle c_7^3 = 1+m_3b_{(0)}^\star +\frac{2(1-m_3)(1-\mu)^3}{X_{3(0)}^3}+ m_3(1-\mu) \hspace*{0.3cm} \times~\left(\frac{\partial a}{\partial X_3}\righ... ...}X_{1(0)}+ m_3\left(\frac{\partial b^\star}{\partial X_3}\right)_{(0)}X_{3(0)},$$

$$\displaystyle c_7^8=2X_{3(0)},\hspace*{0.5cm} c_8^2 = \frac{m_3(1-\mu)a_{(0)}}{X_{3(0)}},\hspace*{0.5cm} c_8^7 =\frac{-2}{X_{3(0)}},$$ (10)

(see Markellos 1981). The characteristic equation of matrix A is then written in the form:

$$\lambda^2(\lambda^2+1)(\lambda^4+\beta\lambda^2+\gamma)=0,$$ (11)

where:

$$\displaystyle\beta=7-{\rm Tr}(P),$$

$$\displaystyle\gamma=9-3{\rm Tr}(P)+c_5^1 c_6^2 + c_5^1 c_7^3 + c_6^2 c_7^3 - c_5^3 c_7^1.$$ (12)

By solving Eq. (11) we determine the eigenvalues of A:

$$\displaystyle\lambda_1=\lambda_2=0, \hspace*{0.7cm} \lambda_{3,4}=\pm i,$$

$$\displaystyle\lambda_{5,6,7,8}=\pm\sqrt{-\frac{\beta}{2}\pm\sqrt{\left(\frac{\beta}{2}\right)^2-\gamma}}.$$ (13)

Important here is the outgoing eigenvector $\delta_5=\left(\delta_5^1,\delta_5^2,\ldots,\delta_5^8\right)^\top$corresponding to the largest real eigenvalue $\lambda_5$, which we normalize by putting $\delta_5^1=1$. The other components of the eigenvector $\delta_5$ are found to be:

 

$$\displaystyle \delta_5^2=\frac{\lambda_5}{c_8^2}\left( \lambda_5 \delta_5^4-c_8^7\delta_5^3\right),$$

$$\displaystyle \delta_5^3=\frac{2c_8^2c_7^8+c_7^1\left(c_6^2-\lamb... ...2c_5^8c_8^2+\left(c_6^2-\lambda_5^2\right) \left[4-(c_7^3-\lambda_5^2)\right]},$$

$$\displaystyle \delta_5^4=\frac{-c_7^1-\left(c_7^3-\lambda_5^2\right)\delta_5^3} {\lambda_5c_7^8},$$

$$\displaystyle \delta_5^{j+4}=\lambda_5 \delta_5^j, \hspace*{0.5cm} j=1,\ldots,4.$$ (14)

Thus, appropriate initial conditions of the asymptotic orbit are:

$$Y_0=X_0+\varepsilon \delta_5,$$ (15)

where X₀ is the equilibrium vector in phase space, namely:

$$X_0=\left(X_{1(0)},0,X_{3(0)},0,0,0,0,1\right)^\top,$$ (16)

and $\varepsilon$ is a small orbital parameter. With the vector Y₀ of initial conditions we proceed to integrate the equations of motion until the orbit intersects the X₁-axis for the $N{\rm th}$ time, $N=1,2,\ldots,$ and the following system must be satisfied in order to obtain a perpendicular intersection with the X₁-axis:

 

$$\displaystyle \dot x (\mu,m_3) = X_5(\mu,m_3) = 0,$$

$$\displaystyle \dot x_2 (\mu,m_3) = X_7(\mu,m_3) = 0.$$ (17)

This follows from the symmetry property of the equations of motion: $x \rightarrow x$, $y \rightarrow -y$, $x_2 \rightarrow x_2$, $\theta \rightarrow -\theta$, $t \rightarrow -t$(see Hadjidemetriou & Christides 1975).

To solve this system we employ a combination of the characteristic bisection and Newton's methods. The characteristic bisection method is based on the characteristic polyhedron concept. It amounts to constructing another refined characteristic polyhedron containing the solution, by "bisecting" a known one in order to calculate the solution with the desired accuracy (for details we refer to Vrahatis 1995). Although this method computes the solution of a system such as (17) with certainty, its convergence is slow compared to that of Newton's method, which is known to be quadratic. Newton's method is applied here as follows. Since System (17) is not satisfied initially we consider corrections $\delta\mu$ and $\delta m_{3}$ of $\mu$ and m₃such that:

 

$$\displaystyle X_5(\mu+\delta\mu,m_{3}+\delta m_{3}) = 0,$$

$$\displaystyle X_7(\mu+\delta\mu,m_{3}+\delta m_{3}) = 0.$$ (18)

We expand Eqs. (18) around $\mu$ and m₃ and by neglecting the nonlinear terms we obtain:

 

$$\displaystyle \frac{\partial X_5}{\partial \mu}\delta\mu +\frac{\partial X_5}{\partial m_{3}}\delta m_{3}= -X_5(\mu,m_{3}),$$

$$\displaystyle \frac{\partial X_7}{\partial \mu}\delta\mu +\frac{\partial X_7}{\partial m_{3}}\delta m_{3}= -X_7(\mu,m_{3}),$$ (19)

where the required derivatives wrt $\mu$ and m₃are computed approximately by means of additional integrations of the equations of motion. The appropriate corrections $\delta\mu$ and $\delta m_{3}$ for a step of Newton's method are then found by solving (19), and the procedure is repeated as required by our accuracy goals.

The procedure employed here is the following. First, the characteristic bisection method is applied to compute a doubly asymptotic orbit with a predetermined mild accuracy $\epsilon$, typically $\epsilon\leqslant 10^{-2}$. Then we utilize the previous estimate of the orbit as an initial approximation for Newton's method. In case a better approximation is not accomplished after one iteration of Newton's method, we reapply the characteristic bisection, starting from its previous stage, in order to achieve a superior estimate, say $\epsilon\leqslant 10^{-3}$. Then we apply again Newton's method and continue with the same procedure until the solution of (17) is obtained with the desired accuracy. This procedure ensures both certainty and speed of convergence (see for details Perdios et al. 2002).

4 Results

We have used the above analysis and numerical procedure, for a suitably small value of the orbital parameter $\varepsilon=\pm 0.00001$, to compute doubly asymptotic orbits to all three collinear equilibrium points for both directions of the outgoing eigenvector in each case, i.e. for both signs of the orbital parameter  $\varepsilon$. The results are presented in Table 1.

  

Table 1: Numerical data for doubly asymptotic orbits at L₁,L₂,L₃.

In this table we give the multiplicity N of the asymptotic orbit and the values of the parameters $\mu$ and m₃. We give the values for $\mu$ and m₃ to five significant figures which are verified to be correct by continuing the iteration until these five figures are stabilized. We have also verified by backward integration that for the value of $\varepsilon$ used the equilibrium point is recaptured with an accuracy of at least five significant figures. The values of the initial conditions are also given in Table 1 for the benefit of the reader although these are obtained from $\mu$ and m₃ by the analysis of Sect. 3. Note that in the table only the orbits with $N\leqslant 9$ are given although more were found corresponding to higher numbers of X₁-axis crossings. All kinds of doubly asymptotic orbits were found and represented in the table except at L₁ for  $\varepsilon >0$.

In Figs. 1-6 the asymptotic orbits of Table 1 marked by asterisks are shown.

	Figure 1: a) Doubly asymptotic orbit at L1 with $\varepsilon <0$, b) magnification of Fig. 1a.
Open with DEXTER

	Figure 2: a) Doubly asymptotic orbit at L1 with $\varepsilon <0$, b) magnification of Fig. 2a.
Open with DEXTER

	Figure 3: Doubly asymptotic orbit at L2 with $\varepsilon <0$.
Open with DEXTER

	Figure 4: Doubly asymptotic orbit at L2 with $\varepsilon >0$.
Open with DEXTER

	Figure 5: Doubly asymptotic orbit at L3 with $\varepsilon <0$.
Open with DEXTER

	Figure 6: Doubly asymptotic orbit at L3 with $\varepsilon >0$.
Open with DEXTER

In these figures the positions of m₁ and m₂ are shown for t=0. The outgoing part of the orbit is represented by the solid line and the incoming part by the dotted line.

We have also found some asymptotic orbits of the restricted three-body problem when within the framework of the general three body problem our numerical procedure converged to a pair of mass parameter values such that $\mu\neq 0$ but m₃=0. For example at L₁ for $\varepsilon >0$ and $\mu = 0.44359409$ we found an asymptotic orbit with initial conditions:

$${X_{10}=0.08020988, \hspace*{0.5cm} X_{20}=-0.00017779,}$$

$${X_{50}=0.00188872, \hspace*{0.5cm} X_{60}=-0.00067159,}$$

	Figure 7: Doubly asymptotic orbits of the restricted problem a) at L1 with $\varepsilon >0$, b) at L2 with $\varepsilon <0$, c) magnification of Fig. 7b.
Open with DEXTER

which has N=6 perpendicular crossings of the outgoing orbit with the X₁-axis. Similarly, we found an asymptotic orbit at L₂ for $\varepsilon <0$ and $\mu=0.01643677$ with initial conditions:

$${ X_{10}= 1.16972353, \hspace*{0.8cm} X_{20}=0.00032052,}$$

$${ X_{50}=-0.00106176, \hspace*{0.5cm} X_{60}=0.00068063,}$$

with N=9 perpendicular crossings of the outgoing orbit with the X₁-axis. For both the above orbits of the restricted problem we have $X_{30}=1-\mu$, X₄₀=X₇₀=0 and X₈₀=1. These orbits are shown in Fig. 7.

Acknowledgements

Thanks are due to Prof. V. V. Markellos for valuable discussions during this work. V. S. Kalantonis acknowledges financial support under a University of Patras "K. Karatheodory" research grant.

References

 

• Danby, J. M. A. 1967, AJ, 72, 198 In the text NASA ADS  
• Danby, J. M. A. 1984, Celes. Mech., 33, 251 In the text NASA ADS  
• Deprit, A., & Henrard, J. 1965, AJ, 70, 271 In the text NASA ADS  
• Hadjidemetriou, J. D., & Christides, Th. 1975, Celes. Mech., 12, 175 In the text NASA ADS  
• Koon, W. S., Lo, M. W., Marsden, J. E., & Ross, S. D. 2000, Chaos, 10(2), 427 In the text NASA ADS  
• Llibre, J., Martínez, R., & Simó, C. 1985, J. Diff. Eqs, 58, 104 In the text  
• Markellos, V. V. 1981, Celes. Mech., 25, 319 In the text NASA ADS  
• Markellos, V. V. 1991, in Predictability, stability and chaos in N-body dynamical systems, ed. A. E. Roy (Plenum press, New York), 413 In the text  
• Markellos, V. V., Perdios, E. A., & Papadakis, K. 1995, in From Newton to chaos, ed. A. E. Roy, & V. A. Steves (Plenum press, New York), 395 In the text  
• Perdios, E. A. 1993, Astrophys. Space Sci., 199, 185 In the text NASA ADS  
• Perdios, E. A. 1996, Astrophys. Space Sci., 240, 141 In the text NASA ADS  
• Perdios, E. A. 1997, Astrophys. Space Sci., 254, 61 In the text NASA ADS  
• Perdios, E. A., Kalantonis, V. S., & Vrahatis, M. N. 2002, Celes. Mech. Dynam. Astron., accepted  
• Perdios, E. A., & Markellos, V. V. 1990, Astrophys. Space Sci., 166, 129 In the text NASA ADS  
• Perdios, E. A., Ragos, O., Perdiou, A. E., & Vrahatis, M. N. 2001, Nonlinear Anal., 47, 3443 In the text  
• Str $\rm\ddot o$mgen, E. 1935, Copenhagen Obs. Publ., No. 100 In the text  
• Szebehely, V., & Nacozy, P. 1967, AJ, 72, 184 In the text NASA ADS  
• Vrahatis, M. N. 1995, J. Comp. Phys., 119, 105 In the text