© The Authors 2026

1 Introduction

Over the past few decades, space agencies have increasingly directed their efforts toward exploring the small bodies that populate our Solar System by dispatching spacecraft to study asteroids and comets, with further missions planned for the future. The shapes of asteroids and comets are often irregular. This poses challenges in gravitational modeling. When space missions are designed to explore these celestial bodies, a critical question is how the gravitational field of these small irregularly shaped objects can be represented accurately (Elipe & Lara 2003; Bartczak & Breiter 2003). To address this challenge, numerous mathematical models have been proposed that each offer potential solutions to better understand and characterize the gravitational effect of these bodies.

Spherical harmonic expansion is commonly employed to model planetary bodies because these objects typically exhibit a nearly spherical shape (Elipe & Lara 2003). This method becomes less effective when it is applied to objects with irregular geometries, however, because the harmonic expansion might no longer be suitable and convergence issues might arise (Riaguas et al. 1999; Jiang & Baoyin 2018; Lan et al. 2017). Traditionally, the gravitational fields of irregular or non-homogeneous mass distributions are modeled using spherical harmonic expansions. However, these series often suffer from high computational costs and poor convergence. These limitations have created a significant bottleneck in current modeling techniques, necessitating the exploration of alternative solutions to the Laplace equation, through the development of closed-form expressions for the potential function (Broucke 1999).

A widely used alternative for representing the gravitational field of irregularly shaped bodies is the polyhedron method (Werner 1994; Werner & Scheeres 1996; Werner 1997; Scheeres et al. 2000; Tsoulis & Petrovi 2001; Tsoulis 2012). This model is regarded as highly accurate and does not present convergence challenges (Scheeres et al. 2000). The polyhedron model has been extensively used in real mission projects to study spacecraft dynamics near asteroids because it is precise. A significant drawback of the polyhedron model is the need for thousands of parameters (vertices and facets) to maintain its high level of accuracy, however. This in turn results in substantial computational demands and extended simulation times. Additionally, the analysis of the general dynamics of spacecraft with respect to the model parameters becomes increasingly complex because the parameters are interrelated and exert a combined effect on the effective gravitational field of irregular bodies.

In light of these limitations, simplified models provide an appealing alternative for studying particle dynamics around complex bodies. They offer a more computationally efficient approach than more elaborate models such as the polyhedron model (Werner 1994) or the mascon model (Geissler et al. 1996; Werner & Scheeres 1996; Amarante et al. 2021). By reducing computational complexity while maintaining reasonable accuracy, simplified models are valuable tools for investigating the gravitational environment around irregularly shaped celestial bodies. Simplified models provide an effective approach for approximating the gravitational potential of irregular celestial bodies. These models are particularly useful in celestial mechanics because they help us to overcome significant challenges related to calculating gravitational fields for bodies with nonuniform mass distributions or irregular shapes, such as decreasing time-step numerical integration costs and/or simplifying analytical analyses. Simplified models can often provide a straightforward analytical expression for calculating the gravitational field of irregular bodies (Yang et al. 2018). An additional advantage of using these models is the ease with which the effects of specific parameters on the dynamics around asteroids can be analyzed. This includes, for example, the stability of equilibrium points (Zeng et al. 2015; Barbosa Torres dos Santos et al. 2017), the distribution of stable periodic orbits (Lan et al. 2017), and the identification of possible hovering regions (Yang et al. 2015; Zeng et al. 2016c). Moreover, simplified models are beneficial in orbit design (Wang et al. 2017; Barbosa Torres dos Santos et al. 2017) and feedback control strategies (Yang et al. 2017). While specific trajectories can be more accurately analyzed with refined models during the final stages of a mission, a general understanding of the gravitational field through a limited set of parameters proves to be invaluable in the early phases of mission planning.

Several simplified models have been introduced to represent irregular asteroids, as discussed by Riaguas et al. (1999, 2001); Bartczak & Breiter (2003); Elipe & Lara (2003), where the dynamics of a particle under the gravitational effect of an asteroid, modeled as a straight segment, were explored. Other studies investigated the motion of a particle around small celestial bodies using simplified models through a simple planar plate (Blesa 2006), a homogeneous cube (Liu et al. 2011b), a triaxial ellipsoid (Gabern et al. 2006), a rotating homogeneous cube (Liu et al. 2011a), double and finite straight segments Jain & Sinha (2014), and a rotating-mass dipole (Kirpichnikov & Kokoriev 1988; Kokoriev & Kirpichnikov 1988). A method for adapting the rotating-mass dipole model to a real asteroid was proposed for the first time by Zeng et al. (2015). This model was subsequently refined by considering the oblateness of the primary bodies (Zeng et al. 2016a,b, 2017) and incorporating the dipolesegment model (Zeng et al. 2018). From then on, studies for the dipole-segment model were presented by Elipe et al. (2021a); Abad et al. (2024) with the aim of analyzing the dynamics around irregular asteroids. In the context of elongated asteroid binary systems, Barbosa Torres dos Santos et al. (2017) applied the restricted synchronous three-body problem model.

Inspired by Zeng et al. (2015), Lan et al. (2017) introduced the rotating-mass tripole model with symmetrical rotation. The authors proposed that small convex bodies are effectively modeled using the mass-tripole model. They demonstrated that by using five parameters (determined with the aid of the polyhedron model), the geometric configuration can be defined, the physical characteristics of a real asteroid can be derived, and its gravitational field can be computed. To extend the work of Lan et al. (2017), Barbosa Torres dos Santos et al. (2020) conducted a semi-analytical investigation into the qualitative dynamics around an asteroid with a convex shape and employed the rotating-tripole model. Santos et al. (2021) proposed a simplified model that considered the spatial distribution of the irregular body instead of a distribution in the xy-plane to describe the gravitational fields of elongated asteroids. This model provided results with a better accuracy than existing axisymmetric models and the mass-point model at low computational cost. Also for a tripolar system, Boaventura et al. (2023) analyzed the dynamics around a simplified tripolar-system model with a segment, based on the study of symmetric periodic orbits, applied to asteroid Holiday. These studies and others demonstrated the value of using simplified models to identify the key parameters that affect the motion of particles in specific asteroid systems.

Numerous models have been proposed to study the gravitational environment around irregular asteroids, with a focus on simplifying the complex dynamics while maintaining a reasonable degree of accuracy. Among these, the dipole-segment model is a commonly used approach for approximating the gravitational field of elongated bodies (e.g., Zhang et al. 2018; Zeng et al. 2018; Elipe et al. 2021b).

To address the challenges associated with the gravitational modeling of elongated irregular bodies, we propose the generalized dipole-segment model (GDSM). This model significantly extends the traditional dipole-segment approach by incorporating five key parameters: the force ratio (k), the pole masses (μ), the connecting rod mass (μ_s), and the flattening coefficients of the poles (A₁ and A₂). To identify these parameters for our proposed model, we employed nonlinear optimization techniques, which are detailed below. The advantages of the model are outlined below.

By advancing beyond traditional dipole-segment models, which simplify the poles as point masses, the proposed approach introduces distinct pole masses and a massive connecting rod and considers the spheroidal shape of the poles. These refinements represent a meaningful step toward a more realistic representation of the gravitational field of elongated bodies.

This paper is structured as follows. Section 2 presents the normalized equation of motion for a particle in the vicinity of elongated asteroids, which are modeled using the Generalized Dipole-Segment model. Section 3 outlines the method we employed to determine the parameters of the simplified model through a nonlinear optimization approach. In Sect. 4, the proposed method and model are applied to two irregular asteroids, Arrokoth and Kleopatra, and to comet 103P/Hartley. Section 5 compares the Generalized Dipole-Segment model with the Dipole-Segment model by evaluating them against the polyhedron model. This comparison highlights the advantages of the newly proposed model. Section 6 presents the numerical simulations, in which the equilibrium points computed from the GDSM are used as reference to construct the heteroclinic orbits around asteroids Arrokoth, Kleopatra, and 103P/Hartley. Finally, Sect. 7 summarizes our conclusions.

2 Dynamical equations of the generalized dipole segment

We considered an inertial reference frame denoted as Oxyz. The body was composed of two masses with a spheroid format m₁ and m₂, connected by a segment with a mass m₃ and length l, with a constant linear mass density ρ. The total mass of the system is M = m₁ + m₂ + m₃. The masses m₁ and m₂ are assumed to be oblate bodies, meaning that they exhibit a flattened shape. Depending on the value of the oblateness coefficient (A), these bodies can either be oblate (flattened at the poles when the coefficient is positive) or prolate (elongated at the poles when the coefficient is negative). The parameter A₁ represents the oblateness coefficient of the body with a mass m₁, and A₂ denotes the oblateness coefficient of the body with a mass m₂.

This body rotates uniformly in the xy-plane with a constant angular velocity ω around its center of mass, denoted by O. We refer to this configuration as the Generalized Dipole-Segment (GDS).

The GDS system rotates, and we therefore introduced a synodic reference frame Oxyz such that the axes OZ and Oz were aligned and the body continuously rotated about the Ox-axis.

We introduced two mass ratios, $$\mu =\frac{m_2}{m_1+m_2},{\mu }_s=\frac{m_3}{m_1+m_2+m_3},$$(1)

which assumed values within the intervals μ ∊ [0,1] and μ_s ∊ [0,1], respectively.

Let (-l₁, 0, 0) and (l₂, 0, 0), where l = l₁+l₂ = 1 in canonical units (c.u.) represent the positions of the masses m₁ and m₂ in the Oxyz reference frame, respectively. l₁ and l₂ are defined as presented in Eqs. (2) and (3), $$l_1=\mu (1-{\mu }_s)+\frac{{\mu }_s}{2},$$(2) $$l_2=(1-\mu )(1-{\mu }_s)+\frac{{\mu }_s}{2}.$$(3)

When the mass ratio μ_s is precisely zero, the system simplifies to two basic scenarios: either a single-point mass, or a conventional dipole configuration. Conversely, when μ _s equals one, the model corresponds to the standard segment configuration. In the canonical unit system, the masses of bodies m₁, m₂, and m₃ are m₁ = (1 - μ)(1 - μ_s), m₂ = μ(1 - μ_s), and m₃ = μ_s, respectively.

We considered the Generalized Dipole-Segment model to rotate with a constant angular velocity. The vector equations governing the motion of a test particle affected by the gravitational field of the Generalized Dipole-Segment model are expressed in Eq. (4), $$\ddot{\mathbf{r}}+2\mathit{\omega }\times \dot{\mathbf{r}}=-\mathrm{\nabla }V,$$(4)

in the body-fixed frame Oxyz, as depicted in Fig. 1. The effective potential V is defined in Eq. (5), $$\begin{array}{r}V={\omega }^2\frac{1}{2}(x^2+y^2)+{\omega }^{2}k[\frac{(1-\mu )(1-{\mu }_s)}{r_1}(1+\frac{A_1(r_1^2-3z^2)}{2r_1^4})\\ +\frac{\mu (1-{\mu }_s)}{r_2}(1+\frac{A_2(r_2^2-3z^2)}{2r_2^4})+{\mu }_s\log \left(\frac{1+r_1+r_2}{-1+r_1+r_2}\right)],\end{array}$$(5)

where r₁ and r₂ are given by $$r_1=\sqrt{(x+l_1{)}^2+y^2+z^2},$$(6) $$r_2=\sqrt{(x-l_2{)}^2+y^2+z^2}.$$(7)

A_i (i = 1, 2) represent the oblateness coefficients of the primary bodies, which are defined in Eq. (8), $$A=\frac{({\rho }_i^e{)}^2-({\rho }_i^p{)}^2}{5l^2},$$(8)

where the parameter ρ denotes the radius of the primary spheroids. The subscript e corresponds to the equatorial radius of the primary body, and the subscript p refers to its polar radius (Santos et al. 2017; Abozaid et al. 2024). The dimensionless scalar parameter k is known as the force ratio, which represents the proportion between gravitational acceleration and centrifugal acceleration, given by Eq. (9) (Santos et al. 2021), $$k=\frac{GM}{{\mathrm{\Omega }}^2{l}^3},$$(9)

where Ω and l are the angular velocity and length of the asteroid, respectively, in the international system of units. The system in Eq. (9) for k = 1 corresponds to the restricted three-body problem. When k < 1, the segment rotates rapidly, whereas k >1 indicates a slower rotation. As noted by Prieto-Llanos & Gomez-Tierno (1994), cohesive forces are required in cases of rapid rotation to prevent the body from disintegrating as a result of inertial forces that are not counterbalanced by self-gravitation. This suggests that scenarios with k << 1 are not particularly realistic (Abad et al. 2024).

The gradients of the effective potential in three-dimensional space (x, y, z) are given by Eqs. (10), (11), and (12), $$\ddot{x}-2\omega \dot{y}=V_x,$$(10) $$\ddot{y}+2\omega \dot{x}=V_y,$$(11) $$\ddot{z}=V_z,$$(12)

where V_x, V_y, and V_z denote the partial derivatives of the effective potential V with respect to x, y, and z, respectively.

Equation (13) possesses the well-known Jacobi integral, $$v^2=2V-C,$$(13)

which establishes the permissible region for possible motions. Equation (13) indicates that for a given value of C, the velocity v depends on the position of the body within the plane of motion (Barbosa Torres dos Santos et al. 2017). The integration constant C is determined by the initial position and velocity of the particle, as stated by McCuskey (1963). Furthermore, Eq. (13) establishes a connection between the squared velocity and the

coordinates of the body with negligible mass in the rotating reference frame, as mentioned by Moulton (1970). Consequently, when the integration constant C is numerically obtained from the initial conditions, Eq. (13) allows us to compute the velocity of the body at any point in space. Conversely, for a given velocity, Eq. (13) defines the geometric regions of permissible positions for the body. Notably, when the velocity is set to zero in Eq. (13), the resulting region corresponds to locations in which the velocity of the particle is also zero, as discussed by Szebehely (1967). Mathematically, this boundary is described by the equation 2V - C = 0, as presented by McCuskey (1963).

Fig. 1 Geometric configuration of the rotating coordinate system for the generalized dipole-segment model, where m1 and m2 represent the masses at the poles, and m3 corresponds to the mass of the central segment.

3 Determining the parameters for the generalized dipole-segment model

The Generalized Dipole-Segment model was presented in the previous section. This model has unknown parameters that need to be determined to provide a complete set of equations for the system. In the simplified model we considered, the total mass and angular velocity of the asteroids are equivalent to the actual values of these asteroids. The Generalized Dipole-Segment model therefore has five unknown parameters, which are the mass ratio (μ), the mass of the segment (μ_s), the spheroidal shape of the poles (A₁ and A₂), and the force ratio (k).

The unknown parameters were obtained using the locations of the external equilibrium points calculated through the polyhedral approach, following a method similar to that used by Zeng et al. (2015, 2016c); Yang et al. (2018). The primary objective was to identify the parameters that produce equilibrium points as close as possible to those derived from the polyhedron model. This approach is suitable for asteroids with external equilibrium points. Certain asteroids lack external equilibrium points, however, such as asteroid 1998 KY26. In these cases, an alternative method that minimizes the errors of the effective potential or its gradient can be applied, as suggested by Zeng et al. (2018). The coordinates of the equilibrium points [◯_E, ŷ_E, ẑ_E] were calculated as indicated in Eq. (14), $$V_x({\hat{x}}_E,{\hat{y}}_E,{\hat{z}}_E)=V_y({\hat{x}}_E,{\hat{y}}_E,{\hat{z}}_E)=V_z({\hat{x}}_E,{\hat{y}}_E,{\hat{z}}_E)=0.$$(14)

The objective was to identify the parameters of the Generalized Dipole-Segment model that minimize the discrepancies between the external equilibrium points of the simplified model and those of the polyhedron model. The procedure we used to determine the parameters for the Generalized Dipole-Segment model is outlined below.

The optimization variables for this model were X = [μ, μ_s, A₁, A₂, k]. Before we performed the optimization, constraints were set for each variable using the lower bounds [μ_min, μ_{s min}, A_{1 min}, A_2min, k_min] and the upper bounds [μ_max, μ_{s max}, A_1max, A_2min, k_max] for these parameters. For the nonlinear optimization problem, the performance index is defined as $$\begin{array}{r}\mathbf{J}(\mu ,{\mu }_s,k,A_1,A_2)=\\ \sum _{i=1}^n\sqrt{({\hat{x}}_{G_i}{d}^{\ast }-x_{P_i}{)}^2+({\hat{y}}_{G_i}{d}^{\ast }-y_{P_i}{)}^2+({\hat{z}}_{G_i}{d}^{\ast }-z_{P_i}{)}^2},\end{array}$$(15)

where [◯_Gi,ŷ_Gi, ẑ_Gi] denotes the positions of the equilibrium points obtained from the simplified models (Generalized DipoleSegment), expressed in canonical units. It is important to note that l can be calculated using Eq. (9). [x_pi,y_pi, z_pi] represent the coordinates of the equilibrium points (measured in meters) as determined by the polyhedron model. The index i corresponds to the ith equilibrium point, and n denotes the total number of external equilibrium points.

The equality constraints (c_eq = 0) for the Generalized Dipole-Segment model are presented in Eqs. (16) and (17), $$c_{eq}=[m_1+m_2+m_3-1],$$(16) $$c_{eq}=[l_1+l_2-1].$$(17)

Mathematically, the performance index can be expressed as a constrained minimization, as illustrated in Eq. (18), $$min\mathbf{J}(\mu ,{\mu }_s,k,A_1,A_2)suchthat{c}_{eq}=0,$$(18)

where J(μ,μ_s, k, A₁, A₂) is a function that produces a scalar output, and c_eq represents functions that yield vector outputs. Nonlinear optimization routines implemented in Matlab were used to obtain the optimal solutions that minimize J (Santos et al. 2022). The optimization problem was formulated and solved using a nonlinear programming approach (NLP).

4 Applications to realistic elongated asteroids

In this section, we apply the Generalized Dipole-Segment (GDS) model to two realistic elongated asteroids, KBO 486958 Arrokoth and Kleopatra, and to comet 103P/Hartley (commonly known as Hartley 2). Using an optimization method, we determine the parameters of the GDS for each body. The main distinction between the dipole-segment (DS) and the Generalized Dipole-Segment (GDS) models lies in the flattening coefficients A₁ and A₂. When A₁ = A₂ = 0, the model reduces to the classical dipole-segment formulation. Conversely, when A₁ and A₂ are nonzero, the Generalized Dipole-Segment model is obtained, which allows us to represent more elongated and asymmetric mass distributions. We demonstrate the improved performance of the Generalized Dipole-Segment model compared to the dipole-segment model below.

4.1 Parameters of the elongated asteroids

The physical and orbital characteristics of asteroids KBO Arrokoth and Kleopatra and of comet 103P/Hartley were obtained from relevant literature sources and are summarized in Table 1. We also employed the polyhedron model to accurately characterize the gravitational fields of these bodies (for further details of this method, we refer to Werner 1994; Werner & Scheeres 1996; Tsoulis & Petrovi 2001; Tsoulis 2012).

4.2 Equilibrium points

We employed the polyhedron model as described by Werner (1994); Werner & Scheeres (1996); Tsoulis & Petrovi (2001); Tsoulis (2012). Using this approach, we determined the positions of the equilibrium points numerically by solving the condition $$\mathrm{\nabla }V=\mathbf{0}.$$(19)

Using the gravitational potential from polyhedron model (Eq. (22)) together with the MINOR-EQUILIBRIA package¹ (Amarante & Winter 2020; Amarante 2020a), we computed the equilibrium points by solving the force balance defined by Eq. (19).

The results are presented in Table 2. We emphasize that only the external equilibrium points were considered because the internal points have no physical significance.

To perform the simulations, it was essential to establish the boundary constraints for each asteroid we studied using the generalized dipole-segment model. These constraints are defined below.

1. Generalized dipole-segment model

   For each asteroid, [μ_min, μ_max] were set to [0.001, 0.999] kg, [μ_{S, min}, μ_{S, max}] were set to [0.001, 0.999], [k_min, k_max] were set to [0.1, 9], [k_min, k_max] were set to [0, 9], [A_1,min, A_1,max] were set to [-4, 4], and [A_2,min, A_2,max] were set to [-4, 4].

2. Dipole-segment model

   For each asteroid, [μ_min, μ_max] were set to [0.001, 0.999] kg, [μ_{S, min}, μ_{S, max}] were set to [0.001, 0.999], [k_min, k_max] were set to [0.1, 9], [k_min, k_max] were set to [0, 9], [A_1,min, A_1,max] were set to [0, 0], and [A_2,min, A_2,max] were set to [0, 0].

The boundary constraints for the geometric parameters μ and μ_S were determined based on the classical canonical unit system definition, where the mass ratios cannot exceed unity because their sum must satisfy the condition μ + μ_s = 1.

The parameter k was initially defined within the range [0.1,9]. When the optimal solution obtained from the simulations approached the upper limit of this range, the interval was adjusted to [5,15], and the optimization was performed again. This process was repeated iteratively until the optimal value of k was well centered within the defined range.

A similar approach was applied to the flattening coefficients A₁ and A₂, which were initially constrained within the interval [-4,4]. When the optimization results indicated values close to the boundaries, the range was adjusted accordingly. This adaptive procedure ensured that the final parameter values remained well within the specified limits.

The initial guesses for the GDSM and DSM models were chosen as described below.

1. GDSM

   For Arrokoth, we used [μ, μ_s, k, A₁, A₂] = [0.1, 0.8, 5, 0.1, 0.1], for Kleopatra, we used [μ, μ_s, k, A₁, A₂] = [0.3, 0.3, 5, 0.1, 0.1], and for 103P∕Hartley, we used [μ, μ_s, k, A₁, A₂] = [0.6, 0.54, 0.5, 0.1, 0.1].

2. DSM

   For Arrokoth, we used [μ, μ_s, k, A₁, A₂] = [0.8, 0.54, 8, 0, 0], for Kleopatra, we used [μ, μ_s, k, A₁, A₂] = [0.3, 0.3, 5, 0, 0], and for 103P∕Hartley, we used [μ, μ_s, k, A₁, A₂] = [0.6, 0.54, 0.5, 0, 0].

The optimal values of the parameters μ, μ_S, k, A₁, and A₂, determined using the optimization method described above, are presented in Tables 3 and 4. These tables also display the performance indices J_max and J_min. The performance indices were computed after we determined the optimal parameters, denoted J. Specifically, the indices J_max and J_min represent the maximum and minimum relative errors for a single equilibrium point and were defined as follows.

J_max refers to the maximum relative error between the equilibrium points calculated using the simplified model and the actual equilibrium points located within the high-fidelity model (the polyhedron model). For four equilibrium points, the error in the location of each equilibrium point was computed for the four external equilibrium points. The largest of these relative errors was designated J_max, and the smallest was identified as J_min. The parameter J_max and J_min depends on the optimization variables (μ,μ_s, k, A₁, A₂), and it was defined as shown in Eqs. (20) and (21), $$\begin{array}{r}{\mathbf{J}}_{\mathbf{m}\mathbf{a}\mathbf{x}}=max\left(\frac{\sqrt{({\hat{x}}_{G_i}{d}^{\ast }-x_{P_i}{)}^2+({\hat{y}}_{G_i}{d}^{\ast }-y_{P_i}{)}^2+({\hat{z}}_{G_i}{d}^{\ast }-z_{P_i}{)}^2}}{L}\right)\\ \times 100\mathrm{\%},i=1,2,3,4,\end{array}$$(20) $$\begin{array}{r}{\mathbf{J}}_{\mathbf{m}\mathbf{i}\mathbf{n}}=min\left(\frac{\sqrt{({\hat{x}}_{G_i}{d}^{\ast }-x_{P_i}{)}^2+({\hat{y}}_{G_i}{d}^{\ast }-y_{P_i}{)}^2+({\hat{z}}_{G_i}{d}^{\ast }-z_{P_i}{)}^2}}{L}\right)\\ \times 100\mathrm{\%},i=1,2,3,4.\end{array}$$(21)

The results show that the values of J are consistently lower when the GDSM is used, which indicates a more accurate representation of the system. For highly elongated bodies with relatively uniform extremities, that is, without significant variations in shape at their poles, the minimized J values yielded comparable results for both models, however. This is particularly evident for asteroid Kleopatra, where the mass parameter μ_s is comparable to the mass at the extremities.

On the other hand, when the body had a noticeable bulge at one of its extremities and deviated from a nearly uniform elongated shape, the GDSM proved to be significantly more effective than the DSM. In these cases, the values of the flattening coefficients A₁ and A2 are substantial, which reinforces the insight that the Generalized Dynamical Shape Model captures the asymmetry of the system better.

The accuracy of the GDSM becomes even more evident for the mass distribution within the system. When the body presents a pronounced protrusion at its poles, most of the mass is concentrated at these extremities, which are typically asymmetric and not perfectly spherical. This structural characteristic enhances the performance of the GDSM because it captures the gravitational effects arising from these irregularities better.

Conversely, for bodies with a more uniform shape, where the poles do not exhibit significant protrusions, the mass distribution along the central structure and at the extremities tends to be on the same order of magnitude. In other words, the mass of the elongated central region is comparable to the mass that is concentrated at the poles. As a result, the system behaves in a nearly homogeneous manner, which leads to comparable optimized values of J for both models, with J remaining below unity. The DSM then also provides satisfactory results, although the GDSM still proves to be the more precise approach.

A clear example is the case of KBO Arrokoth and comet 103P/Hartley, where the mass parameter μ_s is small, which leads to a mass distribution that is concentrated at the poles. Because these poles are not spherical, the GDSM significantly outperforms the DSM in accurately modeling the gravitational field of the system.

5 Comparison between the simplified and the polyhedron models

The polyhedron approach represents an irregularly shaped body with a uniform density (ρ) as a collection of tetrahedra bounded by triangular facets. Based on this geometric description, the gravitational potential and the gradient of gravitational potential can be expressed as (Tsoulis & Petrovi 2001: Tsoulis & Petrovi 2001) $$\begin{array}{r}U(x_1,x_2,x_3)=-\frac{G\rho }{2}\sum _{p=1}^n{\sigma }_p{h}_p[\sum _{q=1}^m{\sigma }_{pq}{h}_{pq}{L}_{N_{pq}}\\ +h_p\sum _{q=1}^m{\sigma }_{pq}{A}_{N_{pq}}+\sin (g{A}_p)],\end{array}$$(22) $$\begin{array}{r}-\frac{\mathrm{\partial }U(x_1,x_2,x_3)}{\mathrm{\partial }{x}_i}=-G\rho \sum _{p=1}^n\cos ({\mathbf{N}}_p,{\mathbf{e}}_i)[\sum _{q=1}^m{\sigma }_{pq}{h}_{pq}{L}_{N_{pq}}\\ +h_p\sum _{q=1}^m{\sigma }_{pq}{A}_{N_{pq}}+\sin (g{A}_p)],\\ (i=1,2,3).\end{array}$$(23)

Auxiliary functions are given by $$L_{N_{pq}}=\ln \left(\frac{s_{2pq}+l_{2pq}}{s_{1pq}+l_{1pq}}\right),$$(24) $$A_{N_{pq}}=\arctan \left(\frac{h_p{s}_{2pq}}{h_{pq}{l}_{2pq}}\right)-\arctan \left(\frac{h_p{s}_{1pq}}{h_{pq}{l}_{1pq}}\right),$$(25)

where the parameters of Eqs. (22)-(25) were defined by Tsoulis & Petrovi (2001) in Sect. 5.

The results indicate that the equilibrium points computed using the GDSM agree better with those derived from the polyhedron model than the points obtained with the DSM, as shown in Table 2. This constitutes the first evidence that supports that the GDSM approximates the reference model better. Next, we analyzed the classifications of the equilibrium points for the three asteroids we studied. Our analysis is consistent with the analysis conducted using the polyhedron model. We then compared the relative errors between the gradients of the pseudo-potential computed with the GDSM and those obtained from the polyhedron model in order to assess the validity of the simplified formulation. Because the gradient of the pseudopotential is a vector quantity, with components associated with the partial derivatives in x, y, and z, the comparison was performed using the magnitude of this vector at each sampled spatial location. The relative error was then calculated as shown in Eq. (26), $${\epsilon }_i=\left(\frac{\mathrm{\nabla }{V}_i^{\text{GDSM}}-\mathrm{\nabla }{V}_i^{\text{poly}}}{\mathrm{\nabla }{V}_i^{\text{poly}}}\right)\times 100\mathrm{\%},$$(26)

where $\mathrm{\nabla }{V}_i^{\text{GDSM}}$ and $\mathrm{\nabla }{V}_i^{\text{poly}}$ denote the gradients of the pseudopotential evaluated at the same spatial position in the rotating reference frame using the GDSM and the polyhedron model, respectively. The mean and maximum values of ε_i over all sampled points were then analyzed to assess the accuracy of the simplified model. The gradient of the gravitational potential from the polyhedron model was computed with the package MinorGravity² (Amarante & Winter 2020; Amarante 2020b) using Eqs. (23)-(25).

Figures 2-4 show the relative errors of the gradients of the pseudo-potential between the simplified models and the polyhedron model for asteroid KBO Arrokoth, asteroid Kleopatra, and comet 103P∕Hartley, respectively. The plots display the relative error obtained when modeling the body as a point mass (red), with the DSM (blue), and with the GDSM (green).

As illustrated in Fig. 2, when the spacecraft approaches KBO Arrokoth, particularly at distances shorter than 20 km from its center, the point-mass approximation becomes inadequate because it produces significant relative errors in the computed gradient of the pseudo-potential. When the body is represented by either the DSM or the GDSM formulation, however, the resulting gradient of the pseudo-potential is strong consistent with the reference model. In regions far from the surface of the asteroid, the relative error obtained using the point-mass model, the DSM and the GDSM becomes comparable, indicating that the irregular shape of the body has little effect at large distances. In regions close to the surface of KBO Arrokoth, the relative error in the gradient of the pseudo-potential computed using the GDSM remains at about 7%. In contrast, the DSM presents errors of almost 15% under the same conditions. This significant difference highlights the superiority of the GDSM in accurately representing the dynamical environment near the asteroid. Consequently, for spacecraft dynamics analyses in the vicinity of KBO Arrokoth, the GDSM provides a more realistic and reliable approximation.

Figure 3 illustrates the relative error in the magnitude of the gradient of the pseudo-potential for three different gravitational representations of comet Harltey2 when compared to the polyhedron model. The point-mass approximation (red) differs most, particularly near the surface, where its simplified structure fails to capture the highly irregular mass distribution of the body. The dipole-segment model (DSM), shown in blue, reduces these errors, but still deviates noticeably between 0 and 2 km from the comet surface. In contrast, the errors of the generalized dipolesegment model (GDSM), shown in green, remain nearly zero at distances greater than 0.5 km from the surface. The GDSM curve remains tightly concentrated around the reference line and agrees excellently with the polyhedron model. This behavior highlights the improved capacity of the GDSM to reproduce the true gravitational environment of elongated bodies while retaining the low computational cost of simplified models.

For the case of asteroid Kleopatra (Fig. 4), a similar trend is observed. At large distances from the surface, all three formulations, the mass-point, the DSM, and the GDSM, provide comparable results, with relative errors that rapidly diminish with increasing distance. As the evaluation points approach the surface of the elongated body, however, the differences between the models become more pronounced. In this regime, the GDSM formulation is a clear improvement over the mass-point and the DSM representations because it maintains significantly lower errors and is markedly superior in capturing the local gravitational environment near the surface.

The results demonstrate that the DSM already performs well near the surfaces of irregular bodies overall. It is a clear improvement over the point-mass model, whose errors become large under strong geometric asymmetries. The GDSM consistently achieves a noticeably higher level of accuracy, however, particularly in the regions in which the gravitational field varies sharply due to complex shape features. The GDSM preserves low errors across a wider range of distances and represents surface-driven irregularities with significantly better fidelity.

Fig. 2 Relative error of the gradient of the pseudo-potential for KBO Arrokoth.

Fig. 3 Relative error of the gradient of the pseudo-potential for comet 103P∕Hartley.

Fig. 4 Relative error of the gradient of the pseudo-potential for asteroid Kleopatra.

6 Heteroclinic orbits

Heteroclinic orbits are trajectories in phase space that connect two different equilibrium points of a dynamical system. In the case of the Restricted Three-Body Problem (RTBP; Gómez et al. 1988), the dipole-segment model (Elipe et al. 2021a; Abad et al. 2024), and similarly for the GDSM, heteroclinic orbits connect the triangular equilibrium points E₂ and E₄ when these points are unstable.

The heteroclinic orbits represent a natural path allowed by the system dynamics in conservative problems to change from one configuration to the next without the need for continuous propulsion. They are important for planning low-fuel missions because they exploit only the natural dynamics of the system.

In the GDSM, the dynamical system is described by the equations of motion, Eqs. (10), (11), and (12), and the equilibrium points are determined according to Eq. (14) and Table 2. We therefore observe that the triangular points are unstable, as indicated in Table 5.

If p is an equilibrium point of a vector field X defined on a certain manifold M, knowing that the stable manifolds (W^s) are the set of states that eventually approach the point and the unstable manifolds (W^u) are the set of states that came from the point before, then we have $$W^s(p)=\{\mathbf{x}\in M|\underset{t\to +\mathrm{\infty }}{lim}{\varphi }_t(\mathbf{t}\mathbf{,}\mathbf{x})=p\}$$(27) $$W^u(p)=\{\mathbf{x}\in M|\underset{t\to -\mathrm{\infty }}{lim}{\varphi }_t(\mathbf{t}\mathbf{,}\mathbf{x})=p\},$$(28)

where φ_t is the flux of the dynamical system. Now, we have the mathematical conditions for the heteroclinic orbit, in the GDSM, for p = E₂ and p = E₄, that is, a heteroclinic orbit is characterized by the intersection of the two manifolds, and any intersection point corresponds to an initial condition that generates a trajectory.

According to Gómez et al. (1988), the intersections between the manifolds on the x-axis, which provide the initial condition of the orbit, can be obtained in the (x, ẋ) plane because y = 0 and ẏ(x, C), with C being the Jacobi constant at the triangular equilibrium point. The graphs define the manifolds at the first crossing with the x-axis.

Next, we present the (x, ẋ) plane for each body we studied for stable and unstable manifolds (see the vertical dashed red lines that limit the size of the body), following by the respective heteroclinic orbits. Table 6 presents the initial conditions (x, C) of each orbit for each body.

Fig. 5 shows for the stable manifold and Fig. 6 for the unstable manifold that they have two initial conditions of heteroclinic orbits for Arrokoth at x = 0 because the region between the vertical dashed red lines lies within the body. The conditions allocated there are physically impossible.

Four heteroclinic orbits were found: U_n and U_p are the unstable manifolds that leave E₂ and arrive at E₄, and S_n and S_p are their respective stable manifolds, which leave E₄ and arrive at E₂, as shown in Fig. 7.

In the sequence, we searched the results for asteroid Kleopatra. The stable and unstable manifolds are presented in Figs. 8 and 9.

Figure 10 shows the four heteroclinic orbits that were found. The result is consistent with Abad et al. (2024) for the dipolesegment model (DSM). The variation in the shape of the extreme poles of the body from circular (DSM) to ellipsoidal (GDSM) caused small variations in the initial conditions of the orbits 6 and their magnitudes.

According to Abad et al. (2024), the Kleopatra heteroclinic orbits combined in stable-unstable pairs form symmetric periodic orbits of the families (see Abad et al. (2024) for the shape of the families),

• S_n - U_p → Family l

• S_n - U_n → Family b

• S_p - U_p → Family a

• S_p — U_n → Family k.

Finally, we have the stable and unstable manifolds at the first crossing with the x-axis for 103P/Hartley (see Figs. 11 and 12) and their respective heteroclinic orbits (Fig. 13).

The heteroclinic orbits for the three bodies modeled by GDSM are similar, with differences in magnitudes and initial conditions. Orbits with smaller amplitudes are found for bodies whose distance from the origin to the extremity of the body is greater, and smaller amplitudes for the opposite (Arrokoth and 103P/Hartley). For Kleopatra, whose distance is more balanced and ranges from (-0.4778, 0.4762), the orbits tend to have a balanced amplitude as well.

Fig. 5 Stable manifold (Ws(E2)) at the first crossing for Arrokoth.

Fig. 6 Unstable manifold (Wu(E2)) at the first crossing for Arrokoth.

Fig. 7 Four heteroclinic orbits of Arrokoth, Un, Up, Sn, and Sp.

Fig. 8 Stable manifold (Ws (E2)) at the first crossing for Kleopatra.

Fig. 9 Unstable manifold (Wu (E2)) at the first crossing for Kleopatra.

Fig. 10 Four heteroclinic orbits of Kleopatra: Un, Up, Sn, and Sp.

Fig. 11 Stable manifold (Ws(E2)) at the first crossing for 103P∕Hartley.

Fig. 12 Unstable manifold (Wu(E2)) at the first crossing for 103P∕Hartley.

Fig. 13 Four heteroclinic orbits of 103P/Hartley: Un, Up, Sn, and Sp.

7 Conclusions

We developed the GDSM to represent the gravitational field of elongated and irregular bodies. The main improvement in the GDSM lies in the representation of the flattened shape of the poles, which differs from the DSM and allows us to approximate the mass distribution of the body better. This improvement leads to a greater reliability, particularly in regions close to the surface. We applied the model to the KBOs Arrokoth and Kleopatra and to comet 103P/Hartley, and the results showed that the GDSM outperforms the point-mass modeling and the dipolesegment modeling. The GDSM provides more accuracy while preserving the advantage of low computational cost. It becomes a powerful and efficient tool for studying the gravitational environment around small celestial bodies and takes a step toward a more realistic scenario while still allowing for a computationally unburdened code in comparison with the expensive polyhedral model. This can be useful for early stages of a mission design for a visit to this type of asteroid. Furthermore, the GDSM can also be adopted in particular for real-time orbital corrections. This needs to be evaluated onboard with the few computational resources.

From the GDSM, we calculated the initial conditions and integrated the heteroclinic orbits that connect the unstable triangular points for Arrokoth, Kleopatra, and comet 103P/Hartley. For each case, four heteroclinic orbits were found, two stable and two unstable orbits. The orbits with smaller amplitudes (Arrokoth and 103P/Hartley) are around bodies whose distance from the origin to the extremity of the body is greater, and smaller amplitudes for the opposite. For Kleopatra, whose distance from the extremities is more balanced, the orbits tend to have balanced amplitudes as well. The combination of each pair of heteroclinic orbits forms a symmetrical planar periodic orbit around the body.

Acknowledgements

The authors wish to express their appreciation for the support provided by grant #309089/2021-2, #443116/2023-7, and #201718/2025-1 from the National Council for Scientific and Technological Development (CNPq), grant #2022/11783-5 from São Paulo Research Foundation (FAPESP). A. Amarante and A. Ferreira thank the financial support of the São Paulo Research Foundation (FAPESP) [grants #2023/11781-5 & #2025/1543 8-9]. Additionally, Santos, L.B.T. would like to thank the support of the Polytechnic School of the University of Pernambuco (POLI-UPE) and the PostGrad Program in Systems Engineering (PPGES). The authors also acknowledge the financial support of the Coordination for the Improvement of Higher Education Personnel (CAPES) - Brazil (finance code 001). This research was also supported by computational resources supplied by the Center for Scientific Computing (NCC/GriduNESP) of the São Paulo State University (UNESP) and the Center for Mathematical Sciences Applied to Industry (CeMEAI), funded by FAPESP [grant #2013/07375-0]. F.M. thanks the fellowship granted by FAPESP (#2024/16260-6). E. Tresaco thanks Grant PID2024-156002NB-I00 funded by MICIU/AEI/10.13039/501100011033/FEDER,UE.

References

All Tables

All Figures

	Fig. 1 Geometric configuration of the rotating coordinate system for the generalized dipole-segment model, where m1 and m2 represent the masses at the poles, and m3 corresponds to the mass of the central segment.
In the text

	Fig. 2 Relative error of the gradient of the pseudo-potential for KBO Arrokoth.
In the text

	Fig. 3 Relative error of the gradient of the pseudo-potential for comet 103P∕Hartley.
In the text

	Fig. 4 Relative error of the gradient of the pseudo-potential for asteroid Kleopatra.
In the text

	Fig. 5 Stable manifold (Ws(E2)) at the first crossing for Arrokoth.
In the text

	Fig. 6 Unstable manifold (Wu(E2)) at the first crossing for Arrokoth.
In the text

	Fig. 7 Four heteroclinic orbits of Arrokoth, Un, Up, Sn, and Sp.
In the text

	Fig. 8 Stable manifold (Ws (E2)) at the first crossing for Kleopatra.
In the text

	Fig. 9 Unstable manifold (Wu (E2)) at the first crossing for Kleopatra.
In the text

	Fig. 10 Four heteroclinic orbits of Kleopatra: Un, Up, Sn, and Sp.
In the text

	Fig. 11 Stable manifold (Ws(E2)) at the first crossing for 103P∕Hartley.
In the text

	Fig. 12 Unstable manifold (Wu(E2)) at the first crossing for 103P∕Hartley.
In the text

	Fig. 13 Four heteroclinic orbits of 103P/Hartley: Un, Up, Sn, and Sp.
In the text