Homoclinic dynamics in a spatial restricted four-body problem: blue skies into Smale horseshoes for vertical Lyapunov families

Abstract

The set of transverse homoclinic intersections for a saddle-focus equilibrium in the planar equilateral restricted four-body problem admits certain simple homoclinic orbits which form the skeleton of the complete homoclinic intersection—or homoclinic web. In the present work, the planar restricted four-body problem is viewed as an invariant subsystem of the spatial problem, and the influence of this planar homoclinic skeleton on the spatial dynamics is studied from a numerical point of view. Starting from the vertical Lyapunov families emanating from saddle-focus equilibria, we compute the stable/unstable manifolds of these spatial periodic orbits and look for intersections between these manifolds near the fundamental planar homoclinics. In this way, we are able to continue all of the basic planar homoclinic motions into the spatial problem as homoclinics for appropriate vertical Lyapunov orbits which, by the Smale tangle theorem, suggest the existence of chaotic motions in the spatial problem. While the saddle-focus equilibrium solutions in the planar problems occur only at a discrete set of energy levels, the cycle-to-cycle homoclinics in the spatial problem are robust with respect to small changes in energy.

Appendices

Appendix A: Obtaining a polynomial field by automatic differentiation of the CRFBP

To facilitate formal series calculations in the CRFBP, we first rewrite the problem as a first-order ordinary differential equation and then introduce a change of variable, often referred to as automatic differentiation, to obtain a polynomial vector field. The problem is recovered via projection, as long as the initial conditions are restricted to an appropriate submanifold. We first set

$$\begin{aligned} u_1 = x, \quad u_2 = \dot{x}, \quad u_3= y, \quad u_4 = \dot{y}, \quad u_5=z,\quad u_6= \dot{z}, \end{aligned}$$

(14)

and obtain a first-order ODE $\dot{u} = f(u)$ given by

$$\begin{aligned} \begin{aligned} \dot{u_1}&= u_2, \\ \dot{u_2}&= 2 u_4 +\varOmega _{u_1},\\ \dot{u_3}&= u_4, \\ \dot{u_4}&= -2u_2 +\varOmega _{u_3},\\ \dot{u_5}&= u_6, \\ \dot{u_6}&= \varOmega _{u_5}, \end{aligned} \end{aligned}$$

(15)

where $\varOmega$ is as previously given but using the new set of variable. This vector field still has singularities introduced by the terms corresponding to the inverse of the distance with the primaries, and we extend our set of variables using the following definitions:

$$\begin{aligned} u_7 = \frac{1}{\sqrt{(x-x_1)^2 +(y-y_1)^2 +(z-z_1)^2}} = \frac{1}{\sqrt{(u_1-x_1)^2 +(u_3-y_1)^2 +(u_5-z_1)^2}}, \end{aligned}$$

(16)

$$\begin{aligned} u_8 = \frac{1}{\sqrt{(x-x_2)^2 +(y-y_2)^2 +(z-z_2)^2}} = \frac{1}{\sqrt{(u_1-x_2)^2 +(u_3-y_2)^2 +(u_5-z_2)^2}}, \end{aligned}$$

(17)

$$\begin{aligned} u_9 = \frac{1}{\sqrt{(x-x_3)^2 +(y-y_3)^2 +(z-z_3)^2}} = \frac{1}{\sqrt{(u_1-x_3)^2 +(u_3-y_3)^2 +(u_5-z_3)^2}}. \end{aligned}$$

(18)

Let $U \subset {\mathbb {R}}^6$ be an open set excluding the primaries. Then, a direct computation provides that for the function $R:U \rightarrow {\mathbb {R}}^9$ given by

$$\begin{aligned} R(u_1,u_2,u_3,u_4,u_5,u_6) = \begin{pmatrix} u_1 \\ u_2 \\ u_3 \\ u_4 \\ u_5 \\ u_6 \\ \frac{1}{\sqrt{(u_1-x_1)^2 +(u_3-y_1)^2 +(u_5-z_1)^2}} \\ \frac{1}{\sqrt{(u_1-x_2)^2 +(u_3-y_2)^2 +(u_5-z_2)^2}} \\ \frac{1}{\sqrt{(u_1-x_3)^2 +(u_3-y_3)^2 +(u_5-z_3)^2}} \end{pmatrix} \end{aligned}$$

(19)

and the polynomial vector field $F :{\mathbb {R}}^9 \rightarrow {\mathbb {R}}^9$ given by

$$\begin{aligned} F(u)= \begin{pmatrix} u_2 \\ 2u_4 +u_1 +m_1(x_1-u_1)u_7u_7u_7 +m_2(x_2-u_1)u_8u_8u_8 +m_2(x_3-u_1)u_9u_9u_9 \\ u_4 \\ -2u_2 +u_3 +m_1(y_1-u_3)u_7u_7u_7 +m_2(y_2-u_3)u_8u_8u_8 +m_2(x_3-u_1)u_9u_9u_9 \\ u_6 \\ m_1(z_1-u_5)u_7u_7u_7 +m_2(z_2-u_5)u_8u_8u_8 +m_2(z_3-u_5)u_9u_9u_9 \\ (x_1-u_1)u_2u_7u_7u_7 +(y_1-u_3)u_4u_7u_7u_7 +(z_1-u_5)u_6u_7u_7u_7 \\ (x_2-u_1)u_2u_8u_8u_8 +(y_2-u_3)u_4u_8u_8u_8 +(z_2-u_5)u_6u_8u_8u_8 \\ (x_3-u_1)u_2u_9u_9u_9 +(y_3-u_3)u_4u_9u_9u_9 +(z_3-u_5)u_6u_9u_9u_9 \end{pmatrix}, \end{aligned}$$

(20)

we have the infinitesimal conjugacy

$$\begin{aligned} DR(u)f(u) = F(R(u)), \quad \forall u \in U. \end{aligned}$$

(21)

Hence, orbits of $u' = F(u)$ have the same dynamics as $x' = f(x)$ after projecting onto the first six components. We note that, as an effect of the change of variable, the new vector field does not have any singularity. Nevertheless, the dynamics of the two are related only on the graph of R, and R caries the singularities of f. The following items formalize the remarks just made.

• 1. Let $\pi :{\mathbb {R}}^9 \rightarrow {\mathbb {R}}^6$ denotes the projection onto the first six coordinates. So that for all $u \in U$ we have $u=\pi (R(u))$ and

  $$\begin{aligned} \pi (F(R(u))) =f(u). \end{aligned}$$

  Therefore, we recover the original problem.

• 2. The orbits of f are mapped onto orbits of F under R and the graph of R is invariant under the flow of F.

• 3. If $H:{\mathbb {R}}^6 \rightarrow {\mathbb {R}}$ is constant along curves solution of the initial system, then $G:{\mathbb {R}}^9 \rightarrow {\mathbb {R}}$ such that $G(R(u))= H(u)$ for all $u \in U$ is constant along curves solution of the extended problem.

It follows from those remarks that it is possible to find periodic orbits of the four-body problem using the vector field F. Our goal is to compute stable and unstable manifolds for a periodic orbit $\gamma (t)$ of ${\dot{u}}=f(u)$, so that we have to show that the associated periodic orbit $\varGamma (t)=R(\gamma (t))$ has the same stability type. This is the object of the next theorem.

Theorem 1

Let $\gamma (t)$ be a periodic orbit of ${\dot{u}}=f(u)$ with Floquet multiplier $\lambda$ associated with the tangent bundle v(t). Then $\lambda$ is a Floquet multiplier of the periodic orbit $\varGamma (t)=R(\gamma (t))$, solution to the system $\dot{x} = F(x)$, moreover $\xi (t)= DR(\gamma (t))v(t)$ is the associated tangent bundle.

Proof

We first note that v(t) will satisfy

$$\begin{aligned} \dot{v} (t)= Df(\gamma (t))v(t) -\lambda v(t) \end{aligned}$$

(22)

and that differentiating (21) provides

$$\begin{aligned} D^2R(u)f(u) +DR(u)Df(u) = DF(R(u))DR(u). \end{aligned}$$

(23)

So that a direct computation provides that

$$\begin{aligned} {\dot{\xi }}(t)&= D^2R(\gamma (t)){\dot{\gamma }}(t)v(t) +DR(\gamma (t))\dot{v}(t) \\&= D^2R(\gamma (t))f(\gamma (t))v(t) +DR(\gamma (t))Df(\gamma (t))v(t) -DR(\gamma (t))\lambda v(t) \end{aligned}$$

where we used the fact that $\gamma (t)$ is a periodic solution of f as well as Eq. (22). Then using (23), $\varGamma (t)=R(\gamma (t))$ and $\xi (t)= DR(\gamma (t))v(t)$, we obtain that

$$\begin{aligned} {\dot{\xi }}(t)= DF(\varGamma (t))\xi (t) -\lambda \xi (t). \end{aligned}$$

This is the desired result. $\square$

It follows from this result that in the extended system six of the multipliers will be known from the usual theory. The other three are all zeros so that the dimension of the stable and unstable manifolds for any orbits remains unchanged.

Appendix B: Orbit data

In this section, we provide several tables of data meant to make the present work more reproducible. Since our calculations of the connecting orbits utilize fairly sophisticated Fourier–Taylor approximations of the local stable/unstable manifolds in the formulation of the two-point boundary value problems, it is unreasonable to think that the casual reader would reimplement these calculations. On the other hand, many readers will have experience in the use of numerical integrators for problems in Celestial Mechanics, and once equipped with the equations of motion it is not unreasonable to think one might want to reproduce some of the periodic orbits and connections discussed in the present work. To this end, we provide accurate initial conditions which can be integrated to reproduce the orbits discussed in the present work. The resulting orbits could also be taken as initial conditions for numerical continuation software packages like AUTO or MatCont.

The table is organized in the same way. In the first column, we give the initial point expressed as a six-dimensional vector representing the initial position and momentum. The coordinates are given in the following order:

$$\begin{aligned} P_0 = ( x, \dot{x}, y, \dot{y}, z, \dot{z}). \end{aligned}$$

Then the second column of the table provides T an approximation of the period of the periodic orbit starting at the point previously given. The third column is n, the number of Floquet multipliers with positive real part. Finally, the last column shows $J(P_0)$, the energy level of the initial data. We note that that case of interest in this paper is when $n=2$ and the multipliers are complex conjugate. To obtain the data, we start by computing the center manifold of each libration point to find an initial guess for $P_0$ and T. To improve the guess, we numerically integrate the approximated periodic orbit and express the result in Fourier coefficients. Then Newton’s method is applied to obtain a guess for the periodic orbit with defect close to machine precision, for all cases covered by the tables it suffices to take 50 Fourier coefficients. The resulting sequence of Fourier coefficients is then a starting point for any continuation method in order to find other members of the family. To construct the table, we used a zeroth-order predictor–corrector algorithm using Newton’s method in the space of Fourier coefficients; in this case, the frequency is an unknown of the system, while the energy level is one of the inputs of the algorithm. The cases of ${\mathcal {L}}_0$ are given at $m_1=0.4$ and $m_2=0.35$, while the case at ${\mathcal {L}}_5$ is given with equal masses.

Table 1 Family at $L_0$

Table 2 Family at $L_5$, this table is computed with equal masses and we recall that periodic orbits at ${\mathcal {L}}_{4,6}$ can be obtained by a rotation of $\pm 120$ degrees. The cases of energy from 2.4 to 2.50 have real Floquet multipliers, while the remaining of the table are complex conjugate

The connecting orbits in Fig. 15 are homoclinic and accumulate to the periodic orbit with initial condition in the row $J=3.2$, given by the table for ${\mathcal {L}}_0$. The initial data with higher accuracy are

$$\begin{aligned} P_0= \begin{pmatrix} 0.134934339930888 \\ 0.003888013139251 \\ 0.117443350170703 \\ 0.000936082833871 \\ 0.216240831347475 \\ 0.101389225000425 \end{pmatrix}, \quad T= 2.998307362412966. \end{aligned}$$

To reproduce the trajectories displayed, one can integrate the following initial values $P_0$ back and forward in time for the given time T. The starting and ending points of the resulting trajectories will lay on the boundary of the parameterized unstable and stable manifolds, respectively.

$$\begin{aligned} P_0&= \begin{pmatrix} -\,0.585194841158983 \\ 0.650674788263036 \\ -\,0.242897059971999 \\ -\,0.809665850514842 \\ 0.015366927308435 \\ 0.645609803647810 \end{pmatrix}, \quad T= 3.9083,\\ P_0&= \begin{pmatrix} -\,0.010028232796882 \\ 0.018905042025788 \\ -\,0.527614375771166 \\ -\,0.278204402447460 \\ 0.066684862193223 \\ -\,0.421303149345866 \end{pmatrix}, \quad T= 3.5848,\\ P_0&= \begin{pmatrix} 0.364232983907004 \\ 0.282004601213298 \\ 0.365277731376544 \\ 0.566929250936993 \\ 0.188719573086921 \\ -\,0.192150484392393 \end{pmatrix}, \quad T= 4.1378. \end{aligned}$$

The three connecting orbits accumulating to the same periodic orbit and member of the families displayed in Fig. 19 can be found using the following initial condition and integration time.

$$\begin{aligned} P_0&= \begin{pmatrix} -\,0.101146445518484 \\ 0.039260485918423 \\ 0.357723970145646 \\ 0.064390124937530 \\ -\,0.215188925518734 \\ -\,0.117478748772784 \end{pmatrix}, \quad T= 2.3112, \\ P_0&= \begin{pmatrix} 0.292042336892103 \\ 0.003508935985276 \\ 0.118262267817677 \\ -\,0.031322268117327 \\ 0.009128811923180 \\ -\,0.481727032516309 \end{pmatrix}, \quad T= 1.7643, \\ P_0&= \begin{pmatrix} -\,0.082031603660355 \\ -\,0.244810917818636 \\ -\,0.371129071110934 \\ -\,0.141117360021255 \\ 0.090892927654736 \\ -\,0.410662336419204 \end{pmatrix}, \quad T= 2.6543. \end{aligned}$$

In the case of $L_5$, the connections computed is at $J= 2.9$, and the initial data for the periodic orbit are given with higher accuracy by

$$\begin{aligned} P_0= \begin{pmatrix} 0.919523300342616 \\ -\,0.018021865086785 \\ -\,0.005720721776858 \\ 0.001586045655911 \\ 0.140748196680255 \\ 0.142288965593486 \end{pmatrix}, \quad T= 5.485186773053060. \end{aligned}$$

The midpoint of each connecting orbit as well as the approximate integrating time needed to reach the boundary of the parameterized manifolds is given by pairs, corresponding to their shape and the figure in which they were presented. The initial data for the connecting orbit displayed in Fig. 16 are given by

$$\begin{aligned} P_0&= \begin{pmatrix} -\,0.221338679671589 \\ 0.290535064047762 \\ 0.807520893403199 \\ -\,0.079212158161279 \\ 0.099168243248453 \\ -\,0.192275633578254 \end{pmatrix}, \quad T= 3.9267,\\ P_0&= \begin{pmatrix} -\,0.027278885368683 \\ -\,0.415243675715196 \\ -\,0.681730123750280 \\ -\,0.689462085636593 \\ 0.002992463395721 \\ 0.060975647933203 \end{pmatrix}, \quad T= 4.1225. \end{aligned}$$

The initial data for the connecting orbit displayed in Fig. 17 are given by

$$\begin{aligned} P_0= \begin{pmatrix} -\,0.093216716467939 \\ 0.539629163160594 \\ 0.029112774416657 \\ 0.487738555586829 \\ 0.028678261003714 \\ -\,0.264834109382665 \end{pmatrix}, \quad T= 5.7103, \\ P_0= \begin{pmatrix} -\,0.369576889105909 \\ -\,0.085470029306448 \\ 0.463741869935280 \\ 0.545543943960734 \\ 0.088979849325104 \\ -\,0.103325348379878 \end{pmatrix}, \quad T= 4.7319. \end{aligned}$$

The initial data for the connecting orbit displayed in Fig. 18 are given by

$$\begin{aligned} P_0&= \begin{pmatrix} 0.497723454800157 \\ 0.803131159532739 \\ 1.346319122336476 \\ -\,0.067118235690442 \\ 0.029291317637547 \\ -\,0.145379858982222 \end{pmatrix}, \quad T= 4.9363, \\ P_0&= \begin{pmatrix} 0.480703865397053 \\ -\,0.766478203112893 \\ -\,1.325717528617470 \\ -\,0.049652453045623 \\ 0.255108298576897 \\ -\,0.003998765894041 \end{pmatrix}, \quad T= 4.9277. \end{aligned}$$