1 Introduction

In recent decades, space exploration has shown special interest in asteroid missions, such as NASA’s OSIRIS-REx, Japan’s Hayabusa2 and ESA’s Hera. To carry out the space missions, it is fundamental to understand the dynamics around these bodies, which present a wide variety of shapes and compositions. The asteroids characteristics depend mainly on their origin, although a common feature found in them is their elongated shape and asymmetric mass distribution.

There have been several approaches in the literature to model the potentials of irregular asteroids. Complex models, such as polyhedra (Werner and Scheeres 1997) and mascons (Geissler et al. 1996), provide high precision, but at the cost of long computing times and numerous parameters to handle. Other simplified models, such as dipole (Bartczak et al. 2006), triaxial homogeneous ellipsoids (Bartczak and Breiter 2003), and segment (Duboshin 1959; Breiter et al. 2005; Riaguas et al. 2001), allow obtaining analytical results on the dynamics. Although these simplified models fail to describe the dynamics close to the asteroid.

Zeng et al. (2018) proposed the dipole-segment model, consisting of a massive segment with two point masses at the ends. In that work, they proved that the dipole-segment model is more realistic when modeling the gravitational fields of asteroids. Later, Elipe et al. (2021) computed families of periodic orbits for this model, but only in the case of symmetric mass distributions.

Scheeres et al. (2000b) developed a method for adjusting the density of an asteroid using tetrahedra. Since then, several studies have estimated the density of celestial bodies using this method from space missions data Ledbetter et al. (2021), Tricarico (2013), Yin et al. (2021), Kanamaru et al. (2019), McMahon et al. (2018), Takahashi and Scheeres (2014). The results of these studies show that asteroids can have a heterogeneous and asymmetric density profile. However, most of the proposed models for analyzing the dynamics focus on the geometry of the asteroid, rather than its internal structure. Therefore, there is a gap in the literature on simplified models of asteroids with variable and continuous density.

In this article we study families of periodic orbits around elongated asteroids using two different models. We first analyze the dipole-segment model (ADS), in which the irregularities of the asteroids are captured by the shape of the asteroid by adding point masses at the ends of the segment. Then, we propose the variable density segment model (VDS), in which the irregularity comes from continuous changes in the internal density of the asteroid. Both models are compared qualitatively and analyzed through the families of periodic orbits.

2 Dynamical model

In this section we explore the dipole-segment and variable density models with bi-parametric mass distributions Throughout the section we will see that the variable segment model can effectively mimic the distribution of a dipole-segment. To achieve this similarity, it is necessary for the density distributions to meet certain properties. In particular, two configurations have been studied: those that have the same center of mass and variance, and those that have the same center of mass and the same amount of mass on each side of this point. These criteria ensure that the qualitative behavior of periodic orbits in both models is equivalent.

2.1 Asymmetric dipole model

Hence, we consider necessary to start with a knowledge of dynamical aspects in some simple cases. In Riaguas et al. (2001) it is analyzed the stability of equilibrium points of a fixed segment, as a first approach of an elongated asteroids. More recently, Zeng et al. (2018) proposed a dipole-segment approach. They analyzed the symmetric case computing the equilibrium points and stability. In this section, we propose to work with a non-symmetric dipole-segment model.

We will consider a model, formed by a massive dipole and a mass segment, denoted ADS, as may be found in Zeng et al. (2018) and Elipe et al. (2021). This model is constituted by a massive dipole formed by two masses \(m_{1}\) and \(m_{2}\) each one at a different end of a massive segment of mass \(m\) and total length \(\ell \). This body rotates with constant angular velocity \(\omega \) around its center of mass. We align the segment with the \(Ox\) axis, make the center of mass of the ADS body and the origin of coordinates coincide and the rotation axis is \(Oz\). Following closely notations in Elipe et al. (2021) and Riaguas et al. (2001) we write

$$\begin{aligned} &M=m+m_{1}+m_{2}, \quad \mu ={\displaystyle \frac{m_{2}}{m_{1}+m_{2}}}, \\ & 1-\mu ={\displaystyle \frac{m_{1}}{m_{1}+m_{2}}}, \quad \mu _{s}={ \displaystyle \frac{m}{M}}. \end{aligned}$$

Therefore, masses attached at the segment ends, \(m_{1}\) y \(m_{2}\) are in points of space \((-\ell _{1}, 0,0)\) and \((\ell _{2},0,0)\) respectively. Then \(\ell _{1}+\ell _{2}=\ell \), and it holds that

$$ {\displaystyle \frac{m_{1}}{M}}=(1-\mu )(1-\mu _{s}),\quad { \displaystyle \frac{m_{2}}{M}}=\mu (1-\mu _{s}). $$

leading to

$$ 0\le \mu _{s}\le 1,\quad 0\le \mu \le 1. $$

Choosing the ADS lying on the \(Ox\) axis, revolving around its center of mass and with rotation axis the \(Oz\) allows us to express the coordinates of the segment ends as it holds that

$$ 0= {\displaystyle \frac{1}{M}}\left ( -m_{1}\ell _{1}+m_{2}\ell _{2}+ \int _{-\ell _{1}}^{\ell _{2}}\lambda x\,\mathrm{d}{x} \right ) $$

with \(\lambda \) the constant linear mass density of the segment \(\lambda \ell =\lambda (\ell _{1}+\ell _{2})=m \). Therefore

$$\begin{aligned} 0&= {\displaystyle \frac{1}{M}}\left ( -m_{1}\ell _{1}+m_{2}\ell _{2}+{ \displaystyle \frac{m}{2(\ell _{1}+\ell _{2})}}(\ell _{2}^{2}-\ell _{1}^{2}) \right ) \\ &= {\displaystyle \frac{1}{M}}\left ( -m_{1}\ell _{1}+m_{2}\ell _{2}+{ \displaystyle \frac{m}{2}}(\ell _{2}-\ell _{1}) \right ) \end{aligned}$$

From here

$$ \textstyle\begin{array}{l} \ell _{2}\left (m_{2}+{\displaystyle \frac{m}{2}}\right )=\ell _{1} \left (m_{1}+{\displaystyle \frac{m}{2}}\right ), \\ {\displaystyle \frac{\ell _{2}}{\ell }}={\displaystyle \frac{1}{M}} \left (m_{1}+{\displaystyle \frac{m}{2}}\right )=(1-\mu )(1-\mu _{s})+{ \displaystyle \frac{1}{2}} \mu _{s}, \\ {\displaystyle \frac{\ell _{1}}{\ell }}={\displaystyle \frac{1}{M}} \left (m_{2}+{\displaystyle \frac{m}{2}}\right )=\mu (1-\mu _{s})+{ \displaystyle \frac{1}{2}} \mu _{s}. \end{array} $$

The gravitational potential due to the ADS can be written as

$$\begin{aligned} U&={-\mathcal{G}} M\bigg( {\displaystyle \frac{(1-\mu )(1-\mu _{s})}{r_{1}}}+{\displaystyle \frac{\mu (1-\mu _{s})}{r_{2}}}\\ &\quad +{\displaystyle \frac{\mu _{s}}{\ell }} \log \left [ {\displaystyle \frac{r_{1}+r_{2}+\ell }{r_{1}+r_{2}-\ell }}\right ]\bigg) \\ &={-\mathcal{G}} M\bigg( {\displaystyle \frac{2\ell _{2}-\ell \mu _{s}}{2 \ell r_{1}}}+{\displaystyle \frac{2\ell _{1}-\ell \mu _{s}}{2 \ell r_{2}}}\\ &\quad +{\displaystyle \frac{\mu _{s}}{\ell }}\log \left [ {\displaystyle \frac{r_{1}+r_{2}+\ell }{r_{1}+r_{2}-\ell }}\right ]\bigg) \\ &={\displaystyle \frac{{-\mathcal{G}} M}{\ell }}\bigg( {\displaystyle \frac{2\ell _{2}-\ell \mu _{s}}{2 r_{1}}}\\ &\quad +{\displaystyle \frac{2\ell _{1}-\ell \mu _{s}}{2 r_{2}}}+\mu _{s}\log \left [ { \displaystyle \frac{r_{1}+r_{2}+\ell }{r_{1}+r_{2}-\ell }}\right ] \bigg) \end{aligned}$$

where \({\mathcal{G}}\) is the gravitational constant and \(r_{1}\) y \(r_{2}\) are given by

$$ \textstyle\begin{array}{l} r_{1}=(x+\ell _{1})^{2}+y^{2}+z^{2}, \\ r_{2}=(x-\ell _{2})^{2}+y^{2}+z^{2}. \end{array} $$

Taking into account that the acceleration of a particle in a rotating frame under the gravitational field of the ADS is

$$ \ddot{\vec{x}}+2\vec{\omega} \times \dot{\vec{x}} +\vec{\omega}\times ( \vec{\omega} \times \vec{x}) +\dot{\vec{\omega} }\times \vec{x}=\nabla _{\vec{x}} U(\vec{x}) , $$

with \(U\) the gravitational potential of the ADS. The usual case for elongated asteroids is to be uniformly rotating about its maximum moment of inertia (Scheeres et al. 2000c). Therefore, if we assume that the reference frame rotates uniformly with angular velocity \(\omega \), the cartesian components of the previous equation are

$$ \textstyle\begin{array}{l} \ddot{x}-2\omega \dot{y} = \omega ^{2} x + U_{x}, \\ \ddot{y}+2\omega \dot{x} = \omega ^{2} y + U_{y}, \\ \ddot{z} = U_{z}. \end{array} $$
(1)

Then, the equations of motion in the rotating frame become time invariant, and in consequence, the system has a first integral, the Jacobi constant,

$$ C=-2U(x,y,z)+(x^{2}+y^{2})-({\dot{x}}^{2}+{\dot{y}}^{2}+{\dot{z}}^{2}) $$
(2)

As usual, we introduce the effective potential \(W\) as

$$ W(x,y,z)= \frac {\omega ^{2}}{2}(x^{2}+y^{2})-U(x,y,z). $$
(3)

Doing a change of variables, we can always take the length of the asteroid as the reference unit. Therefore, from now on we will consider that \(\ell = 1\).

The reader may notice the ADS is similar to the Restricted Three Body Problem, but with the difference that the space located between the primaries is now occupied by the asteroid. Equilibria are found by zeroing the right-hand side of Eq. (1). Four points are obtained. Two of them, named the collinear points, are located in the \(Ox\)-axis where the asteroid is lying, and are denoted by \(L_{1}\) and \(L_{2}\). The other two, named the triangular points, form an equilateral triangle with the extreme positions of the asteroids, they are denoted \(L_{3}\) and \(L_{4}\). The coordinates of collinear and triangular equilibrium points depend on the parameters of the system. In Fig. 1 it is depicted the location of the equilibria for the values of the parameters \(\mathcal{G M = 1}\), \(\omega = 1\), \(\mu =0.3\) and \(\mu _{S}=0.5\), which are the values that correspond to an asymmetric distribution that will be used in the numerical analysis presented in Sect. 3. Although the case \(\mu =0.3\) is used, the topology is the same for all values of the mass parameter, except for the degenerated case \(\mu =0\) (two-body problem).

Fig. 1
figure 1

Equilibria distribution for an asymmetric dipole-segment, \(\mu =0.3\), \(\mu _{s}=0.5\)

Full size image

If \(C_{L_{i}}\) is the value of Jacobi constant given by Eq. (2), at the equilibrium points we get that

$$ C_{L_{1}}>C_{L_{1}}>C_{L_{3}}=C_{L_{4}}. $$

Figure 2 represents the zero velocity curves in the equatorial plane decreasing from a high to a low Jacobi constant level. The motion is forbidden in the blue regions (\(C-W<0\)), the white area is where the asteroid is located, and the yellow area is the region where the motion is allowed (\(C-W>0\)). Figure 2(a) shows that a body with this energy is captured in the gravity of the asteroid mass, or if it is located outside the system, it is prohibited from entering by the centrifugal force of the rotating asteroid. The energy levels of the equilibrium points display different topology of the zero velocity curves, and a transition from one to two regions of possible motions occurs.

Fig. 2
figure 2

Zero velocity surfaces for different values of the Jacobi constant

Full size image

2.2 Segment model with variable density

The model under consideration in this section involves a mass segment with variable density (VDS), represented by a two degree polynomial linear density, rotating around its center of mass. This model is an alternative to the dipole segment to study the dynamics around irregular asteroids. Assuming the segment is of unit length, its linear density can be modeled as

$$ \sigma \left ( \nu \right ) = M\left ( a_{0} + a_{1} \nu + a_{2}\nu ^{2} \right ) \quad \nu \in \left [ 0,1 \right ], $$
(4)

where \(M>0\) is the total mass of the segment, \(a_{1},a_{2} \in \mathbb{R}\) and \(a_{0} = \frac{1}{6}(6 - 3 a_{1} - 2 a_{2})\) such that

$$ \int _{0}^{1} \sigma (\nu ) d \nu =M. $$
(5)

In addition, the density must be positive, which translates into the following conditions for \(\sigma \):

  • \(\sigma (0) > 0 \Rightarrow a_{0}>0\),

  • \(\sigma (1) > 0 \Rightarrow a_{0} + a_{1} + a_{2}>0\),

  • If \(a_{2} \neq 0\) and \(\dfrac{-a_{1}}{2a_{2}} \in \left ( 0,1 \right )\) then \(a_{0}-\dfrac{a_{1}^{2}}{4 a_{2}} > 0\).

These conditions outline a feasibility domain for the parameters \(a_{1}\) and \(a_{2}\), denoted as \(\Omega \), and illustrated in Fig. 3. Using a degree two polynomial provides a good balance between simplicity and adaptability in modeling the density of the rod, capturing symmetry and convexity properties. The density function is symmetric only when the coefficients of \(\sigma \) satisfy \(a_{1} + a_{2} = 0\). If \(a_{2} > 0\), the density function is convex, indicating the mass distribution of the rod accumulates at the ends. On the other hand, if \(a_{2} < 0\), the density function is concave, indicating mass accumulation in the middle of the segment.

Fig. 3
figure 3

Domain of parameters \((a_{1},a_{2})\) resulting in positive density

Full size image

Following the line of the ADS, we consider the rotational motion of the VDS with constant angular velocity \(\omega \) around its center of mass \(\ell _{1} = \frac{1}{12}\left ( 6+a_{1}+a_{2} \right )\). To establish the effective potential, the segment is placed at \(x \in \left [ -\ell _{1} ,1-\ell _{1} \right ]\), \(y=z=0\), ensuring that the center of mass lies at the origin of coordinates. Subsequently, the effective potential is derived as follows:

$$ W=\frac{\omega ^{2}}{2}\left (x^{2}+y^{2}\right )+\frac{\mathcal{G M}}{24} \left (f+g \cdot \log \left (\frac{s+1}{s-1}\right )\right ), $$
(6)

where

  • \(s=r_{1}+r_{2}\),

  • \(d=r_{1}-r_{2}\),

  • \(f=6 a_{2}\left (1-3 d^{2}\right ) s-24\left (a_{1}+a_{2}\right ) d\),

  • \(g=24+a_{2}-3 a_{2} d^{2}-\frac{1}{2} s f\).

Since the problem configuration is the same as in ADS, both problems share the same Lagrangian structure. Therefore, the equations of motion and Jacobi integral of VDS are formally the same. The VDS model, as seen in the uniform rod Riaguas et al. (2001) and dipole models Zeng et al. (2018); Elipe et al. (2021), has four equilibrium points. It is easy to show that these equilibria are located in the \(Oxy\) plane due to the symmetry of equations. However, given the complexity of the equations of motion, it is unfeasible to obtain the position of these equilibria analytically. A numerical approach involving an exhaustive scan of \(\Omega \) was employed, leading to the verification that, for any allowed values of \(a_{1}\) and \(a_{2}\), there exist two equilibrium points collinear with the rod, alongside two additional triangular equilibria. It should be noted that only in the case of symmetric density, the triangular equilibria are on the \(Oy\) axis.

2.3 Qualitative equivalence between models

This section aims to replicate the dynamic behavior around the ADS via the VDS. The applied approach involves using density parameters \(a_{1}\) and \(a_{2}\) to fit mass distribution properties of the ADS. Additionally, a case study will be examined to comprehensively explore the qualitative behavior of the VDS model related to the ADS.

Let’s assume both segments have the same length, total mass, and angular velocity. Thus, the problem of fitting the variable density can be interpreted as searching for values of \(a_{1}\) and \(a_{2}\) that replicate the properties of an asymmetric dipole mass distribution. We considered three properties to establish a resemblance between the VDS and ADS: mass on both sides of the center of mass, location of the center of mass, and variance. Because both models have the center of mass located at the origin, the first property refers to the mass portion between \(-\ell _{1}\) and 0 (equivalently, between 0 and \(1-\ell _{1}\)). This condition translates into the following equation:

$$\begin{aligned} &(1-\mu ) (1-\mu _{s})+\mu _{s} \left (\mu (1-\mu _{s})+ \frac{\mu _{s}}{2}\right ) \\ &\quad = \frac{a_{2} (a_{1}+a_{2}+6)^{3}}{5184}+ \frac{1}{288} a_{1} (a_{1}+a_{2}+6)^{2} \\ &\qquad +\frac{1}{72} (-3 a_{1}-2 a_{2}+6) (a_{1}+a_{2}+6). \end{aligned}$$
(C0)

The second property refers to the location of the center of mass in the rod, that is, the value of \(\ell _{1}\) in each model. The condition for these values to coincide is

$$ \mu (1-\mu _{s})+\frac{\mu _{s}}{2} = \frac{1}{12} (a_{1}+a_{2}+6) $$
(C1)

Finally, the variance is defined as the mean value of the square deviation from the center of mass. If we adjust the variance of both rods to match, we obtain the third condition:

$$\begin{aligned} & \frac{1}{720} \left (4 (a_{2}+15)-5 (a_{1}+a_{2})^{2}\right ) \\ &\quad = (1- \mu ) (1-\mu _{s}) \left (\mu (1-\mu _{s})+\frac{\mu _{s}}{2}\right )^{2} \\ &\qquad +\mu (1-\mu _{s}) \left (-\mu (1-\mu _{s})-\frac{\mu _{s}}{2}+1 \right )^{2} \\ &\qquad +\mu _{s} \left (\frac{1}{3} \left (-\mu (1-\mu _{s})- \frac{\mu _{s}}{2}+1\right )^{3}\right . \\ &\qquad \left .-\frac{1}{3} \left (-\mu (1-\mu _{s})-\frac{\mu _{s}}{2} \right )^{3}\right ) \end{aligned}$$
(C2)

Notice the properties under consideration involve distribution moments of orders 0, 1, and 2. The moment of order 0 represents the total mass, the moment of order 1 represents the center of mass, and the moment of order 2 represents the variance or spread. With only two free parameters, it is impossible to simultaneously satisfy all three conditions above. Therefore, a decision must be made regarding the pair of conditions that should be prioritized. Given the computational complexity involved, the focus of this study was limited to exploring the pairs C0-C1 and C1-C2.

Once the pair of conditions has been chosen, the next step involves solving a system of two equations - two unknowns for \(a_{1}\) and \(a_{2}\). Note that even after obtaining the values of \(a_{1}\) and \(a_{2}\), there is no guarantee that these values will fall within the feasible region. Lets call \(\Psi = \left [0,1\right ] \times \left [0,1\right ]\) the parametric domain of the ADS mass distribution, that is,

$$ (\mu , \mu _{s}) \in \Psi . $$

Depending on the chosen pair of conditions, each point in \(\Psi \) corresponds to a value of \((a_{1}, a_{2})\) in the plane, obtained after solving the corresponding system of equations. For the system C0-C1, we have

$$ \left \{ \textstyle\begin{array}{l} a_{1} = 6 \bigg( \dfrac{2 (\mu -1)}{(2 \mu -1) (2 \mu (\mu _{s}-1)-\mu _{s})}\\ \hphantom{a_{1} =} -2 \mu \mu _{s}-\dfrac{2 \mu }{(2 \mu -1) (2 \mu (\mu _{s}-1)-\mu _{s}+2)}\\ \hphantom{a_{1} =}+2 ( \mu +1)+\mu _{s}\bigg) \\ [3ex] a_{2} = \dfrac{6 \left (-12 \mu ^{2} (\mu _{s}-1)^{2}+12 \mu (\mu _{s}-1)^{2}+\mu _{s} (8-3 \mu _{s})-4\right )}{(2 \mu (\mu _{s}-1)-\mu _{s}) (2 \mu (\mu _{s}-1)-\mu _{s}+2)}. \end{array}\displaystyle \right . $$
(7)

Note that this expression is not defined for \(\mu = 0.5\). This is because for this case, the dipole-segment has a symmetric mass distribution, and therefore any symmetric density (\(a_{1} + a_{2} = 0\)) satisfies conditions C0 and C1. On the other hand, when solving the C1-C2 system, we obtain

$$ \left \{ \textstyle\begin{array}{l} a_{1} = -12 (\mu -3) (\mu _{s}-1) \\ [1.5ex] a_{2} = -30 (\mu _{s}-1). \end{array}\displaystyle \right . $$
(8)

However, we are only interested in solutions that take values within \(\Omega \). This means that not every mass distribution of the ADS is fordable by the parameters of the VDS using the aforementioned properties. In Fig. 4, we illustrate this situation. We can see that in \(\Psi \), there are two sub-regions corresponding to the parameter values resulting in a feasible density (\(a_{1}, a_{2} \in \Omega \)) when solving the systems C0-C1 or C1-C2. By representing the solutions obtained from these regions, we see that \(\Omega \) is almost entirely covered. That is, the variable density model is almost always satisfying C0-C1 or C1-C2 for some mass distribution in the ADS.

Fig. 4
figure 4

Representation of parameters \((\mu , \mu _{s})\) that can be fitted by the variable density model. The blue region on the left corresponds to the parameters in \(\Psi \) that result in a positive density after solving the C0-C1 system, and the orange region with the C1-C2 system. On the right, we can see the respective solution in \(\Omega \) \((a_{1}, a_{2})\)

Full size image

3 Numerical computation of periodic orbits

As explained in the previous section, we chose a rather simple mathematical model of the potential function to illustrate relevant structures of the phase space and gain understanding of the dynamic behavior surrounding these elongated objects. Due to the time-invariant nature of the dynamical problem, the relevant structures are equilibrium points and periodic orbits. Periodic orbits are part of larger families, which can be employed to investigate the structure of the dynamical system. The examination of equilibrium points and periodic orbits, along with their stability as solutions, helps to understand the motion of a small object (orbiter) around an elongated asteroid. To identify the families, we employ various methods for computing and continuing periodic orbits.

3.1 Methodology

In this work, we computed the natural families of periodic orbits for variations of the energy integral taken as a parameter. As it was explained before, we apply different techniques. Typically, asteroid potential models involve the use of spherical harmonics, which result in asymmetric potential functions, leading to three-dimensional periodic orbits without symmetries. In consequence, it is not possible to simplify the calculation of periodic orbits using the symmetry of the system. This fact has led to the development of simpler dynamical models that retain a symmetry. Both models, introduced in Sect. 2, have symmetry with respect to the \(Oy\) axis, that allow us to use a symmetric procedure: the Grid-search method. This method, detailed in 3.1.1, provide us with a map of planar symmetric orbits. These orbits are the initial conditions for the continuation method (see 3.1.2), which aims to evaluate their stability and investigate new bifurcating families. In the following paragraphs we briefly describe the numerical techniques used in the search and characterization of families of periodic orbits.

3.1.1 The grid-search method

Based on the symmetry of the system equation \((x, y, \dot{x}, \dot{y}; t) \mapsto (x, -y, -\dot{x}, \dot{y}; -t)\), it is easy to compute periodic orbits that are symmetric with respect to the \(Ox\)-axis Markellos et al. (1974). The grid-search method applies to two-dimensions Hamiltonian formulations, so it can be used in our dynamical problem to search for planar symmetric orbits. Due to the existence of the Jacobi constant \(C\), Eq. (2), we built a grid on the plane: \((x_{i}, C_{i})\), \(i\in \mathbb{N}\). For each pair \((x_{i},C_{i})\) we have a set of initial conditions at \(t_{0} = 0\): \((x_{0},0,0,\dot{y}_{0}(C))\) that is integrated checking perpendicular crossings with the \(Ox\)-axis. When two crossings are found then it will represent a periodic symmetric and planar orbit. Details of the implementation can be reached in Elipe et al. (2007), Barrio and Blesa (2009). This method has been widely applied to study the families of symmetric periodic orbits of the Three Body Problem in terms of the parameter \(\mu \) that rules the dynamics. Following the similarities of the dipole-segment model with the Restricted Three Body Problem (Elipe et al. 2021), in the following we apply the grid-search technique to search initial conditions of symmetric orbits on the equatorial plane of the asteroid.

3.1.2 Deprit-Henrard continuation method

Natural parameter continuation is based in the smooth evolution of the parameters in the family of periodic orbits. The method proposed is a predictor-corrector (pseudo arc-length) continuation method. It is based on the Deprit and Henrard approach Deprit and Henrard (1967) and later expanded upon by Lara and Peláez (2002), and it is designed to continue families of periodic orbits in conservative systems. First, a periodic orbit is computed using a simple differential corrections process. Then, one of the natural parameters is changed by a small amount. Using the newly converged periodic orbit as the new guess, the process is repeated to evolve new members of the periodic family. The Deprit-Henrard approach follows a predictor-corrector scheme but based in the formulation of variational equations in intrinsic variables to separate normal and tangential displacements, and therefore, avoid trivial eigenvalues that spoils the continuation when working with Hamiltonian systems. The implemented code is a three-dimensional vectorial formulation of the previous scheme whose details can be found in Tresaco et al. (2013).

3.1.3 Stability analysis

The stability of periodic orbits provides information about the motion in their neighborhood. It is possible to estimate if the motion will remain in the vicinity of an orbit for some time, or it will rapidly depart on an unstable trajectory. This fact is retained even for the computation using simplified dynamical models, see Scheeres et al. (2000a). Linear stability is computed using the eigenvalues of the monodromy matrix, which is defined as the solution to the state transition matrix for one period. It is known that for Hamiltonian time-invariant system, the six eigenvalues of monodromy matrix appear in reciprocal pairs. For periodic orbits, a pair of trivial eigenvalues is always equal to one. Hence, the behaviour of a periodic orbit is given by the eigenvalues: \(\{1,1,\lambda _{1} \lambda _{1}^{-1},\lambda _{2},\lambda _{2}^{-1} \}\). The criteria for linear stability is based in the definition of two stability indexes: \(k_{i}=\lambda _{i}+\lambda _{i}^{-1}\) for \(i=1,2\). If both stability indexes satisfy \(|k_{i}|<2\), the orbit is linearly stable, and if one of the index satisfies \(|k_{i}|>2\) it is unstable. The qualitative change in the stability within a periodic family is a bifurcation, which implies that a new family may arise at the bifurcation point.

There are two types of bifurcations during the continuation of the families that predominate in this problem depending on the direction of change of the stability index: tangent bifurcations, when the eigenvalues cross the unit circle at 1, i.e., they jump from the unit circle to the real axis or vice versa and period multiplying bifurcations, when eigenvalues cross at −1.

We also compute the intrinsic stability indexes, \(k_{n}\) and \(k_{b}\), which respectively measure the stability in the normal and binormal directions of the orbit at the initial conditions. This values coincide with the \(k_{1}\) and \(k_{2}\) stability indexes in the case of planar orbits Lara and Peláez (2002).

3.2 Families of periodic orbits

This study investigate the periodic solutions to the Hamiltonian systems described in Sects. 2.1 and 2.2. First, the value of the parameters to carry out numerical experiments is fixed. The characterization of the non-symmetric asteroid for the dipole-segment is given by the values of the parameters: \(\mu =0.3\), \(\mu _{s}=0.5\), and \(l = \omega = \mathcal{G M = 1}\). For the segment with variable density, along Sect. 2.2 we chose the characterization that matches the same ADS configuration with conditions C0-C1: \(a_{1} = -1.95\), \(a_{2} = 0.75\). For these parameters, all the four equilibrium points are linearly unstable for both models. Furthermore, as shown previously, the VDS has a similar qualitative behavior regarding families of periodic orbits.

In the following, we apply the numerical methods of orbital continuations explained in Sect. 3 to compute relevant families of periodic orbits. The outline follows these steps:

  • A grid-search method is applied for a global search for periodic symmetric orbits in the equatorial plane.

  • The periodic orbits encountered in the former step are numerically continued into a family, by varying the energy system. The continuation is carried out by the predictor-corrector scheme.

  • With the numerical continuation, we determine the monodromy matrix for each periodic orbit in the family i.e. the stability of the family, identifying the bifurcation points of these orbits.

These methods require a numerical integrator of differential equations. In our case we have used the Runge-Kutta method RK5(4)7M described in Dormand and Prince (1980), with absolute error tolerance of \(10^{-15}\).

As it can be seen in Fig. 5, we obtain five families of symmetric orbits in the equatorial plane. The notation of these families follows the Strömgren’s notation later used by Henon Hénon (1965) for the analysis of the Restricted Three Body Problem (RTBP).

Fig. 5
figure 5

Families of symmetric periodic orbits. Each point on the plane \((x,C)\) corresponds to a periodic orbit in the equatorial plane

Full size image

We observe three zones: the upper left part, named \(Z_{l}\), that contains families \(l\) and \(ub\); the upper right part, named \(Z_{r}\), that contains families \(k\) and \(va\); and the lower part, named \(Z_{d}\) containing family \(m\). These families are well-know for the restricted three-body problem for the symmetric case, so we follow the same notation for the families.

Families \(l\) (red curve) and \(m\) (green one) correspond to the families obtained taking initial conditions of direct and retrograde circular periodic orbits respectively that start far way from the body. On the other hand, families \(va\) (cyan curve) and \(ub\) (blue curve) are retrograde orbits that encircle the collinear equilibrium \(L_{1}\) and \(L_{2}\) respectively, known as the Lyapunov orbits. These orbits are a composition of known families in the Restricted 3-Body Problem, namely \(u\) and \(b\) on the one hand, and \(v\) and \(a\) on the other Hénon (1965). The merging process of these families is also reported in Elipe et al. (2021), when they study the symmetric dipole model for variation of the parameter \(\mu _{s}\). In addition to this, we obtain family \(k\) (pink family) that encircles the triangular equilibria. Notice that the region from \(x = -0.4\) to \(x = 0.6\) is an unfeasible region since the segment is occupying this part of the \(Ox\) axis.

It is also observed two asymptotic points. Asymptotic or spiral orbits in the RTBP appear when triangular Lagrangian points are unstable, which is the case for the parameters chosen in this paper. These spiral points are the end of families of periodic orbits and represent asymptotic orbits i.e. heteroclinic orbits that connect the two triangular Lagrangian points Henrard (1973). In order to reach the equilibrium point, the asymptotic orbit must have a Jacobi constant equal to the value obtained for the triangular Lagrange points, which is plotted by a black line in Fig. 5. Coordinate \(x\) of the spiral points depends on the values of the parameters that rule the system, see Abad et al. (2023) for the RTBP. Further, the two black points in the figure represent the \((x,C)\) position of the collinear points \(L_{1}\) and \(L_{2}\), and the grey areas are the forbidden regions of motion. Families \(l\) and \(ub\) ends in asymptotic orbits spiralling into the triangular equilibrium \(L_{3}\), and families \(k\) and \(va\) ends in asymptotic orbits spiralling into the other equilibrium \(L_{4}\).

Now, we take advantage of the grid search-method to get a wide set of initial conditions of symmetric periodic orbits to continue. These initial condition of symmetric planar orbits are continued in families using procedures that are not restricted to symmetric orbits, so, families may evolve to non-symmetric orbits, as we will observe in the following analysis. Let us now analyze in detail the behaviour of each family of Fig. 5. Along with the orbit shape of each family we depict the stability evolution of each family orbit depending on the parameter chosen for continuation, the Jacobi constant. By detecting bifurcations, it is possible to identify new periodic families.

3.2.1 Family l

We start with a circular equatorial direct orbit where the satellite is to a distance of \(a=2\ell \). As the Jacobi constant \(C\) tends to \(+\infty \), direct equatorial orbits with large radii are full stable (horizontal and vertically) far from the asteroid. Note that, as happens in this family, it is possible that the continuation of a periodic orbit can always be conducted. In this case, the periodic orbit will converge to a nearly circular planar periodic orbit. Nevertheless, as the family comes closer to the body it becomes unstable. Figure 6 depicts the continuation of the family \(l\) in terms of the Jacobi constant.

Fig. 6
figure 6

Stability evolution along the continuation of family \(l\)

Full size image

We detect three critical points in this family. Two for the horizontal stability index \(k_{n}\), and one for the vertical one \(k_{b}\). Horizontal critical orbits may correspond to an intersection to another family of planar orbits or to an extreme in the Jacobi constant. Both situation occur for this family. For \(k_{n}=2\) there is an extreme of \(C\) that corresponds to a new shape of the orbits Markellos (1977). Near this critic value of \(k_{n}\) the family starts to encircle the triangular equilibrium points. On the other hand, two other critical points close to \(C=3.4\) are observed, one for \(k_{n}=-2\) and other for \(k_{b}=-2\). Horizontal critical orbits, obtained for \(k_{n}=-2\), originates a new planar family of periodic orbits from family \(l\) with doubling-period. These bifurcations are described in Sect. 3.3.

3.2.2 Family \(m\)

Family \(m\) is composed of large quasi-circular retrograde orbits around the asteroid. As it can be observed in Fig. 7, family \(m\) ends in orbits infinitely close to the dipoles. This family is the only one stable for most of its range, establishing a practical region for searching minor bodies orbiting the asteroid. For example, in Lan et al. (2017) it is shown that the moon Dactyl orbits the asteroid Ida in a quasi-periodic motion close to the family \(m\).

Fig. 7
figure 7

Circular retrograde orbits of family \(m\) and the evolution of its stability

Full size image

The Jacobi constant tends to \(-\infty \), and orbits with large radii are full stable. As the family comes very close to the body they become unstable, but note that there is a large interval where orbits are stable in plane but not in space. This fact highlights the importance of considering vertical stability, otherwise a slight perturbation in the out-of-plane component will make the particle to move away from the planar periodic orbit in a spiralling movement away the asteroid.

3.2.3 Families \(ub\) and \(va\)

The families \(va\) and \(ub\) are retrograde orbits that encircle the collinear equilibrium \(L_{1}\) and \(L_{2}\) respectively. \(L_{1}\) and \(L_{2}\) are center-center-saddle equilibrium points, therefore, by Lyapunov’s center manifold Theorem, two families of periodic orbits arise from these equilibrium points. These families are called Lyapunov’s orbits. In the planar family studied in this section, as we increase the radius of the Lyapunov orbits, they become more kidney-bean shaped, see Fig. 8. Families of Lyapunov orbits also exist in the vicinity of the equilateral points Lan et al. (2017). If we observe the evolution of the families, it can be observed that the most external orbits also encircle the triangular equilibria. Note that these planar Lyapunov orbits around the collinear equilibria, are also the termination orbits of a bifurcated family of the Vertical Lyapunov orbit, see Fig. 25.

Fig. 8
figure 8

Evolution of the orbits of the Lyapunov families \(ub\) and \(va\) emanating respectively from \(L_{1}\) and \(L_{2}\)

Full size image

It is found that the planar Lyapunov families contain certain horizontally stable orbits, but all the family is almost vertically stable while normally unstable, see Fig. 9. Three-dimensional Halo orbits also exist which originate from a tangent bifurcations of the Lyapunov families, and will be discussed in the Sects. 3.3.1 and 3.3.2.

Fig. 9
figure 9

Stability evolution of the Lyapunov family \(ub\) and \(va\)

Full size image

Also observe that the difference between the families \(ub\) and \(va\) that can be seen in Fig. 9 is due to the asymmetry of the model considered for the asteroid, \(L_{2}\) is influenced by the bigger dipole but \(L_{1}\) suffers the effect of the smaller dipole.

3.2.4 Family \(k\)

Finally, family \(k\) is composed of direct symmetric periodic orbits that encircle the asteroid. In Fig. 10 orbits belonging to family \(k\) and its stability are represented. This family converges asymptotically to the same heteroclinic orbit as family \(va\). It presents bifurcation points that are analyzed in Sect. 3.3.

Fig. 10
figure 10

Evolution of the orbits of family \(k\) and their stability

Full size image

3.3 Bifurcations of planar families

In this section the bifurcation of the former planar orbits are studied. New bifurcated families are presented along with their stability. The bifurcations are organized depending on the character of the critical orbit: horizontal (normal) and vertical (binormal) bifurcations of simple and double period. Table 1 and Fig. 11 summarize the bifurcation points analyzed that led to new families of periodic orbits.

Fig. 11
figure 11

Representation of the bifurcation points in the plane \((x,C)\)

Full size image
Table 1 Bifurcation points obtained for families of planar symmetric orbits
Full size table

3.3.1 Normal bifurcations

In this subsection, bifurcations due to \(k_{n}=2\) are analyzed. A horizontal critical point \(k_{n}=2\) occurred in the evolution of families \(l\) and \(ub\). As it was pointed in Henon (1967), horizontal critical orbits may correspond to an extrema in \(C\) or to an intersection with another family of planar orbits.

First we have the bifurcation at \(k_{n} =2\) of the family \(l\) (see Table 1). If we apply a perturbation to the initial condition of this orbit in the normal direction, we find a family of asymmetric periodic orbits intersecting the family \(l\), as shown in Fig. 12. This family only appears for a short interval of \(C\), and has a period similar to the bifurcating orbit of \(l\).

Fig. 12
figure 12

Asymmetric family bifurcated out of family \(l\)

Full size image

On the other hand, in the stability diagram of family \(ub\) depicted in Fig. 9, it is observed an horizontal critical point \(k_{n}=2\). A new family of planar non-symmetric orbit emerged through a \(k_{n}=2\) critical point. This bifurcated branch exhibits a vertically stable family of non-symmetric periodic orbits with respect to the \(Ox\) axis, orbits of this family are shown in Fig. 13.

Fig. 13
figure 13

Asymmetric family bifurcated out of Lyapunov family \(ub\)

Full size image

3.3.2 Binormal bifurcations

In this subsection, bifurcations due to \(k_{b}=2\) are analyzed. This critical value indicates an intersection with another family of vertical orbits of single period. A vertical critical point \(k_{b}=2\) occurred in the evolution of all the families presented in 5.

In the stability diagram of family \(m\), Fig. 7, we had a transversal bifurcation \(k_{b}=2\). The orbit shape of this 3-D family is shown in Fig. 14. As the family evolves, one of the extreme of this family gets closer to the dipole segment for increments of the Jacobi constant, and the family seems to terminate in collision orbits.

Fig. 14
figure 14

Orbits and stability of the vertical family, bifurcated out of family \(m\)

Full size image

The same situation occurs in family \(l\). However, in this case the bifurcated vertical family is unstable throughout its entire range, as seen in Fig. 15.

Fig. 15
figure 15

Simple vertical bifurcation of family \(l\)

Full size image

Additionally, the stability diagrams of Lyapunov horizontal families \(ub\) and \(va\) showed critical vertical bifurcations (\(k_{b}=2\)), see Fig. 9. These vertical orbits provided another example of periodic orbits in the vicinity of the collinear equilibria: the Halo orbits. These three-dimensional orbits were first introduced in Farquhar and Kamel (1973) and they have served as the basis for spacecraft mission trajectories. Halo orbits bifurcate from the family of Lyapunov orbits. To compute the Halo orbit from the planar Lyapunov one, a perturbation from the initial condition is taken in the \(z\)-direction. Depending on the perturbation sign, two families of periodic orbits are obtained in the neighborhood of each libration point, labeled as northern and southern. The northern Halo families, presented in Figs. 16 and 17, extend for positive \(z\), while the southern Halo families expand for the negative \(z\) direction, and they are symmetric with respect the \(xy\)-plane to their norther counterparts.

Fig. 16
figure 16

Northern Halo orbits obtained from a bifurcation in family \(ub\) with \(k_{b} = 2\)

Full size image
Fig. 17
figure 17

Northern Halo orbits obtained from a bifurcation in family \(va\) with \(k_{b} = 2\)

Full size image

Finally, family \(k\) also presented a vertical critical orbit for \(k_{b}=2\). This out-of-plane bifurcation gives birth to a new three-dimensional family of periodic orbits, whose shape is shown in Fig. 18.

Fig. 18
figure 18

Three-dimensional family bifurcated out of family \(k\)

Full size image

3.3.3 Doubling period normal bifurcations

We observed during continuation that four of the five families of planar symmetric periodic orbits: \(l\), \(k\), \(ub\) and \(va\), have bifurcation points when \(k_{n}=-2\). This value implies a bifurcation with another planar orbits of double period. Multiple-period orbits are specific periodic orbits interesting for mission design purpose due to their repeatable orbit pattern. In this section we want to highlight the computation of periodic orbits with doubling period by means of the application of the grid-search method.

The grid search method looks for symmetric periodic orbits. As it is explained in Sect. 3.1.1, if the initial conditions of an orbit at \(t=0\) is \((x_{0},0,0,\dot{y}_{0})\) then at half period \(t=T/2\) the orbit conditions will be \((x_{T/2},0,0,\dot{y}_{T/2})\). Then, to find a symmetric orbit, it is checked if the first perpendicular crossing of the orbit with the \(Ox\)-axis satisfies \(y=0\) and \(\dot{x}=0\) (simple orbits). This is how Fig. 5 was obtained. If we are interested in orbits that are close to resonance \(1:2\), we may search for the second perpendicular crossing to the \(Ox\)-axis. Orbits in resonance \(1:2\) can be found inside these bifurcated families of double period. Applying the grid-search procedure to find symmetric planar orbits with multiplicity 2, the map depicted in Fig. 19 is obtained. It presents a map of initial conditions \((x_{i},C)\) of the families of multiplicity 1 (black curves) and the intersection with new families of multiplicity 2 (coloured ones).

Fig. 19
figure 19

Families of multiplicity 1 (black curves) and their corresponding period-doubling branches in zones \(Z_{l}\) and \(Z_{d}\). Grid \((x,C)\)

Full size image

Note that in Fig. 19 we have only presented the families that have the second crossing to \(Ox\)-axis perpendicular and intersect to the planar families of multiplicity 1; there exist more families of multiplicity 2, but along this paper we are only interested in the analysis of bifurcations of planar symmetric orbits. The intersection points occur for the simple-period families at stability index value \(k_{n}=-2\), which means that the coloured curves are the period-doubling bifurcated families of the simpler ones: \(l\), \(k\), \(ub\) and \(va\). The intersections among horizontal (and vertical) critical orbits with families of planar ones were first reported in Henon (1967).

Let us analyze each family separately. We first highlight the family in red color. This family is composed of the doubling-period bifurcated orbits of family \(l\). The evolution of this family is shown in Fig. 20. This family starts at the critical point, \(k_{n}=-2\), and is composed of circular orbits around the dipole, and its evolution include orbits that encircle the triangular equilibrium with large period and high instability.

Fig. 20
figure 20

Orbit evolution and stability diagram of planar family from family \(l\) with doubling period

Full size image

Blue curve in Fig. 19 corresponds to period-doubling orbits for family \(ub\), both families collide through a period-doubling bifurcation. This doubling-period family evolves to collision orbits for decreasing values of the Jacobi constant, see Fig. 21. Similarly to family \(ub\) in \(Z_{l}\), the family obtained for multiplicity 2 coloured in cyan in Fig. 19 bifurcated out of family \(va\) (Fig. 8) by a period doubling phenomena. Its evolution is depicted in Fig. 22.

Fig. 21
figure 21

Planar doubling-period family bifurcated from \(ub\) and its stability

Full size image
Fig. 22
figure 22

Planar doubling-period family bifurcated from \(va\) and its stability

Full size image

Finally, the magenta family of multiplicity 2 is composed of the doubling-period bifurcated orbits of family \(k\), see Fig. 23. This family evolves to orbits that approaches to the triangular libration points with high instability, and it terminates in collision orbits.

Fig. 23
figure 23

Continuation of the normal family emanating from \(k\) and its stability

Full size image

Note that family \(m\) does not present a planar period doubling bifurcation. Its stability diagram, see Fig. 7, showed that the horizontal stability index \(k_{n}\) is tangent to the value −2, so we did not expect any family of multiplicity 2 in this zone.

3.3.4 Doubling period binormal bifurcations

In this subsection, bifurcations due to \(k_{b}=-2\) are analyzed. This critical value indicates an intersection with another family of vertical orbits of double period. A vertical critical point \(k_{b}=-2\) occurred in the evolution of family \(l\) and family \(m\).

Family \(l\) presents a vertical critical orbit of this type. For the out-of-plane stability index \(k_{b}=-2\) observed in Fig. 6 we have a family of three-dimensional periodic orbits which intersects the family of planar orbits twice. Its trajectory and stability diagram is presented in Fig. 24.

Fig. 24
figure 24

Vertical bifurcated family from family \(l\) with doubling period

Full size image

Family \(m\) presented a bifurcation of this type, and emerges a three-dimensional family of periodic orbits which intersects the family of planar orbits twice. This family is composed of eight-shaped orbits near a collinear equlibrium point, known as Vertical Lyapunov orbits, see Fig. 25. This family is unstable and shrinks in size until its termination at one of the collinear points.

Fig. 25
figure 25

\(L_{1}\) vertical Lyapunov family and its stability evolution

Full size image

Although we have not systematically followed bifurcations due to \(|k_{i}|=2\) in the bifurcated families analyzed in this Sect. 3.3, we highlight a bifurcated case of this unstable \(L_{1}\) vertical family. In Fig. 25 it is observed a bifurcation point for stability index value \(k_{2}=2\). A new family is originated. This family is interesting because it is composed of unstable eight-shaped orbits that are twisting for increments of the Jacobi constant until its termination in the Horizontal Lyapunov orbit. Members of the vertical family appear in Fig. 26, observe that the planar orbit corresponds to the planar Lyapunov orbit given in Fig. 8. Remind that, according to Lyapunov’s center theorem, there are two families of periodic orbits that emanate from the collinear equilibrium points, the planar Lyapunov orbit and the vertical Lyapunvov orbit. Lyapunov vertical orbit connects to the planar Lyapunov orbit and to the circular retrograde plane orbit exterior to the equilibrium points, which is of family \(m\). These bridges connecting the two Lyapunov families are well know in literature Gómez and Mondelo (2001).

Fig. 26
figure 26

Orbit evolution of the bifurcated family of \(L_{1}\) Vertical Lyapunov family

Full size image

4 Conclusions

This paper aimed to the search for periodic orbits around irregular asteroids using simplified models of the gravitational potential. In particular, this research studies elongated asteroids with asymmetric density profile by means of two different modelling approaches. The first one, the dipole-segment model, describes the gravitational potential of elongates bodies focusing on the geometry of the different density components. On the other hand, the variable density segment model is proposed in this paper as an alternative to the dipole-segment by implementing a continuous parametrization of the internal density of the asteroid.

Both dynamical models are presented and studied in terms of the equations of motion, zero-velocity curves and equilibrium points. Using a second degree polynomial allows the variable density model to fit mass distribution properties of the dipole-segment, increasing the resemblance between models. In such scenario, there is no qualitative differences in the computation of periodic trajectories and stability areas, ensuring the applicability of both models to the orbital dynamic study around elongated asteroids. Nevertheless, the variable density model could be extended implementing a higher order density function if the irregular asteroid presents more lobes with different mass concentrations.

Along the paper we have shown a wide set of periodic orbits. The symmetric, planar and 1-multiplicity orbits served as starting point for the computation of a large number of new families emanating from the bifurcation process. Most of these orbits are unstable, except periodic orbits near the equatorial plane and, in particular, resonant quasi-planar orbits that are full stable. The analysis of the stability indexes also reveals that multiple-period planar bifurcations can be obtained searching for symmetric orbits with multiplicity greater than one in the grid-search method. When searching for multiplicity 2 orbits, we encountered four families of periodic orbits corresponding to doubling-period bifurcations of families \(l\), \(k\), \(ub\) and \(va\). This procedure can be repeated to compute families of multiplicity greater than 2, and also to obtain resonant orbits. This would be tackled in an oncoming work. In addition, we have presented three-dimensional families originated as a superposition of a planar periodic motion and a vertical impulse. For instance, Halo orbits emerged as bifurcation from the planar Lyapunov orbits, being a good desirable trajectory to study the asteroid if the orbit is selected with a large enough vertical amplitude. As it can be observed in the stability diagrams of vertical families, their evolution reveal new critical points corresponding to new bifurcated three-dimensional families that occur at \(k_{i}=-2\). Then, these vertical critical orbits can be used again as starting points for finding new families of three-dimensional periodic orbits. The sample orbits listed along the paper are not all the periodic orbits around elongated asteroids, but are an example of the variety of orbits that can be followed for a particle in orbital motion about it, and help us to a better understanding of the dynamics around these elongated bodies.