Hill four-body problem with oblate tertiary: an application to the Sun-Jupiter-Hektor-Skamandrios system

TL;DR

This paper analyzes a restricted four-body problem involving an oblate tertiary, deriving Hill's approximation, equilibrium points, and stability, with application to the Sun-Jupiter-Hektor-Skamandrios system.

Contribution

It introduces a Hill's approximation for a four-body system with an oblate tertiary, extending classical models to more realistic celestial configurations.

Findings

Triangular central configurations are isosceles due to oblateness.

Derived Hill's equations for the infinitesimal body near the oblate tertiary.

Identified equilibrium points and analyzed their stability.

Abstract

We consider a restricted four-body problem with a precise hierarchy between the bodies: two point-mass bigger bodies, a smaller one with oblate shape, and an infinitesimal body in the neighborhood of the oblate body. The three heavy bodies are assumed to move in a plane under their mutual gravity, and the fourth body moves under the gravitational influence of the three heavy bodies, but without affecting them. We start by finding the triangular central configurations of the three heavy bodies; since one body is oblate, the triangle is isosceles, rather than equilateral as in the point mass case. We assume that the three heavy bodies are in such a central configuration and we perform a Hill's approximation of the equations of motion describing the dynamics of the infinitesimal body in a neighborhood of the oblate body. Through the use of Hill's variables and a limiting procedure, this approximation amounts to sending the two other bodies to infinity. Finally, for the Hill approximation, we find the equilibrium points of the infinitesimal body and determine their stability. As a motivating example, we consider the dynamics of the moonlet Skamandrios of Jupiter's Trojan asteroid Hektor.

Figures

Equations

$$V(r,\phi,\lambda)={{\mathcal{G}m_{H}}\over r}\ \sum_{n=0}^{\infty}\Bigl{(}{R_{H}\over r}\Bigr{)}^{n}\ \sum_{m=0}^{n}P_{nm}(\sin\phi)\ (C_{nm}\cos m\lambda+S_{nm}\sin m\lambda)\ ,$$

$$\begin{array}{rl}{P}_n(x)& =\frac{1}{2^{n}n!}\frac{d^n}{d{x}^n}(x^2-1{)}^n\\ P_{nm}(x)& =(1-x^2{)}^{\frac{m}{2}}\frac{d^m}{d{x}^m}{P}_n(x),\end{array}$$

$$S_{n,m}$$

$$C_{2p+1,2q}$$

$$C_{2p,2q+1}$$

$$C_{2p,2q}$$

$$C_{20}$$

$$C_{22}$$

$$C_{20}=-0.476775$$

$$V(r,\phi,\lambda)={{\mathcal{G}m_{H}}\over r}\ \sum_{n=0}^{2}\Bigl{(}{R_{H}\over r}\Bigr{)}^{n}\ \sum_{m=0}^{n}P_{nm}(\sin\phi)\ C_{nm}\cos(m(\lambda+\Theta t)),$$

$$\begin{array}{rl}V(x,y,z)& =\frac{m_3}{r}-\frac{m_3}{r}{\left(\frac{R_3}{r}\right)}^2\left(\frac{J_2}{2}\right)\left(3{\left(\frac{z}{r}\right)}^2-1\right)\\ & =\frac{m_3}{r}+\frac{m_3}{r}{\left(\frac{R_3}{r}\right)}^2\left(\frac{C_{20}}{2}\right)\left(3\sin {\varphi }^2-1\right),\end{array}$$

$$V(q)=\frac{m_3}{r}+\frac{C{m}_3}{r^3},$$

$$C=R_3^2{J}_2/2>0.$$

$$\begin{array}{rl}{m}_1{\ddot{q}}_1& =m_1{m}_2(q_2-q_1)\frac{1}{\Vert q_2-q_1{\Vert }^3}+m_1{m}_3(q_3-q_1)\left[\frac{1}{\Vert q_3-q_1{\Vert }^3}+\frac{3C}{\Vert q_3-q_1{\Vert }^5}\right],\\ m_2{\ddot{q}}_2& =m_2{m}_1(q_1-q_2)\frac{1}{\Vert q_1-q_2{\Vert }^3}+m_2{m}_3(q_3-q_2)\left[\frac{1}{\Vert q_3-q_2{\Vert }^3}+\frac{3C}{\Vert q_3-q_2{\Vert }^5}\right],\\ m_3{\ddot{q}}_3& =m_3{m}_1(q_1-q_3)\frac{1}{\Vert q_1-q_3{\Vert }^3}+m_3{m}_2(q_2-q_3)\frac{1}{\Vert q_2-q_3{\Vert }^3},\end{array}$$

$$\mathbf{M}\ddot{\mathbf{q}}=\nabla U(\mathbf{q}),$$

$$U(\mathbf{q})=m_1{m}_2\frac{1}{r_{12}}+m_1{m}_3\left(\frac{1}{r_{13}}+\frac{C}{r_{13}^3}\right)+m_2{m}_3\left(\frac{1}{r_{23}}+\frac{C}{r_{23}^3}\right)$$

$$\mathbf{M}\mathbf{q}=\sum _{i=1}^3{m}_i{q}_i=0.$$

$$\left(\begin{array}[]{rr}\cos(\theta)&-\sin(\theta)\\ \sin(\theta)&\cos(\theta)\\ \end{array}\right)\in SO(2).$$

$$\mathbf{M}\left(\ddot{\mathbf{z}}+2\omega J\dot{\mathbf{z}}-{\omega }^2\mathbf{z}\right)=\nabla U,$$

$$\left(\begin{array}[]{rr}0&-1\\ 1&0\\ \end{array}\right).$$

$$\nabla U(\mathbf{z})+{\omega }^2\mathbf{M}\mathbf{z}=0.$$

$$I(\mathbf{z})=\sum _{1\le i<j\le 3}{m}_{ij}\Vert z_i-z_j{\Vert }^2:=\sum _{1\le i<j\le 3}{m}_{ij}{r}_{ij}^2.$$

$$\nabla f(\mathbf{z})=0,I(\mathbf{z})-\bar{I}=0,$$

$$\begin{array}{r}-\frac{1}{r_{12}^2}+{\omega }^2{r}_{12}=0,\\ -\frac{1}{r_{13}^2}-\frac{3C}{r_{13}^4}+{\omega }^2{r}_{13}=0,\\ -\frac{1}{r_{23}^2}-\frac{3C}{r_{23}^4}+{\omega }^2{r}_{23}=0,\\ m_1{m}_2{r}_{12}^2+m_1{m}_3{r}_{13}^2+m_2{m}_3{r}_{23}^2=\bar{I}.\end{array}$$

$$\omega =\sqrt{\frac{1}{r_{12}^3}}=\sqrt{\frac{1}{r_{13}^3}+\frac{3C}{r_{13}^5}}.$$

$$r_{12}:=v={\left(\frac{u^5}{u^2+3C}\right)}^{1/3}.$$

$$m_1{m}_2{\left(\frac{u^5}{u^2+3C}\right)}^{2/3}+(m_1{m}_3+m_2{m}_3)u^2=\bar{I},$$

$$\frac{m_1^3{m}_2^3{u}^{10}}{{\left(u^2+3C\right)}^2}={\left(\bar{I}-(m_1{m}_3+m_2{m}_3)u^2\right)}^3.$$

$$\frac{a{z}^5}{(z+b{)}^2}=(\bar{I}-cz{)}^3.$$

$$k^{\prime }(z)=\frac{3a{z}^6+8ab{z}^5+5a{b}^2{z}^4}{(z+b{)}^4}>0$$

$$l^{\prime }(z)=-3c(\bar{I}-cz{)}^2<0.$$

$$k(0)=0\text{ and }\underset{z\to +\infty }{\lim }k(z)=+\infty ,$$

Full text

We start by finding the triangular central configurations of the three heavy bodies; since one body is oblate, the triangle is isosceles, rather than equilateral as in the point mass case. We assume that the three heavy bodies are in such a central configuration and we perform a Hill’s approximation of the equations of motion describing the dynamics of the infinitesimal body in a neighborhood of the oblate body. Through the use of Hill’s variables and a limiting procedure, this approximation amounts to sending the two other bodies to infinity. Finally, for the Hill approximation, we find the equilibrium points of the infinitesimal body and determine their stability. As a motivating example, we consider the dynamics of the moonlet Skamandrios of Jupiter’s Trojan asteroid Hektor.

1. Introduction

The discovery of binary asteroids has led to considering dynamical models formed by four bodies, two of them being typically the Sun and Jupiter. Among possible four-body models (see also [HS86, Sch98, GJ03, SB05, ARV09, BGD13b, BGD13a, BG16, KMJ17]), a relevant role is played by the models in which three bodies lie on a triangular central configuration. Given that asteroids have often a (very) irregular shape, it is useful to start by investigating the case in which one body has an oblate shape. Among the different questions that this model may rise, we concentrate on the existence of equilibrium points and the corresponding linear stability analysis. Within such framework, we consider a four-body simplified model and we concentrate on the specific example given by the Trojan asteroid 624 Hektor, which is located close to the Lagrangian point $L_{4}$ of the Sun-Jupiter system, and its small moonlet.

In our model Sun, Jupiter and Hektor form an isosceles triangle (nearly equilateral) whose shape remains unchanged over time. The small body represents Hektor’s moonlet Skamandrios. Obviously, we could also replace the moonlet by a spacecraft orbiting Hektor. The system Sun-Jupiter-Hektor-Skamandrios plays a relevant role for several reasons. Indeed, Hektor is the largest Jupiter Trojan, it has one of the most elongated shapes among the bodies of its size in the Solar system, and it is the only known Trojan to possess a moonlet (see, e.g., [DLZ12] for stability regions of Trojans around the Earth and [LC15] for dissipative effects around the triangular Lagrangian points).

The study of asteroids with satellites presents a special interest for planetary dynamics, as they provide information about constraints on the formation and evolution of the Solar system.

Another motivation to study the dynamics of a small body near a Trojan asteroid comes from astrodynamics, as NASA prepares the first mission, Lucy, to the Jupiter’s Trojans, which is planned to be launched in October 2021 and visit seven different asteroids: a Main Belt asteroid and at least five Trojans.

As a model for the Sun-Jupiter-Hektor dynamics, we consider a system of three bodies of masses $m_{1}\geq m_{2}\geq m_{3}$, which move in circular orbits under mutual gravity, and form a triangular central configuration. We refer to these bodies as the primary, the secondary, and the tertiary, respectively. We assume that the first two bodies of masses $m_{1},m_{2}$ are spherical and homogeneous, so they can be treated as point masses, while the third body of mass $m_{3}$ is oblate. We describe the gravitational potential of $m_{3}$ in terms of spherical harmonics, and we only retain the most significant ones. We show the existence of a corresponding triangular central configuration, which turns out to be an isosceles triangle; if the oblateness of the mass $m_{3}$ is made to be zero, the central configuration becomes the well-known equilateral triangle Lagrangian central configuration. We stress that when $m_{3}$ is oblate, the central configuration is not the same as in the non-oblate case, since the overall gravitational field of $m_{3}$ is no longer Newtonian; it is well known that central configurations depend on the nature of the gravitational field (see, e.g., [CLPC04, APC13, DSZ, MS17]). We note that there exist papers in the literature (e.g., [APHS16]), which consider systems of three bodies, with one of the bodies non-spherical, which are assumed to form an equilateral triangle central configuration. Such assumption, while it may lead to very good approximations, is not physically correct.

The moonlet Skamandrios is represented by a fourth body, of infinitesimal mass, which moves in a vicinity of $m_{3}$ under the gravitational influence of $m_{1},m_{2},m_{3}$, but without affecting their motion. We consider the motion of the infinitesimal mass taking place in the three-dimensional space; it is not restricted to the plane of motion of the three heavy bodies. This situation is referred to as the spatial circular restricted four-body problem, and can be described by an autonomous Hamiltonian system of $3$-degrees of freedom.

We ‘zoom-in’ onto the dynamics in a small neighborhood of $m_{3}$ by performing a Hill’s approximation of the restricted four-body problem. This is done by a rescaling of the coordinates in terms of $m_{3}^{1/3}$, writing the associated Hamiltonian in the rescaled coordinates as a power series in $m_{3}^{1/3}$, and neglecting all the terms of order $O(m_{3}^{1/3})$ in the expansion, since such terms are small when $m_{3}$ is small. This yields an approximation of the motion of the massless particle in an $O(m_{3}^{1/3})$-neighborhood of $m_{3}$, while $m_{1}$ and $m_{2}$ are ‘sent to infinity’ through the rescaling. This model is an extension of the classical lunar Hill problem [Hil78]. Since the tertiary is assumed to be oblate, and the corresponding central configuration formed by the three heavy bodies is not an equilateral triangle anymore, this model also extends Hill’s approximation of the restricted four-body problem developed in [BGG15].

The Hill approximation is more advantageous to utilize for this system than the restricted four-body problem, since it allows for an analytical treatment, and yields more accurate numerical implementations when realistic parameters are used. The main numerical difficulty in the restricted four-body problem is the large differences of scales among the relevant parameters, i.e. the mass of Hektor is much smaller than the masses of the other two heavy bodies. The rescaling of the coordinates involved in the Hill approximation reduces the difference of scales of the parameters to more manageable quantities; more precisely, in normalized units the oblateness effect in the restricted four-body problem is of the order $O(10^{-15})$, while in the Hill approximation is of the order $O(10^{-7})$ (see Section 3.2 for details).

Once we have established the model for the Hill four-body problem with oblate tertiary, we study the equilibrium points and their linear stability. We find that there are $2$ pairs of symmetric equilibrium points on each of the $x$-, $y$-, and $z$-coordinate axes, respectively. The equilibrium points on the $x$- and $y$-coordinate axes are just a continuation of the corresponding ones for the Hill four-body problem with non-oblate tertiary [BGG15]. The equilibrium points on the $z$-coordinate axis constitute a new feature of the model. In the case of Hektor, these equilibrium points turn out to be outside of the body of the asteroid but very close to the surface, so they are of potential interest for low altitude orbit space missions, such as the one of NASA/JPL’s Dawn mission around Vesta ([Del11]).

This work is organized as follows. In Section 2 we describe in full details the restricted four-body model in which the tertiary is oblate; in particular, we describe the isosceles triangle central configuration of three bodies in which two bodies are point masses and the third is oblate. Hill’s approximation is introduced in Section 3. The determination of the equilibria and their stability is given in Section 4.

2. Restricted four-body problem with oblate tertiary

In this section we develop a model for a restricted four-body problem, which consists of two bigger bodies (e.g., the Sun and Jupiter), a smaller body – called tertiary – with oblate shape (e.g. an asteroid), and an infinitesimal mass (e.g., moonlet) around the tertiary.

As mentioned in Section 1, we consider the three masses $m_{1}\geq m_{2}\geq m_{3}$ as moving under the mutual gravitational attraction; the bodies with masses $m_{1}$ and $m_{2}$ are considered as point masses, while $m_{3}$ is the oblate body. We normalize the units of mass so that $m_{1}+m_{2}+m_{3}=1$.

We assume that the bodies with masses $m_{1}$, $m_{2}$, $m_{3}$ move on a triangular central configuration, which will be determined in Section 2.3, once the gravitational field of the oblate body has been discussed in Section 2.2. We will concentrate on the specific example given by the asteroid Hektor and its moon Skamandrios, where Hektor moves on a central configuration with Jupiter and the Sun. Orbital and physical values are given in Section 2.1. The positions of the three man bodies in the triangular central configuration is computed in Section 2.4, while the equations of motion of the moonlet – with infinitesimal mass moving in the vicinity of $m_{3}$ – are given in Section 2.5.

2.1. Data on the Sun-Jupiter-Hektor-Skamandrios system

The models which we will develop below will be applied to the case of the Sun-Jupiter-Hektor-Skamandrios system. We extract the data for this system from [JPL, MDCR*+*14, Des15].

Hektor is approximately located at the Lagrangian point $L_{4}$ of the Sun-Jupiter system. According to [Des15], Hektor is approximately $416\times 131\times 120$ km in size, and its shape can be approximated by a dumb-bell figure; the equivalent radius (i.e., the radius of a sphere with the same volume as the asteroid) is $R_{H}=92$ km111Note that [Des15] claims that there are some typos in the values reported in [MDCR*+*14]..

Hektor spins very fast, with a rotation period of approximately $6.92$ hours (see the JPL Solar System Dynamics archive [JPL]).

The moonlet Skamandrios orbits around Hektor at a distance of approximately $957.5$ km, with an orbital period of $2.965079$ days; see [Des15]. Its orbit is highly inclined, at approximately $50.1^{\circ}$ with respect to the orbit of Hektor, which justifies choosing as a model the spatial restricted four-body problem rather than the planar one; see [MDCR*+*14].

We also note that the inclination of Hektor is approximately $18.17^{\circ}$ (see [JPL]). Although a more refined model should include a non-zero inclination, we will consider that Sun-Jupiter-Hektor move in the same plane, an assumption that is needed in order for the three bodies to form a central configuration. We will further assume that the axis of rotation of Hektor is perpendicular to the plane of motion.

For the masses of Sun, Jupiter and Hektor we use the values of $m_{1}=1.989\times 10^{30}$ kg, $m_{2}=1.898\times 10^{27}$, and $m_{3}=7.91\times 10^{18}$ kg, respectively. For the average distance Sun-Jupiter we use the value $778.5\times 10^{6}$ km.

In Figure 1 we provide a comparison between the strength of the different forces acting on the moonlet: the Newtonian gravitational attraction of Hektor, Sun, Jupiter, and the effect of the non-spherical shape of the asteroid, limited to the the so-called $J_{2}$ coefficient, which will be introduced in Section 2.2.

2.2. The gravitational field of a non-spherical body

We first consider that the tertiary body, representing Hektor, has a general (non-spherical) shape. The gravitational potential, relative to a reference frame centered at the barycenter of the tertiary and rotating with the body, is given in spherical coordinates $(r,\phi,\lambda)$ by (see, e.g., [CG18]):

[TABLE]

where $\mathcal{G}$ is the gravitational constant, $m_{H}$ is the mass of Hektor, $R_{H}$ is its average radius, $P_{nm}$ are the Legendre polynomials defined as

[TABLE]

and $C_{nm}$ and $S_{nm}$ are the spherical harmonics coefficients.

In the case of an ellipsoid of semi-axes $a\geq b\geq c$, we have the following explicit formulas ([Boy97]):

[TABLE]

In particular, $C_{20}$ and $C_{22}$ turn out to be given by simple expressions

[TABLE]

For the Sun-Jupiter-Hektor data, we take $a=208$ km, $b=65.5$ km, $c=60$ km, $R_{H}=92$ km, following [Des15] (see Section 2.1), we calculate the following coefficients:

[TABLE]

Notice that each term $C_{2p,2q}$ is multiplied in (2.2) by the factor $R_{H}^{2p}/r^{2p+1}$. For $r$ equal to the average distance from the moonlet to the asteroid, we have $R_{H}/r\approx 0.096$. Therefore, in the following we will ignore the effect of the coefficients $C_{2p,2q}$, with $p\geq 2$, which are at least of $O(R_{H}^{4}/r^{5})$.

Note that the value of $C_{20}$ computed above is significantly bigger in absolute value than the one reported in [MDCR*+*14], which equals $-0.15$. The reason is that we use different estimates for the size of Hektor, following [Des15] (see Section 2.1).

If we consider a frame centered at the barycenter of the tertiary, and which rotates with the angular velocity of the tertiary about the primary, the time dependent gravitational potential takes of the form

[TABLE]

where $\Theta$ represents the frequency of the spin of Hektor.

For $n=2$, $m=0$ the corresponding term $C_{nm}\cos(m(\lambda+\Theta t))$ in the summation (2.2) is equal to $C_{20}$ and is independent of time; for $n=2$, $m=2$ the corresponding term $C_{nm}\cos(m(\lambda+\Theta t))$ is equal to $C_{22}\cos(2(\lambda+\Theta t))$, so is time-dependent. We do not consider the other terms in the sum (2.2).

Since the ratio of the rotation period of Hektor to the orbital period of the moonlet is relatively small, approximately $0.09740991$, in this paper we will only consider the average effect of $C_{22}\cos(2(\lambda+\Theta t))$ on the moonlet, which is zero.

In conclusion, in the model below we will only consider the effect of $C_{20}:=-J_{2}<0$, which amounts to approximating Hektor as an oblate body (i.e., an ellipsoid of revolution obtained by rotating an ellipse about its minor axis); the dimensionless quantity $J_{2}$ is referred to as the zonal harmonic in the gravitational potential. The term corresponding to $C_{20}$ is the larger one, followed by that corresponding to $C_{22}$; however, since the $C_{22}$ term introduces a time dependence, thus further complicating the model, we start by disregarding it and plan to study its effect in a future work.

2.3. Central configurations for the three-body problem with one oblate body

We now consider only the three heavy bodies, of masses $m_{1}\geq m_{2}\geq m_{3}$, with the body of mass $m_{3}$ being oblate, in which case we only take into account the term corresponding to $C_{20}=-J_{2}$ in (2.2). We write the approximation of the gravitational potential of the tertiary (2.2) in both Cartesian and spherical coordinates (in the frame of the tertiary and rotating with the body):

[TABLE]

where $m_{3}$ is the normalized mass of Hektor (the sum of the three masses is the unit of mass), $R_{3}:=R_{H}$ is the average radius of Hektor in normalized units (the distance between Sun and Hektor is the unit of distance), the gravitational constant is normalized to $1$, and $\sin\phi=z/r$.

We want to find the triangular central configurations formed by $m_{1}$, $m_{2}$, $m_{3}$; we will follow the approach in [APC13]. Since for a central configuration the three bodies lie in the same plane, in the gravitational field (2.3) of $m_{3}$ we set $\phi=0$, obtaining

[TABLE]

where $q=(x,y)$ is the position vector of an arbitrary point in the plane, $r=\|q\|$ is the distance from $m_{3}$, and we denote

[TABLE]

Let $q_{i}$ be the position vector of the mass $m_{i}$, for $i=1,2,3$, in an inertial frame centered at the barycenter of the three bodies.

The equations of motion of the three bodies are

[TABLE]

where the last terms in the first two equations are due to (2.4), and the gravitational constant is normalized to $\mathcal{G}=1$. Denote $r_{ij}=\|q_{i}-q_{j}\|$, for $i\neq j$, ${\bf q}=(q_{1},q_{2},q_{3})$, and ${\bf M}=\textrm{diag}(m_{1},m_{1},m_{2},m_{2},m_{3},m_{3})$ the $6\times 6$ matrix with $2$ copies of each mass along the diagonal. Then (2.6) can be written as

[TABLE]

where

[TABLE]

is the potential for the three body problem with oblate $m_{3}$.

Let us assume that the center of mass is fixed at the origin, i.e.,

[TABLE]

We are interested in relative equilibrium solutions for the motion of the three bodies, which are characterized by the fact they become equilibrium points in a uniformly rotation frame.

Denote by $R(\theta)$ the $6\times 6$ block diagonal matrix consisting of $3$ diagonal blocks the form

[TABLE]

Substituting ${\bf q}(t)=R(\omega t){\bf z}(t)$ for some $\omega\in\mathbb{R}$ in (2.7), where ${\bf z}=(z_{1},z_{2},z_{3})\in\mathbb{R}^{6}$, we obtain

[TABLE]

where $J$ is the block diagonal matrix consisting of $3$ diagonal blocks the form

[TABLE]

The condition for an equilibrium point of (2.3) yields the algebraic equation

[TABLE]

A solution ${\bf z}$ of the three-body problem satisfying (2.11) is referred to as a central configuration. This is equivalent to $\ddot{z}_{i}=-\omega^{2}z_{i}$, for $i=1,2,3$, meaning that the accelerations of the masses are proportional to the corresponding position vectors, and all accelerations are pointing towards the center of mass. Thus, the solution ${\bf q}(t)$ is a relative equilibrium solution if and only if ${\bf q}(t)=R(\omega t){\bf z}(t)$ with ${\bf z}(t)$ being a central configuration solution, and the rotation $R(\omega t)$ being a circular solution of the Kepler problem.

Let $I({\bf z})={\bf z}^{T}M{\bf z}=\sum_{i}m_{i}\|z_{i}\|^{2}$ be the moment of inertia. It is easy to see that this is a conserved quantity for the motion, that is, $I({\bf z}(t))=\bar{I}$ for some $\bar{I}$ at all $t$. Using Lagrange’s second identity (see, e.g., [GN12]), and that ${\bf M}{\bf z}=0$, normalizing the masses so that $\sum_{i=1}^{3}m_{i}=1$, the moment of inertia can be written as:

[TABLE]

Thus, central configurations correspond to critical points of the potential $U$ on the sphere ${\bf z}^{T}M{\bf z}=1$, which can be obtained by solving the Lagrange multiplier problem

[TABLE]

where $f({\bf z})=U({\bf z})+\frac{1}{2}\omega^{2}(I({\bf z})-\bar{I})$. In the above, we used the fact that $\nabla I({\bf z})=2{\bf M}{\bf z}$.

We solve this problem in the variables $r_{ij}=\|z_{i}-z_{j}\|$ for $1\leq i<j\leq 3$, since both $U$ and $I$ can be written in terms of these variables. This reduces the dimension of the system (2.13) from $7$ equations to $4$ equations. Denote ${\bf r}=(r_{12},r_{13},r_{23})$, and let $\tilde{f}({\bf r})$ be the function $f$ expressed in the variable ${\bf r}$, that is $\tilde{f}({\bf r}({\bf z}))=f({\bf z})$. By the chain rule, $\nabla_{r}\tilde{f}\cdot\left(\frac{\partial{\bf r}}{\partial{\bf z}}\right)=\nabla_{\bf z}f({\bf z})$. It is easy to see that the rank of the matrix $\left(\frac{\partial{\bf r}}{\partial{\bf z}}\right)$ is maximal provided that $z_{1},z_{2},z_{3}$ are not collinear (for details, see [CLPC04, APC13]). As we are looking for triangular central configurations, this condition is satisfied. Thus, $\nabla_{r}\tilde{f}({\bf r})=0$ if and only if $\nabla_{\bf z}f({\bf z})=0$. In other words, we can now solve the system (2.13) in the variable ${\bf r}$. We obtain

[TABLE]

Note that the function $h(r)=-\frac{1}{r^{3}}-\frac{3C}{r^{5}}$ has positive derivative $h^{\prime}(r)=\frac{3}{r^{4}}+\frac{15C}{r^{6}}$ with $C,r>0$, hence it is injective. Thus, the second and third equation in (2.14) yield $r_{13}=r_{23}:=u$. Solving for $\omega$ in the first and second equation we obtain:

[TABLE]

Solving for $r_{12}$ yields:

[TABLE]

Notice that $C>0$ implies $0<v<u$.

The condition $I({\bf r})=\bar{I}$ yields

[TABLE]

which is equivalent to

[TABLE]

To simplify the notation, let $u^{2}=z$, $m_{1}^{3}m_{2}^{3}=a$, $m_{1}m_{3}+m_{2}m_{3}=c$, $3C=b$, obtaining

[TABLE]

The function $k(z)=\displaystyle\frac{az^{5}}{(z+b)^{2}}$ has derivative

[TABLE]

and the function $l(z)=(\bar{I}-cz)^{3}$ has derivative

[TABLE]

Since

[TABLE]

and

[TABLE]

the equation $k(z)=l(z)$ has a unique solution with $z>0$. We conclude that for each fixed $\bar{I}$, there is a unique solution of (2.14). Thus, we have proved the following result.

Proposition 2.1**.**

In the three-body problem with one oblate primary, for every fixed value $\bar{I}$ of the moment of inertia there exists a unique central configuration, which is an isosceles triangle.

We note that, while [APC13] studies central configurations of three oblate bodies (as well as of three bodies under Schwarzschild metric), the isosceles central configuration found above is not explicitly shown there (see Theorem 4 in [APC13]).

To put this in quantitative perspective, when we use the data from Section 2.1 in (2.5), we obtain $C=3.329215\times 10^{-15}$. If we set $u=r_{13}=r_{23}=1$, from (2.16) we obtain $v=r_{12}=0.9999999999999967=1.0-3.3\times 10^{-15}$. In terms of the Sun-Jupiter distance $r_{12}=778.5\times 10^{6}$ km, the distance $r_{13}=r_{23}$ differs from the corresponding distance in the equilateral central configuration by $2.7\times 10^{-6}$ km. Practically, this isosceles triangle central configuration is almost an equilateral triangle.

2.4. Location of the bodies in the triangular central configuration

We now compute the expression of the location of the three bodies in the triangular central configuration, relative to a synodic frame that rotates together with the bodies, with the center of mass fixed at the origin, and the location of $m_{1}$ on the negative $x$-semi-axis. We assume that the masses lie in the $z=0$ plane. Instead of fixing the value $\bar{I}$ of the moment of inertia, we fix $u=r_{13}=r_{23}=1$ and have $v=r_{12}<1$ given by (2.16). Then, we obtain the following result.

Proposition 2.2**.**

In the synodic reference frame, the coordinates of the three bodies in the triangular central configuration, satisfying the constraints

[TABLE]

are given by

[TABLE]

Proof.

Denote by $A=x_{1}-x_{2}$ and $B=x_{1}-x_{3}$, so $x_{2}=x_{1}-A$, $x_{3}=x_{1}-B$, and $x_{3}-x_{2}=A-B$. Substituting these in (2.20) we obtain $(m_{1}+m_{2}+m_{3})x_{1}=m_{2}A+m_{3}B$. From (2.22) it follows $x_{1}=m_{2}A+m_{3}B=m_{3}({\bar{\mu}}A+B)$, where we denoted ${\bar{\mu}}:=m_{2}/m_{3}$. From (2.23) and (2.21) we have $y_{3}=-(m_{2}/m_{3})y_{2}=-{\bar{\mu}}y_{2}$, and $y_{3}-y_{2}=-(1+{\bar{\mu}})y_{2}$. From (2.17), we can solve for $y_{2}$ in terms of $A$ (see (2.25) below). From now on, the objective is to solve for $A$ and $B$, which in turn will yield $x_{1},x_{2},x_{3},y_{2},y_{3}$.

Equations (2.17), (2.18), (2.19) become:

[TABLE]

Adding (2.25) and (2.26), and subtracting (2.27) yields

[TABLE]

From (2.25), $y_{2}^{2}=v^{2}-A^{2}$, so (2.28) yields

[TABLE]

Substituting in (2.26) and solving for $A$ yields

[TABLE]

Substituting $y_{2}^{2}=v^{2}-A^{2}$ in (2.28) and solving for $B$ yields

[TABLE]

with $\textrm{sign}(A)=\textrm{sign}(B)$.

Substituting $A$, $B$ and ${\bar{\mu}}$ in $x_{1}=m_{2}A+m_{3}B$, and choosing the negative sign for $x_{1}$ to agree with our initial choice that $x_{1}<0$, after simplification, we first obtain the value of $x_{1}$ below. Then, substituting in $x_{2}=x_{1}-A$, $x_{3}=x_{1}-B$, $y_{2}^{2}=v^{2}-A^{2}$, $y_{3}=-(m_{2}/m_{3})y_{2}$, we compute $x_{2}$, $x_{3}$, $y_{2}$, $y_{3}$, obtaining (2.24).

For future reference, we note that if we let $m_{3}\to 0$ in (2.24), we obtain

[TABLE]

∎

Remark 2.3*.*

In the case when the oblateness coefficient $J_{2}$ of $m_{3}$ is made equal to zero, then $v=1$, and in (2.24) we obtain the Lagrangian equilateral triangle central configuration, with the position given by the following equivalent formulas (see, e.g., [BP13]):

[TABLE]

where $K=m_{2}(m_{3}-m_{2})+m_{1}(m_{2}+2m_{3})$.

Notice that the equations (2.33) are expressed in terms of $m_{1},m_{2},m_{3}$, while (2.24) are expressed in terms of $m_{2},m_{3}$; we obtain corresponding expressions that are equivalent when we substitute $m_{1}=1-m_{2}-m_{3}$ in (2.33). One minor difference is that in (2.33) the position of $x_{1}$ is not constrained to be on the negative $x$-semi-axis, as we assumed for (2.24); the position of $x_{1}$ in (2.33) depends on the quantity $\textrm{sign}(K)$; when $\textrm{sign}(K)>0$ we have $|K|/K=1$, and the equations (2.24) become equivalent with the equation (2.33).

We remark that when $m_{3}\to 0$, the limiting position of the three masses in (2.33) is given by:

[TABLE]

with $(x_{1},y_{1})$ and $(x_{2},y_{2})$ representing the position of the masses $m_{1}$ and $m_{2}$, respectively, and $(x_{3},y_{3})$ representing the position of the equilibrium point $L_{4}$ in the planar circular restricted three-body problem.

2.5. Equations of motion for the restricted four-body problem with oblate tertiary

Now we consider the dynamics of a fourth body in the neighborhood of the tertiary. This fourth body represents the moonlet Skamandrios orbiting around Hektor. We model the dynamics of the fourth body by the spatial, circular, restricted four-body problem, meaning that the moonlet is moving under the gravitational attraction of Hektor, Jupiter and the Sun, without affecting their motion which remains on circular orbits and forming a triangular central configuration as in Section 2.3. As before, we assume that Hektor has an oblate shape, with the gravitational potential given by (2.3).

The equations of motion of the infinitesimal body relative to a synodic frame of reference that rotates together with the three bodies is given by

[TABLE]

where the effective potential ${\tilde{\Omega}}$ is given by

[TABLE]

with $(x_{i},y_{i},z_{i})$ representing the $(x,y,z)$-coordinates in the synodic reference frame of the body of mass $m_{i}$, $r_{i}=\left((x-x_{i})^{2}+(y-y_{i})^{2}+z^{2}\right)^{\frac{1}{2}}$ is the distance from the moonlet to the mass $m_{i}$, for $i=1,2,3$, and $\omega$ is the angular velocity of the system of three bodies around the center of mass given by (2.15). The perturbed mean motion of the primaries $\omega$ in the above equation depends on the oblateness parameter. The coordinates $(x_{1},y_{1})$, $(x_{2},y_{2})$, $(x_{3},y_{3})$ of the bodies $m_{1}$, $m_{2}$, $m_{3}$, respectively, are given by (2.24), while $z_{i}=0$, for $i=1,2,3$.

For the choice that we have made $r_{13}=r_{23}=u=1$ and $r_{12}=v$ satisfying (2.16), we have

[TABLE]

We remark that if we set $m_{2}=0$ we obtain the restricted three-body problem with one oblate body and $\omega=\sqrt{1+3R_{H}^{2}J_{2}/2}$ agrees with the formulas in [McC63, SR76, AGST12]. If $m_{3}$ has no oblateness, i.e., $J_{2}=C_{20}=0$, then $\omega=1$.

We rescale the time so that in the new units the angular velocity is normalized to $1$, obtaining

[TABLE]

with

[TABLE]

The equations of motion (2.37) have the total energy $H$ defined below as a conserved quantity:

[TABLE]

We switch to the Hamiltonian setting by considering the system of symplectic coordinates $(x,y,z,p_{x},p_{y},p_{z})$ with respect to the symplectic form $\varpi=x\wedge p_{x}+y\wedge p_{y}+z\wedge p_{z}$, and making the transformation $\dot{x}=p_{x}+y$, $\dot{y}=p_{y}-x$ and $\dot{z}=p_{z}$. We obtain:

[TABLE]

where we denote $C^{\prime}={R_{3}}^{2}C_{20}/2$. Thus, the equations of motion (2.35) are equivalent to the Hamilton equations for the Hamiltonian given by (2.38).

3. Hill four-body problem with oblate tertiary

In this section we describe a model which is obtained by taking a Hill’s approximation of the spatial, circular, restricted four-body problem with oblate tertiary. The procedure goes as follows: we shift the origin of the coordinate system to $m_{3}$, perform a rescaling of the coordinates depending on $m_{3}^{1/3}$, write the associated Hamiltonian in the rescaled coordinates as a power series in $m_{3}^{1/3}$, and neglect all the terms of order $O(m_{3}^{1/3})$ in the expansion, since such terms are small when $m_{3}$ is small. Through this procedure the masses $m_{1}$ and $m_{2}$ are ‘sent to infinite distance’. This model is an extension of the classical Hill’s approximation of the restricted three-body problem, with the major differences that our model is a four-body problem, and takes into account the effect of the oblateness parameter $C_{20}=-J_{2}$ of the tertiary; compare with [Hil78, MS82, BGG15].

3.1. Hill’s approximation of the restricted four-body problem with oblate tertiary in shifted coordinates

The main result is the following:

Theorem 3.1**.**

Let us consider the Hamiltonian (2.38), let us shift the origin of the reference frame so that it coincides with $m_{3}$ and let us perform the conformal symplectic scaling given by

[TABLE]

Accordingly, we rescale the average radius of the tertiary as $R_{3}=m_{3}^{1/3}\rho_{3}$. Expanding the resulting Hamiltonian as a power series in $m_{3}^{1/3}$ and neglecting all the terms of order $O(m_{3}^{1/3})$ in the expansion, we obtain the following Hamiltonian describing the Hill’s four-body problem with oblate tertiary:

[TABLE]

where $v=\left(\frac{1}{1-\frac{3}{2}R_{3}^{2}C_{20}}\right)^{1/3}$, $\mu=\frac{m_{2}}{m_{1}+m_{2}}$, and $c:=\rho_{3}^{2}C_{20}/2=m_{3}^{-\frac{2}{3}}R^{2}_{3}C_{20}/2$.

Proof.

We start by shifting the origin of the coordinate system $(x,y,z)$ to the location of the mass $m_{3}$ (representing Hektor), via the change of coordinates

[TABLE]

The Hamiltonian corresponding to (2.38) becomes

[TABLE]

where $\bar{r}_{i}^{2}=(\xi-\bar{x}_{i})^{2}+(\eta-\bar{y}_{i})^{2}+\zeta^{2}=(\xi+x_{3}-x_{i})^{2}+(\eta+y_{3}-y_{i})^{2}+\zeta^{2}$, with $\bar{x}_{i}=x_{i}-x_{3}$, $\bar{y}_{i}=y_{i}-y_{3}$. Note that $\bar{r}_{3}=r_{3}$. Since $-\frac{1}{2}(x_{3}^{2}+y_{3}^{2})$ is a constant term, it plays no role in the Hamiltonian equations and it will be dropped in the following calculation.

Since $\sin{\phi}=\frac{z}{\bar{r}_{3}}$, we have

[TABLE]

We now perform the following conformal symplectic scaling with multiplier $m_{3}^{-2/3}$, given by (with a little abuse of notation, we call again the new variables $x$, $y$, $z$, $p_{x}$, $p_{y}$, $p_{z}$):

[TABLE]

Consistently with this scale change, we also introduce the scaling transformation of the average radius of the smallest body

[TABLE]

The choice of the power of $m_{3}$ is motivated by the fact that in this way the gravitational force becomes of the same order of the centrifugal and Coriolis forces (see, e.g., [MS82]).

The purpose of the subsequent calculation is that, after the above substitutions, we expand the resulting Hamiltonian as a power series in $m_{3}^{1/3}$ and neglect all the terms of order $O(m_{3}^{1/3})$ in the expansion.

We expand the terms $\displaystyle\frac{1}{\bar{r}_{1}}$ and $\displaystyle\frac{1}{\bar{r}_{2}}$ in Taylor series around the new origin of coordinates, obtaining

[TABLE]

where $P_{k}^{j}(x,y,z)$ is a homogeneous polynomial of degree $k$, for $j=1,2$.

The resulting Hamiltonian takes the form:

[TABLE]

In the following, we will disregard in (3.6) all terms which are of order of $m_{3}^{1/3}$, as in classical Hill’s theory of lunar motion ([MS82]).

We compute the first-degree polynomials $P^{i}_{1}$, for $i=1,2$,

[TABLE]

where $r_{i3}=\sqrt{(x_{i}-x_{3})^{2}+(y_{i}-y_{3})^{2}}=u=1$ represents the distance between the mass $m_{i}$ and $m_{3}$.

We now compute the contribution of the different terms in (3.6). Using that $m_{1}+m_{2}+m_{3}=1$, and that $m_{1}x_{1}+m_{2}x_{2}+m_{3}x_{3}=m_{1}y_{1}+m_{2}y_{2}+m_{3}y_{3}=0$, from (2.20), (2.21), (2.22), we obtain that

[TABLE]

where we used (2.36) and (3.5). The expression in (3.8) is $O(m_{3}^{1/3})$ so it will be omitted in the Hill approximation.

We compute the second degree polynomials $P^{i}_{2}$, for $i=1,2$,

[TABLE]

since $r_{13}=r_{23}=u=1$.

The corresponding terms in (3.6) yield

[TABLE]

Using (2.32) and that

[TABLE]

omitting the terms of order $O(m_{3}^{1/3})$, the expression (3.10) becomes

[TABLE]

The expression in (3.6) contains terms of the form

[TABLE]

which can be written in terms of positive powers of $m_{3}^{1/3}$, which are omitted in the Hill approximation.

The remaining terms in (3.6) are

[TABLE]

and they do not depend on $m_{3}$.

Therefore, when we omit all terms of order $O(m_{3}^{1/3})$ in (3.6), and taking into account that $\frac{1}{\omega^{2}}=1+O(m_{3}^{1/3})$, we obtain the following Hamiltonian

[TABLE]

where we used $\mu=m_{2}/(m_{1}+m_{2})$ and $c:=\rho_{3}^{2}C_{20}/2=m_{3}^{-\frac{2}{3}}R^{2}_{3}C_{20}/2$. ∎

We remark that a similar strategy was adopted in [MRPD01], where a Hill’s three body problem with oblate primaries has been considered.

We refer to the Hamiltonian in (3.11) as the Hill’s approximation. It can be thought of as the limiting Hamiltonian, when the primary and the secondary are sent at an infinite distance, and their total mass becomes infinite. It provides an approximation of the motion of the infinitesimal particle in an $O(m_{3}^{1/3})$ neighborhood of $m_{3}$. Remarkably, the angular velocity $\omega$ does not appear in the limiting Hamiltonian.

We introduce the gravitational potential as

[TABLE]

and the effective potential as

[TABLE]

The equations of motion associated to (3.11) can thus be written as:

[TABLE]

Remark 3.2*.*

In the case when $C_{20}=0$, we have that $v=1$ and the Hamiltonian in (3.11) is the same as the one obtained in [BGG15]. Also, its quadratic part coincides with the quadratic part of the expansion of the Hamiltonian of the restricted three-body problem centered at the Lagrange libration point $L_{4}$. Moreover, in the special case $\mu=0$ we obtain the classical lunar Hill’s problem after some coordinate transformation (see Section 3.2).

3.2. Hill’s four-body model applied to the Sun-Jupiter-Hektor system

In the case of the Sun-Jupiter-Hektor system, using the data from Section 2.1 in the above equations we obtain $\mu=m_{2}/(m_{1}+m_{2})=1.898\times 10^{27}/(1.989\times 10^{30}+1.898\times 10^{27})=0.0009533386$. Also, for Hektor we have $C_{20}=-0.476775$, the average radius of Hektor is $R_{3}=92$ km, and the mass of Hektor is $m_{H}=7.91\times 10^{18}$ kg. In the normalized units, where we use the average distance Sun-Jupiter $778.5\times 10^{6}$ km as the unit of distance, and the mass of Sun-Jupiter-Hektor $7.91\times 10^{18}+1.989\times 10^{30}+1.898\times 10^{27}=1.990898\times 10^{30}$ kg as the unit of mass, we have that the normalized average radius of Hektor is $R_{3}=92/(778.5\times{10}^{6})=1.18176\times 10^{-7}$ and the normalized mass of Hektor is $m_{3}=7.91\times 10^{18}/1.990898\times 10^{30}=3.97308\times 10^{-12}$. Hence, we obtain

[TABLE]

Also, $\rho_{3}=m_{3}^{-\frac{1}{3}}R_{3}=0.000746$.

We note that if we consider the restricted four-body problem (without the Hill’s approximation) described by the Hamiltonian (2.38), the oblateness effect is given by the coefficient $C^{\prime}=R_{3}^{2}C_{20}/2=-3.32921544\times 10^{-15}$, which is much smaller then $c$ in (3.14). As expected, the Hill’s approximation is acting like a ‘magnifying glass’ of the dynamics in a neighborhood of Hektor.

3.3. Hill’s approximation of the restricted four-body problem with oblate tertiary in rotated coordinates

In this section we write the Hamiltonian of Hill’s approximation of the restricted four-body problem with oblate tertiary in a rotating reference frame in which the primary and the secondary will be located on the horizontal axis.

Corollary 3.3**.**

The Hamiltonian (3.1) is equivalent, via a rotation of the coordinate axes that places the primary and the secondary on the $x$-axes, to the Hamiltonian

[TABLE]

where $\lambda_{2}$ and $\lambda_{1}$ are the eigenvalues corresponding to the rotation transformation in the $xy$-plane, and $\rho_{3}=m_{3}^{-1/3}R_{3}$.

Proof.

We perform a rotation on the $xy-$plane and re-write the Hamiltonian in (3.11) in the framework of the rotating coordinates, which are more suitable for the subsequent analysis. Since the rotation will be performed on the plane, we restrict the computations to the planar case. The planar effective potential restricted to the $xy-$plane is given by

[TABLE]

which can be written in matrix rotation as

[TABLE]

where $w=(x,y)^{T}$ and

[TABLE]

Notice that the matrix $M$ is symmetric, so its eigenvalues are real, the eigenvectors $v_{1}$ and $v_{2}$ are orthogonal, and the corresponding orthogonal matrix $C=\textrm{col}(v_{2},v_{1})$ defines a rotation in the $xy$-plane. We find the eigenvalues of $M$ by solving the characteristic equation:

[TABLE]

Since $0<\mu\leq\frac{1}{2}$, $\mu-\mu^{2}\leq 1/4$, and we have $1>1-v^{2}(4-v^{2})(\mu-\mu^{2})\geq 1-\frac{1}{4}v^{2}(4-v^{2})=\left(1-\frac{v^{2}}{2}\right)^{2}>0$. Thus $\lambda_{1},\lambda_{2}>0$ and $\lambda_{1}\neq\lambda_{2}$.

We notice that when $\mu=\frac{1}{2}$, the matrix $M$ is already a diagonal matrix, and the corresponding eigenvalues are $\lambda_{1}=\frac{3v^{2}}{4}$ and $\lambda_{2}=\frac{3(4-v^{2})}{4}$. Therefore, below we consider the case $\mu\neq\frac{1}{2}$ for which we proceed to compute the eigenvectors associated to $\lambda_{1}$ and $\lambda_{2}$.

The eigenvector $v_{1}$ such that $Mv_{1}=\lambda_{1}v_{1}$ and $\lVert v_{1}\rVert=1$ is given by

[TABLE]

where

[TABLE]

Similarly, the eigenvector $v_{2}$ such that $Mv_{2}=\lambda_{2}v_{2}$ and $\lVert v_{2}\rVert=1$ is given by

[TABLE]

where

[TABLE]

The equations of motion for the planar case can be written as

[TABLE]

where

[TABLE]

Consider the linear change of variable $w=C\bar{w}$ with $\bar{w}=(\bar{x},\bar{y})^{T}$. By substituting the new variable and multiplying $C^{-1}$ from the left, we obtain

[TABLE]

Notice that $D=C^{-1}MC$ is the diagonal matrix $D=\textrm{diag}(\lambda_{2},\lambda_{1})$, that is $\lVert C\bar{w}\rVert^{3}=\lVert\bar{w}\rVert^{3}$. Therefore the equation becomes

[TABLE]

Recall that $v_{1}=(v_{11},v_{12})^{T}$, $v_{2}=(v_{21},v_{22})^{T}$ and $C=\textrm{col}(v_{2},v_{1})$. Since $C$ is unitary, we have $C^{-1}=C^{T}$ and moreover

[TABLE]

A direct computation shows that $v_{12}v_{21}-v_{11}v_{22}=1$, which implies $C^{-1}\mathcal{J}C=\mathcal{J}$. Since $C^{-1}\mathcal{J}C=C^{T}\mathcal{J}C=\mathcal{J}$, the matrix $C$ is symplectic by definition. Therefore, the change of coordinates is symplectic. Thus, the equations of motion can be written as

[TABLE]

For $\mu\in[0,\frac{1}{2})$, we obtain the equations

[TABLE]

with

[TABLE]

From the expressions for $\bar{\Omega}_{\bar{x}}$ and $\bar{\Omega}_{\bar{y}}$, we notice the symmetry properties:

[TABLE]

Using these properties, we see that the equations (3.17) are invariant under the transformations $\bar{x}\rightarrow\bar{x}$, $\bar{y}\rightarrow-\bar{y}$, $\dot{\bar{x}}\rightarrow-\dot{\bar{x}}$, $\dot{\bar{y}}\rightarrow\dot{\bar{y}}$, $\ddot{\bar{x}}\rightarrow\ddot{\bar{x}}$ and $\ddot{\bar{y}}\rightarrow-\ddot{\bar{y}}$.

If we now go back to the spatial problem, we need to replace $\bar{\Omega}$ by

[TABLE]

writing $\bar{\Omega}(\bar{x},\bar{y},\bar{z})=\frac{1}{2}\bar{x}^{2}+\frac{1}{2}\bar{y}^{2}+\bar{U}(\bar{x},\bar{y},\bar{z})$, we can define $\bar{U}$ as

[TABLE]

In conclusion, the Hamiltonian in these new coordinates is given by the following expression (we omit the bars for $x$ and $y$ to simplify the notation):

[TABLE]

which coincides with (3.15). ∎

Remark 3.4*.*

If we let $C_{20}=0$ and $\mu=0$ in (3.15), we obtain the Hamiltonian for the classical lunar Hill problem, see, e.g., [MS82].

4. Linear stability analysis of the Hill four-body problem with oblate

tertiary

In this section we determine the equilibrium points associated to the potential in (3.18) and we analyze their linear stability.

4.1. The equilibrium points of the system

Considering the potential (3.18) (again we omit the bars for a simplified notation),

[TABLE]

we compute its derivatives as

[TABLE]

where $r=(x^{2}+y^{2}+z^{2})^{\frac{1}{2}}$.

To find the equilibrium points we have to solve the system

[TABLE]

In the above expressions $A$ and $B$ cannot simultaneously equal to [math] since $\lambda_{1}\neq\lambda_{2}$. Also, $A$ and $C$, or $B$ and $C$, cannot simultaneously equal to [math], since $A-C=\lambda_{2}+1-\frac{6c}{r^{5}}>0$ because $c<0$; a similar argument holds for $B$ and $C$. This implies that, for example, if $A=0$, then $B\neq 0$ and $C\neq 0$, so $y=z=0$ and $x$ is given by the equation $A=0$; the same reasoning applies for the other combinations of variables. Thus, all equilibrium points must lie on the $x$-, $y$-, $z$-coordinate axes. Precisely, we have the following results.

$i)$** Equilibrium points on the $x$-axis: **

In the case $A=0$, $B\neq 0$, $C\neq 0$, we must have $y=z=0$. From $A=0$ and $z=0$ we infer $\displaystyle h_{A}(r):=\lambda_{2}-\frac{1}{r^{3}}+\frac{3c}{r^{5}}=0$. We have $\displaystyle h^{\prime}_{A}(r)=\frac{3}{r^{4}}-\frac{15c}{r^{6}}>0$, since $c<0$; also, $\lim_{r\to 0}h_{A}(r)=-\infty$ and $\lim_{r\to\infty}h_{A}(r)=\lambda_{2}>0$. Hence, the equation $h_{A}(r)=0$ has a unique solution $r^{*}_{x}>0$, yielding the equilibrium points $(\pm r^{*}_{x},0,0)$.

$ii)$** Equilibrium points on the $y$-axis: **

In the case $B=0$, $A\neq 0$, $C\neq 0$, we must have $x=z=0$. From $B=0$ and $z=0$ we infer $\displaystyle h_{B}(r):=\lambda_{1}-\frac{1}{r^{3}}+\frac{3c}{r^{5}}=0$. We have $\displaystyle h^{\prime}_{B}(r)=\frac{3}{r^{4}}-\frac{15c}{r^{6}}>0$, since $c<0$; also, $\lim_{r\to 0}h_{B}(r)=-\infty$ and $\lim_{r\to\infty}h_{B}(r)=\lambda_{1}>0$. Hence, the equation $h_{B}(r)=0$ has a unique solution $r^{*}_{y}>0$, yielding the equilibrium points $(0,\pm r^{*}_{y},0,0)$.

$iii)$** Equilibrium points on the $z$-axis: **

In the case $C=0$, $A\neq 0$, $B\neq 0$, we must have $x=y=0$, so $z=\pm r$. Hence $C=0$ implies $-1-\frac{1}{r^{3}}+\frac{9c}{r^{5}}-\frac{15cr^{2}}{r^{7}}=-1-\frac{1}{r^{3}}-\frac{6c}{r^{5}}=(-r^{5}-r^{2}-6c)/r^{5}=0$. Let $h_{C}(r)=-r^{5}-r^{2}-6c$. We have $h^{\prime}_{C}(r)=-5r^{4}-2r<0$; also, $\lim_{r\to 0}h_{C}(r)=-6c>0$ and $\lim_{r\to+\infty}h_{C}(r)=-\infty$. Hence, the equation $h_{C}(r)=0$ has a unique solution $r^{*}_{z}>0$, yielding the equilibrium points $(0,0,\pm r^{*}_{z})$.

In the case of the Sun-Jupiter-Hektor system, in normalized units, we obtain $\lambda_{1}=0.0021444999866622183$, $\lambda_{2}=2.997855500013338$, and the equilibrium points location are given by the following figures.

[TABLE]

We remark that in the case of the Hill’s four body problem without a non-oblate tertiary, the $x$-axis equilibria and the $y$-axis equilibria also exist, see [BGG15]; their locations, in the case of Hektor, are very close to the ones in the case of an oblate tertiary. Precisely, we have the following results.

[TABLE]

This result leads us to conclude that the $x$-axis equilibria and the $y$-axis equilibria for the Hill’s problem with oblate tertiary are continuations of the ones for the Hill’s problem with non-oblate tertiary. On the other hand, the $z$-axis equilibria do not exist for the Hill’s problem with non-oblate tertiary, so they are new features of the Hill’s problem with oblate tertiary.

To summarize, the Hill’s three-body problem has 2 equilibrium points, Hill’s four-body problem has 4 equilibrium points, and the Hill’s four-body problem with oblate tertiary has 6 equilibrium points.

In Fig. 3 we plot the dependence on the $c$ of the distance from $z$-axis equilibrium point to the origin (in km), when we let the parameter $C_{20}$ range between -0.001 and -0.95. The estimates on the axes of Hektor are $a=208$ km, $b=65.5$ km, $c=60$ km ([Des15]). Using the ellipsoid model, we find $C_{20}=-0.476775$ and hence $z\simeq 100$ km, which lies outside the asteroid. Using $C_{20}=-0.15$ as provided by [MDCR*+*14], we obtain $x\simeq 62$ km, which basically coincides with the surface of the asteroid. We notice that the vertical equilibrium is not outside the Brillouin sphere, which corresponds to the region where the spherical harmonic series expansion is convergent. Inside the Brillouin sphere the series is divergent, if the shape is an ellipsoid. Given that the true shape of the asteroid is unknown, we cannot determine precisely the region where the spherical harmonic series is convergent or divergent and, hence, we cannot decide whether the $z$-equilibrium is indeed real, but just postulate its existence.

To convert to real units, the distances from the equilibrium points to the center need to be multiplied by $m_{3}^{1/3}$ – due to the rescaling involved in the Hill procedure –, and by the unit of distance which in this case is the distance Sun-Jupiter. It follows that the $x$-axis equilibrium points are at a distance of $85,512.774$ km from Hektor, the $x$-axis equilibrium points are at a distance of $956,149.406$ km, and the $z$-axis equilibrium points are at a distance of $110.028$ km. As the smallest semi-minor axis of Hektor is $60$ km, the $z$-axis equilibrium points are outside the body of the asteroid. It seems though that for many other asteroids the $z$-axis equilibrium points are located inside their body.

Remark 4.1*.*

These $z$-axis equilibria also appear in the case of the motion of a particle in a geopotential field (see, e.g., [CG14]). From that model it can be derived that the distance from the $z$-axis equilibrium points to the center is given by $\hat{r}_{z}=R_{3}(-3C_{20})^{1/2}$. When we apply this formula in the case of Hektor, the numerical result is very close to the one found above. This formula can be derived from the equation $h_{C}(r)=0$ if we drop the term $r^{5}$.

4.2. Linear stability of the equilibrium points

We study the linear stability of the equilibrium points in the case of Hektor.

The Hamiltonian (3.15) yields the following system of equations

[TABLE]

where $\Omega$ is the effective potential given by (3.18) (again, we omit the overline bar on the variables).

The second order derivatives of $\Omega$ are given by

[TABLE]

The Jacobian matrix describing the linearized system is

[TABLE]

Since the equilibria are of the form $(\pm r^{*}_{x},0,0)$, $(0,\pm r^{*}_{y},0)$, $(0,0,\pm r^{*}_{z})$, the mixed second order partial derivatives $\Omega_{xy}$, $\Omega_{xz}$, $\Omega_{yz}$ vanish at each of the equilibrium points.

Hence Jacobian matrix (4.3) evaluated at the equilibria is of the form:

[TABLE]

The characteristic equation of (4.4) is

[TABLE]

The stability of the equilibria depends on the signs of $\Omega_{zz}$ and of $A$, $B$, and $D$.

We find the following stability character of the equilibrium positions in the case of the Sun-Jupiter-Hektor system:

$i)$** Eigenvalues of $x$-axis equilibria at $(\pm 0.6935267570,0,0)$: **

[TABLE]

Stability type: center $\times$ center $\times$ saddle.

$ii)$** Eigenvalues of $y$-axis equilibria at $(0,7.7545747196,0)$: **

[TABLE]

Stability type: center $\times$ center $\times$ center.

$iii)$** Eigenvalues of $z$-axis equilibria at $(0,0,\pm 0.0008923544)$: **

[TABLE]

Stability type: center $\times$ complex saddle.

Note that for the $z$-axis equilibria the imaginary part of the ‘Krein quartet’ of eigenvalues of $z$-axis equilibria is approximately $\pm 1$. This means that the infinitesimal motion around the equilibrium point is close to the $1:1$ resonance with the rotation of the primary and the secondary. In Fig. 4 we show that for a range of $r^{*}_{z}$ values between $z=0.000892354498497342$ (corresponding to the $c$ value for Hektor) and $z=0.01$ (corresponding to $c=-1.666668333\times 10^{-5}$), the real part and the imaginary part of the ‘Krein quartet’ of eigenvalues; the imaginary part stays close to $\pm 1$.

In Section 4.2.1 we will provide an analytic argument that the real part of the ‘Krein quartet’ of eigenvalues is always non-zero, and the imaginary part is close to $\pm 1$ for $r^{*}_{z}$ sufficiently small; this result will help us to explain the behavior observed in Fig. 4.

In the sequel we give a more detailed analysis of the linear stability of all equilibria, for a wide range of parameters $\mu$ and $c$.

4.2.1. Linear stability of the equilibria on the $z$-axis

The $z$-axis equilibrium points are of the form $(0,0,\pm r^{*}_{z})$, with

[TABLE]

which yields

[TABLE]

Evaluating $\Omega_{xx}$, $\Omega_{yy}$, $\Omega_{zz}$ at the equilibrium point yields:

[TABLE]

Substituting in (4.7) we have

[TABLE]

Using (3.16) and denoting $d:=\sqrt{1-v^{2}(4-v^{2})(\mu-\mu^{2})}$ we can write

[TABLE]

Also for $c=0$ we have $d_{0}=\sqrt{1-3(\mu-\mu^{2})}$ and

[TABLE]

This is in agreement with the results in [BGG15].

For future reference, we expand $d$ as a power series in the parameter $c$ as

[TABLE]

where the coefficient $d_{1}$ can be obtained can be obtained from the Taylor’s theorem around $c=0$ as

[TABLE]

From the characteristic equation (4.5), we obtain that the pair of eigenvalues $\rho_{1,2}=\pm(\Omega_{zz})^{1/2}$ is purely imaginary, since by (4.8), $\Omega_{zz}<0$.

The ‘Krein quartet’ eigenvalues are given by

[TABLE]

where

[TABLE]

Then we have

[TABLE]

Since $-A>0$ and $D<0$, we obtain that the eigenvalues $\rho_{3,4,5,6}$ are complex numbers, non-real, non-purely-imaginary, for all parameter values. Let $\rho=a+ib$ be such that $\rho^{2}=-\frac{A}{2}\pm\frac{\sqrt{4B-A^{2}}}{2}i:=\alpha+i\beta$. We have

[TABLE]

To show that $b$ is approximately $\pm 1$, or $b^{2}\approx 1$, for $r^{*}_{z}\approx 0$, note that

[TABLE]

where $\Upsilon:=v^{2}(4-v^{2})(\mu-\mu^{2})$. Since

[TABLE]

we have that $\lim_{r^{*}_{z}\to 0}b^{2}=-\frac{3}{4}+\frac{7}{4}=1$, so $b^{2}\approx 1$ for $r^{*}_{z}\approx 0$, as in the case of Hektor.

We obtained the following result:

Proposition 4.2**.**

Consider the equilibria on the $z$-axis. For $\mu\in(0,1/2]$, $\Omega_{zz}$, $A$ and $D$ are negative. Consequently, one pair of eigenvalues is purely imaginary, and the two other pairs of eigenvalues are complex conjugate, with the imaginary part close to $\pm i$ for $c$ negative and sufficiently small. The linear stability is of center $\times$ complex-saddle type.

4.2.2. Linear stability of the equilibria on the $y$-axis

The $y$-axis equilibrium points are of the form $(0,\pm r^{*}_{y},0)$, with

[TABLE]

which yields

[TABLE]

Evaluating $\Omega_{xx}$, $\Omega_{yy}$, $\Omega_{zz}$ at the equilibrium point yields:

[TABLE]

Substituting $c$ from (4.15) we obtain

[TABLE]

We also expand $r^{*}_{y}$ as a power series in the parameter $c$ as

[TABLE]

where $\pm r_{y0}$ is the position of the $y$-equilibrium in the case when $c=0$, which is given by $r_{y0}^{3}=1/\lambda_{10}$; this is in agreement with [BGG15]. The computation of $r_{y1}$ yields

[TABLE]

with $d_{1}$ as in (4.12).

We will also need $\frac{1}{(r^{*}_{y})^{3}}$ as a power series in the parameter $c$

[TABLE]

and a simple calculation yields

[TABLE]

For $\mu=1/2$, we have $d_{0}=\frac{1}{2}$, $\lambda_{10}=\frac{3}{4}$, $d_{1}=-m_{3}^{2/3}$, $r_{y0}=\left(\frac{4}{3}\right)^{1/3}$. It is easy to see that dominant part $d_{0}$ of $d$ is a strictly decreasing function with respect to $\mu\in(0,1/2]$ and takes values in $[1/2,1)$. The dominant part $\lambda_{10}$ of $\lambda_{1}$ is increasing with respect to $\mu\in(0,1/2]$ and takes values in $(0,3/4]$. Also, the dominant part $r_{y0}$ of $r^{*}_{y}$ is a strictly decreasing function in $\mu\in(0,1/2]$, where $r_{y0}(1/2)=\sqrt[3]{4/3}$ and $r_{y0}\rightarrow\infty$ when $\mu\rightarrow 0$; as a consequence the values of $r_{y0}$ are in the interval $[\sqrt[3]{4/3},\infty)$.

From (4.16) we have

[TABLE]

since $r^{*}_{y}>0$ and $c$ is negative. Therefore, $\Omega_{zz}<0$ for all admissible values of $\mu$.

For $A=4-\Omega_{xx}-\Omega_{yy}$, using (4.10) and the expansions (4.11) and (4.20) we obtain

[TABLE]

for $c$ small.

For $B=\Omega_{xx}\Omega_{yy}$ using (4.10) and the expansions (4.11) and (4.20) we obtain

[TABLE]

for $c$ small.

For $D=A^{2}-4B$, using (4.10) and the expansions (4.11) and (4.20) we have

[TABLE]

For $\mu\approx 0$ we have $D\approx 1+O(c)$ and for $\mu=1/2$ we have $D=-\frac{215}{16}+O(c)$. The intermediate value theorem implies that $D$ changes its sign from positive to negative for $\mu\in(0,1/2]$, provided $c$ is small. We have thus proved the following result:

Proposition 4.3**.**

Consider the equilibria on the $y$-axis. For $\mu\in(0,1/2]$ and for the parameter $c$ negative and small enough, $\Omega_{zz}$ is always negative, the coefficients $A$ and $B$ are always positive, and the value of the discriminant $D$ changes from positive to negative values. Consequently, one pair of eigenvalues is always purely imaginary, and there exists $\mu_{*}$, depending on $c_{20}$, where the other two pairs of eigenvalues change from being purely imaginary to being complex conjugate. The linear stability changes from center $\times$ center $\times$ center type to center $\times$ complex-saddle type.

4.2.3. Linear stability of the equilibria on the $x$-axis

The $x$-axis equilibrium points are of the form $(\pm r^{*}_{x},0,0)$, with

[TABLE]

which yields

[TABLE]

Evaluating $\Omega_{xx}$, $\Omega_{yy}$, $\Omega_{zz}$ at the equilibrium point yields:

[TABLE]

Substituting $c$ from (4.23) we obtain

[TABLE]

We expand $r^{*}_{x}$ as a power series in the parameter $c$ as

[TABLE]

where $\pm r_{x0}$ is the position of the $x$-equilibrium in the case when $c=0$, which is given by $r_{x0}^{3}=1/\lambda_{20}$; see [BGG15]. The computation of $r_{x1}$ yields

[TABLE]

We will also need $\frac{1}{(r^{*}_{x})^{3}}$ as a power series in the parameter $c$

[TABLE]

and a simple calculation yields

[TABLE]

From (4.24) we have

[TABLE]

since $r^{*}_{x}>0$ and $c<0$. Therefore, $\Omega_{zz}<0$ for all admissible values of $\mu$.

For $A=4-\Omega_{xx}-\Omega_{yy}$, using (4.10) and the expansions (4.11) and (4.28), we obtain

[TABLE]

for $c$ small.

For $B=\Omega_{xx}\Omega_{yy}$, using (4.10) and the expansions (4.11) and (4.28). we obtain

[TABLE]

for $c$ small.

For $D=A^{2}-4B$, using (4.10) and the expansions (4.11) and (4.28) we have

[TABLE]

for $c$ small.

We have proved the following result:

Proposition 4.4**.**

Consider the equilibria on the $x$-axis. For $\mu\in(0,1/2]$ and for parameter $c$ negative and small enough, $\Omega_{zz}$ is negative, $A$ and $B$ are negative, and the value of the discriminant $D$ is always positive. Consequently, two pairs of eigenvalues are purely imaginary, and one pair of eigenvalues are real (one positive and one negative). The linear stability is of center $\times$ center $\times$ saddle type.

Acknowledgements

This material is based upon work supported by the National Science Foundation under Grant No. DMS-1440140 while A.C. and M.G. were in residence at the Mathematical Sciences Research Institute in Berkeley, California, during the Fall 2018 semester.

This research was carried out (in part) at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration and funded through the Internal Strategic University Research Partnerships (SURP) program.

A.C. was partially supported by GNFM-INdAM and acknowledges the MIUR Excellence Department Project awarded to the Department of Mathematics, University of Rome Tor Vergata, CUP E83C18000100006. M.G. and W-T.L. were partially supported by NSF grant DMS-0635607 and DMS-1814543.

We are grateful to Rodney Anderson, Edward Belbruno, Ernesto Perez-Chavela, and Pablo Roldán for discussions and comments.