p_{2}\\\ p_{3}\\\ \vdots\\\ p_{N+1}\end{pmatrix}$

The matrix appearing is a Vandermonde matrix which is invertible, so each
function $a_{i}$ is actually a linear combination of the polynomials
$p_{1},\dots,p_{N}$. Therefore, (B.1) shows that the function is a polynomial.
∎

##  Appendix C Centralizers of ergodic homogeneous actions

######  Theorem C.1.

Let $G$ be a simply connected Lie group,
$\phi_{t}:\mathrm{G}/\Gamma\to\mathrm{G}/\Gamma$ be an ergodic homogeneous
flow generated by an $\mathbb{R}$-semisimple element. Let
$Z_{\operatorname{Lip}}$ denote the set of Lipschitz transformations commuting
with $\phi_{t}$. Then $Z_{\operatorname{Lip}}=Z_{\mathrm{Aff}}$, the group of
affine transformations commuting with $\phi_{t}$.

The principle tool in proving Theorem C.1 is a result of Zeghib, which we
summarize here. If $X$ is a $C^{\infty}$ manifold, a subset $N\subset X$ is
called rectifiable if it is the Lipschitz image of a bounded subset of
$\mathbb{R}^{n}$ for some $n$. If $N$ is a rectifiable subset, it has a well-
defined Hausdorff measure, which we denote by $\mu_{N}$.

######  Theorem C.2 (Théorèm A, [56]).

Let $G$ be a simply connected Lie group, $\phi_{t}:G/\Gamma\to G/\Gamma$ be an
ergodic homogeneous flow, generated by an $\mathbb{R}$-semisimple element. If
$N\subset X$ is a rectifiable, $\phi_{t}$ invariant set, then $\mu_{N}$ is
$\phi_{t}$-invariant, and the ergodic components of $\mu_{N}$ are the Haar
measures on closed $H$-orbits in $G/\Gamma$, where $H$ is some fixed closed
subgroup $H\subset G$.

######  Remark C.3.

In fact, Zeghib claims more, by working with bi-homogeneous flows, but we will
only use the homogeneous flow version here.

###### Proof of Theorem C.1.

Let $f\in Z_{\operatorname{Lip}}$, and consider the graph

$N=\left\\{(x,f(x)):x\in G/\Gamma\right\\}\subset(G\times
G)/(\Gamma\times\Gamma).$

By construction, since $f$ commutes with $\phi_{t}$, $N$ is
$\phi_{t}\times\phi_{t}$-invariant. Furthermore, since $f$ is Lipschitz, $N$
is rectifiable. Therefore, by Theorem C.2, $\mu_{N}$, the Hausdorff measure on
$N$, is $\phi_{t}$-invariant. $\phi_{t}$ is ergodic on $G/\Gamma$,
$\phi_{t}\times\phi_{t}$ must be ergodic on $N$ as well.

Therefore, there exists a unique subgroup $H\subset G\times G$ such that
$N=(g,h)\cdot H/(\Gamma\times\Gamma)$ for some $(g,h)\in G\times G$. Since
$(e\Gamma,f(e\Gamma))\in N$, we may without loss of generality assume $g=e$.

Finally, consider the Lipschitz transformation $\pi:N\to G/\Gamma$ defined by
$\pi(x,y)=x$. Since $N$ is the graph of a Lipschitz transformation, $\pi$ is a
Lipschitz homeomorphism. This further implies that $\pi$ is a diffeomorphism,
since $N$ is a coset of $H$. Hence, $\dim(H)=\dim(G)$, and for each
$X\in\operatorname{Lie}(G)$, there exists a unique
$\bar{F}(X)\in\operatorname{Lie}(G)$ such that
$(X,\bar{F}(X))\in\operatorname{Lie}(H)$. It follows immediately that
$\bar{F}$ is a Lie algebra homomorphism since $H$ is a subgroup. Let $F:G\to
G$ denote the lift of $\bar{F}$ to $G$.

It is immediate that $H=\left\\{(g,F(g)):g\in G\right\\}$. Therefore, from the
definition of $N$, we get that $f(g\Gamma)=hF(g)\Gamma$, and $f$ is affine, as
claimed. ∎

##  Appendix D Lifting $\mathrm{G}$-actions

######  Theorem D.1.

Let $r=\infty$ or $2$ and $G$ be a real semisimple group such that every
simple factor has real rank at least 2. Let $G\curvearrowright X$ be a $C^{r}$
action which is totally Anosov with respect to some split Cartan subgroup $A$,
and assume that

  * •

$\pi:Y\to X$ is a $C^{r}$ fiber bundle over $X$,

  * •

the action of $C_{G}(A)$ lifts to a $C^{r}$ action on $Y$,

  * •

for each Weyl chamber $\mathfrak{C}$ of $A\subset G$ there exists some
distinguished $a\in\mathfrak{C}$ such that the lifted action of $a$ on $Y$ is
partially hyperbolic with respect to $E_{y}=T(\pi^{-1}(\pi(y))\oplus
T(C_{G}(A)\cdot y)$, the sum of the tangent bundle to the fibers of $\pi$ and
orbits of $C_{G}(A)$.

Then there exists some continuous action of $\tilde{G}$, the universal cover
of $G$, on $Y$ which is a lift of the $G$-action in the sense the if
$p:\tilde{G}\to G$ is the canonical projection, $\pi(g\cdot
x)=p(g)\cdot\pi(x)$. Furthermore, if there exists a continuous metric on $Y$
for which $da|_{T\pi^{-1}(\pi(y))}$ is isometric for every $a\in A$ and $y\in
Y$, the $\tilde{G}$-action is $C^{r}$.

The proof of this theorem follows the scheme introduced by the first author
and Katok in [14], extended by Zhenqi Wang in [52, 53] and made fully general
by the third author in [51]. We summarize some important definitions before
proceeding with the proof.

Let $\Delta_{G}$ denote the set of roots $\chi$ of $G$ such that $\chi/2$ is
not a root, and for each $\chi\in\Delta_{G}$, let $U_{\chi}$ denote
corresponding coarse Lyapunov subgroup. As described in Section 6.2, we
consider the free product $\mathcal{P}$ of the groups $U_{\chi}$, where $\chi$
ranges over $\Delta_{G}$, and note that since the groups $U_{\chi}$ generate
$G$, there exists a projection $\pi:\mathcal{P}\to\tilde{G}$ such that the
kernel $\mathcal{C}$ is exactly the expressions in $\mathcal{P}$ which yield
contractible cycles on $G$.

Since $G$ acts on $X$, any relations among the $U_{\chi}$ on $G$ hold on $X$
as well. Importantly, the commutator relations hold: if
$\chi_{1},\chi_{2}\in\Delta_{G}$ are linearly independent, $u\in U_{\chi_{1}}$
and $v\in U_{\chi_{2}}$, then
$[u,v]\in\displaystyle\prod_{\begin{subarray}{c}\chi=s\chi_{1}+t\chi_{2}\\\
s,t\in\mathbb{Z}_{+}\end{subarray}}U_{\chi}$. Hence, the commutator relations
hold on $X$. The following lemma becomes of particular interest:

######  Lemma D.2 (Lemma 4.7, [51]).

Let $\mathcal{P}$ denote the free product of the coarse Lyapunov subgroups of
$G$, and $\mathcal{C}_{S}\subset\mathcal{P}$ denote the normal closure (in
$\mathcal{P}$) of the group generated by the relations
$v^{-1}*u^{-1}*v*u*\rho^{-1}\in\mathcal{P}$, where $\rho$ is any presentation
of $[u,v]$ in the group
$\displaystyle\prod_{\begin{subarray}{c}\chi=s\chi_{1}+t\chi_{2}\\\
s,t\in\mathbb{Z}_{+}\end{subarray}}U_{\chi}$. Then $\mathcal{C}_{S}$ is a co-
abelian subgroup of $\mathcal{C}$, and every continuous action of
$\mathcal{C}/\mathcal{C}_{S}$ on a finite dimensional space is trivial.

###### Proof of Theorem D.1.

Let $a\in A$ be one of the distinguished elements of the Weyl chambers which
lift to parially hyperbolic maps. Consider $W^{s}_{a,X}$ and $W^{s}_{a,Y}$,
the stable and unstable foliations of $a$ as it acts on $X$ and $Y$,
respectively. Since the fibers of $\pi$ are contained in the center foliation,
it follows that $\pi|_{W^{s}_{a,Y}(y)}$ is a $C^{r}$ diffeomorphism. For every
$\chi$, we may therefore build an action $U_{\chi}\curvearrowright Y$ by
choosing some $a\in\mathbb{R}^{k}$ such that $\chi(a)<0$ and letting $u\cdot
y$ denote the unique element of $W^{s}_{a,Y}(y)$ which projects to $u\cdot
x\in W^{s}_{a,X}(\pi(y))$. Note that the action does not depend on the choice
of $a$ or Weyl chamber which $a$ belongs to. By the universal property of
(topological) free products, we may construct an action of $\mathcal{P}$, the
free product of the groups $U_{\chi}$, on $Y$ as well.

Fix any pair of linearly independent roots $\chi_{1},\chi_{2}\in\Delta_{G}$.
Then, by linear independence, there exists $a\in\mathbb{R}^{k}$ such that
$\chi_{1}(a),\chi_{2}(a)<0$. If $u\in U_{\chi_{1}}$ and $v\in U_{\chi_{2}}$,
then by Lemma D.2, there exists a $w$ written as a product of elements from
the groups $U_{\chi}$, $\chi=s\chi_{1}+t\chi_{2}$, $s,t\in\mathbb{Z}_{+}$ such
that $[u,v]w^{-1}=e$, as elements in $G$. Since the relations hold on $G$,
they act trivially on $X$. Furthermore, since the orbits of $U_{\chi_{1}}$,
$U_{\chi_{2}}$ and every group appearing in the presentation of $w$ are
contained in $W^{s}_{a,X}$, the relation holds when the groups act on $Y$ as
well.

Therefore, the restriction of the $\mathcal{P}$-action on $Y$ to
$\mathcal{C}_{S}$ is trivial. Hence there exists a well-defined continuous
action of $\mathcal{C}/\mathcal{C}_{S}$ on $Y$. This action must be trivial by
Lemma D.2. Finally, since the $\mathcal{C}$ action is trivial, the action of
$\mathcal{P}$ induces an action of $\tilde{G}=\mathcal{P}/\mathcal{C}$ on $Y$,
which by construction must cover the action on $X$.

To prove the regularity of the action, notice that the lift of the actions of
the groups $U_{\chi}$ to $Y$ is exactly determined by stable holonomies.
Since, restricted to $T\pi^{-1}(\pi(y))$, the action is isometric, the stable
holonomies will have arbitrarily good pinching properties. By now-standard
arguments in partial hyperbolicity (see, eg, [43, Theorem 6.1]), it follows
that the holonomies, and hence action of each of the groups $U_{\chi}$ on $Y$
is $C^{r}$. Finally, since the $U_{\chi}$ generate $\tilde{G}$, every element
of $\tilde{G}$ is a $C^{r}$ transformation of $Y$. By Theorem 5.5, the action
$\tilde{G}\curvearrowright Y$ is a $C^{r}$ group action. ∎

######  Remark D.3.

The smoothness of the lift can be obtained from a weaker assumption, namely
center bunching. To obtain a $C^{\infty}$ lift, one requires a more
restrictive form of center bunching than that laid out in Section 5.1, which
allows for no exponents on the fiber whatsoever (i.e., something imitating
unipotent behavior). Since our application is to a compact group extension,
and the action of a semisimple Lie group will always be conformal, we use this
simpler version of the statement.

## References

  * [1] A. Avila, J. Santamaria and M. Viana, 2013. Holonomy invariance: rough regularity and applications to Lyapunov exponents. _Astérisque_ , 358, pp.13-74.
  * [2] A. Avila, M. Viana, (2010). Extremal Lyapunov exponents: an invariance principle and applications. _Inventiones mathematicae_ , 181(1), 115-178.
  * [3] Y. Benoist, P. Foulon and F. Labourie, (1992). Flots d’Anosov à distributions stable et instable différentiables. Journal of the American Mathematical Society, 5(1), 33-74.
  * [4] A. Brown, F. Rodriguez Hertz, and Z. Wang. Smooth ergodic theory of $Z^{d}$-actions. arXiv:1610.09997.
  * [5] A. Brown, D. Fisher and S. Hurtado. Zimmer’s conjecture for actions of ${\mathrm{S}L}(m,\mathbb{Z})$, _Inventiones mathematicae_ , 221 (3), 1001–1060 (2020).
  * [6] A. Brown, D. Fisher and S. Hurtado. Zimmer’s conjecture: Subexponential growth, measure rigidity, and strong property (T), to appear in _Annals of Mathematics_.
  * [7] A. Brown, D. Fisher and S. Hurtado. Zimmer’s conjecture for non-uniform lattices and escape of mass. arXiv:2105.14541
  * [8] A. Brown, F. Rodriguez Hertz, and Z. Wang. Global smooth and topological rigidity of hyperbolic lattice actions. _Ann. of Math._(2), 186(3):913–972, 2017.
  * [9] K. Burns and A. Wilkinson, 2010. On the ergodicity of partially hyperbolic systems. _Annals of Mathematics._ , pp. 451-489.
  * [10] C. Butler. Rigidity of equality of Lyapunov exponents for geodesic flows. Journal of Differential Geometry 109, no. 1 (2018): 39-79.
  * [11] P. Chernoff and J. Marsden (1970). On continuity and smoothness of group actions. Bulletin of the American Mathematical Society, 76(5), 1044-1049.
  * [12] D. Damjanović, A. Wilkinson and D. Xu, D. (2021). Pathology and asymmetry: centralizer rigidity for partially hyperbolic diffeomorphisms. Duke Mathematical Journal, 170(17), 3815-3890.
  * [13] D. Damjanović and D. Xu. On classification of higher rank Anosov actions on compact manifold, _Isr. J. Math._ 238, 745–806 (2020).
  * [14] D. Damjanović and A. Katok, Local Rigidity of Partially Hyperbolic Actions. II: The Geometric Method and Restrictions of Weyl Chamber Flows on $SL(n,\mathbb{R})/\Gamma$. International Mathematics Research Notices. 2011; 2011(19):4405-30.
  * [15] V. Deodhar, On Central Extensions of Rational Points of Algebraic Groups. American Journal of Mathematics, Vol. 100, No. 2 (Apr., 1978), pp. 303-386.
  * [16] R. Feres and F. Labourie. Topological superrigidity and Anosov actions of lattices, _Ann. Sci. École Norm. Sup._ 31 (4), 1998, 599–629.
  * [17] D. Fisher. Recent developments in the Zimmer program. _Notices Amer. Math. Soc._ 67 (2020), no. 4, 492–499.
  * [18] D. Fisher. Groups acting on manifolds: around the Zimmer program, Geometry, rigidity, and group actions, Chicago Lectures in Math., Univ. Chicago Press, Chicago, IL, 2011, pp. 72–157, DOI 10.7208/chicago/9780226237909.001.0001.
  * [19] E. Goetze and R. J. Spatzier. On Livšic’s theorem, superrigidity, and Anosov actions of semisimple Lie groups. _Duke Mathematical Journal_ , 88 (1), 1–27, 1997.
  * [20] Gleason, Andrew M.; Palais, Richard S. On a class of transformation groups. _A_ mer. J. Math. 79 (1957), 631–648.
  * [21] E. Goetze and R. J. Spatzier. Smooth classification of Cartan actions of higher rank semisimple Lie groups and their lattices. _Ann. of Math._ (2), 150(3):743–773, 1999.
  * [22] M. I. Graev. On free products of topological groups. _Izvestiya Akad. Nauk SSSR. Ser. Mat._ , 14:343–354, (1950).
  * [23] Humphreys, James E. Introduction to Lie algebras and representation theory. Second printing, revised. _Graduate Texts in Mathematics_ , 9. Springer-Verlag, New York-Berlin, 1978.
  * [24] S. Hurder, 1992. Rigidity for Anosov actions of higher rank lattices. Annals of Mathematics, 135(2), pp.361-410.
  * [25] J. L. Journé, (1988). A regularity lemma for functions of several variables. Revista Matematica Iberoamericana, 4(2), 187-193.
  * [26] B. Kalinin, Non-stationary normal forms for contracting extensions, to appear in Dynamics and related topics, 2020 vision and agenda.
  * [27] B. Kalinin and V. Sadovskaya, 2016. Normal forms on contracting foliations: smoothness and homogeneous structure. Geometriae Dedicata, 183(1), pp.181-194.
  * [28] B. Kalinin and V. Sadovskaya. Cocycles with one exponent over partially hyperbolic systems. _Geom Dedicata_ 167, 167–188 (2013).
  * [29] B. Kalinin and V. Sadovskaya, 2007. On classification of resonance-free Anosov $\mathbb{Z}^{k}$ actions. Michigan Mathematical Journal, 55(3), pp.651-670.
  * [30] B. Kalinin and R. J. Spatzier. On the classification of Cartan actions. _Geometric And Functional Analysis_ , (2007), 17(2), 468–490.
  * [31] A. Katok and R. J. Spatzier. Differential rigidity of Anosov actions of higher rank abelian groups and algebraic lattice actions. Tr. Mat. Inst. Steklova, 216(Din. Sist. i Smezhnye Vopr.):292–319, 1997.
  * [32] A. Katok A, J. Lewis, R. Zimmer, Cocycle superrigidity and rigidity for lattice actions on tori. Topology. 1996 Jan 1;35(1):27-38.
  * [33] A. Katok and J. Lewis. Local rigidity for certain groups of toral automorphisms. Israel J. Math., 75(2-3):203–241, 1991.
  * [34] A. Katok and J. Lewis. Global rigidity results for lattice actions on tori and new examples of volume-preserving actions. Israel J. Math., 93:253–280, 1996.
  * [35] W. Kryszewski and S. Plaskacz. Topological methods for the local controllability of nonlinear systems. SIAM J. Control Optim. 32 (1994), no. 1, 213–223.
  * [36] F. Ledrappier. Positivity of the exponent for stationary sequences of matrices. In Lyapunov exponents (Bremen, 1984), volume 1186 of Lect. Notes Math., 56–73. Springer, 1986.
  * [37] F. Ledrappier and L. S. Young. The metric entropy of diffeomorphisms: part ii: relations between entropy, exponents and dimension. _Annals of Mathematics_ , Vol. 12,(3), p. 540-574, 1985).
  * [38] G. A. Margulis, 2000. Problems and conjectures in rigidity theory. Mathematics: frontiers and perspectives. Amer. Math. Soc., Providence, RI.
  * [39] G. A. Margulis, 1991. Discrete subgroups of semisimple Lie groups (Vol. 17). Springer Science $\&$ Business Media.
  * [40] G. A. Margulis. Discrete groups of motions of manifolds of nonpositive curvature. Proceedings of the International Congress of Mathematicians (Vancouver, B.C., 1974), Vol. 2, 21–34.
  * [41] G. A. Margulis. Arithmeticity of the irreducible lattices in the semisimple groups of rank greater than $1$. _Inventiones Mathematicae_ , 93–120, 76, (1984).
  * [42] D. Montgomery and L. Zippin. Topological transformation groups. Courier Dover Publications, 2018.
  * [43] C. Pugh, M. Shub and A. Wilkinson, Hölder foliations, Duke Math. J., 86 (1997), no. 3, 517–546
  * [44] E. T. Ordman. Free products of topological groups which are $k_{\omega}$-spaces. _Trans. Amer. Math. Soc._ , 191:61– 73, 1974.
  * [45] D. Repovs, E. Scepin. and A. Skopenkov, $C^{1}$-homogeneous compacta in $\mathbb{R}^{n}$ are $C^{1}$-submanifolds of Rn. Proc. Amer. Math. Soc., 124:1219–1226, 1996.
  * [46] B. Schmidt. Weakly hyperbolic actions of Kazhdan groups on tori. _Geom. Funct. Anal._ 16 (5):1139–1156, 2006.
  * [47] S. Smale. Differentiable dynamical systems. _Bull. Amer. Math. Soc._ , 73:747–817, 1967.
  * [48] R. J. Spatzier and K. Vinhage. Cartan Actions of Higher Rank Abelian Groups and their Classification. _JAMS_ , DOI: 10.1090/jams/1033.
  * [49] T. Tao. Hilbert’s fifth problem and related topics, volume 153 of _Graduate Studies in Mathematics. American Mathematical Society_ , Providence, RI, 2014.
  * [50] M. Taylor, Existence and regularity of isometries. _Trans. Amer. Math. Soc._ , 358 (6), 2415–2423, 2006.
  * [51] K. Vinhage. On the rigidity of Weyl chamber flows and Schur multipliers as topological groups. _J. Mod. Dyn._ , 9:25–49, 2015.
  * [52] Z. J. Wang, 2010. Local rigidity of partially hyperbolic actions. Journal of Modern Dynamics, 4(2), p.271.
  * [53] Z. J. Wang, (2010). New cases of differentiable rigidity for partially hyperbolic actions: Symplectic groups and resonance directions. Journal of Modern Dynamics, 4(4), 585.
  * [54] A. Wilkinson. The cohomological equation for partially hyperbolic diffeomorphisms. _Astérisque_ , 358 (2013).
  * [55] Edward N. Wilson. Isometry groups on homogeneous nilmanifolds. _Geometriae Dedicata_ , volume 12, pages 337–346 (1982).
  * [56] A. Zeghib. Ensembles invariants des flots géodésiques des variétés localement symétriques. _Ergodic Theory Dynam. Systems_ 15 (1995), no. 2, 379–412.
  * [57] R. J. Zimmer. Ergodic Theory of Semisimple Groups, _Birkhäuser_ , Monographs in Mathematics (1984).
  * [58] R.J. Zimmer. Actions of semisimple groups and discrete subgroups. In _Proceedings of the International Congress of Mathematicians_ , Vol. 1, 2 (Berkeley, Calif., 1986), pages 1247–1258. Amer. Math. Soc., Providence, RI, 1987.
  * [59] R. J. Zimmer. Arithmetic groups acting on compact manifolds. Bull. Amer. Math. Soc. (N.S.), 8(1):90– 92, 1983.
  * [60] R. J. Zimmer. Volume preserving actions of lattices in semisimple groups on compact manifolds. _Inst. Hautes Études Sci. Publ. Math._ , (59):5–33, 1984.
  * [61] R. J. Zimmer. Unpublished notes.


# A Particle-Linkage Model for Elongated Asteroids with Three-Dimensional Mass
Distribution

L. B. T. Santos1, L. O. Marchi1, S. Aljbaae1, P. A. Sousa-Silva3, D. M.
Sanchez1 E-mail<EMAIL_ADDRESS>A. F. B. A. Prado1
1Division of Space Mechanics and Control. National Institute for Space
Research, INPE. São Paulo, Brazil.
2São Paulo State University, UNESP. São João da Boa Vista, Brazil

(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

The goal of the present paper is to develop a simplified model to describe the
gravitational fields of elongated asteroids. The proposed model consists of
representing an elongated asteroid using a triple-particle-linkage system
distributed in the three-dimensional space and it is an extension of previous
planar models. A non-linear optimisation method is used to determine the
parameters of our model, minimising the errors of all the external equilibrium
points with respect to the solutions calculated with a more realistic
approach, the Mascon model, which are assumed to give the real values of the
system. The model considered in this article is then applied to three real
irregular asteroids: 1620 Geographos, 433 Eros, and 243 Ida. The results show
that the current triple-particle-linkage three-dimensional model gives better
accuracy when compared to the axisymmetric triple-particle-linkage model
available in the literature, and provides an advantage in terms of accuracy
over the mass point model, while keeping computational time low. This model is
also used to carry out simulations to characterise regions with solutions that
remain bounded or that escape from around each asteroid under analysis. We
investigated initial inclinations of $0^{\circ}$ (direct orbits) and
$180^{\circ}$ (retrograde orbits). We considered the gravitational field of
the asteroid, the gravitational attraction of the Sun, and the SRP. Our
results are then compared to the results obtained using the Mascon
gravitational model, based on the polyhedral shape source. We found good
agreement between the two models.

###### keywords:

gravitation fields– celestial mechanics – minor planets – asteroids:
individual: 1620 Geographos – asteroids: individual: 243 Ida – asteroids:
individual: 433 Eros.

††pubyear: 2020††pagerange: A Particle-Linkage Model for Elongated Asteroids
with Three-Dimensional Mass Distribution–A Particle-Linkage Model for
Elongated Asteroids with Three-Dimensional Mass Distribution

##  1 Introduction

For the past few decades, space agencies have focused their attention on
visiting the small objects that make up the Solar System, sending spacecraft
to asteroids and comets and having plans to send more in the future.

Most asteroids and comets have an irregular shape. When planning a space
mission to visit these objects, a question is raised: “how to effectively
represent the gravitational field of these small irregular bodies?”(Elipe &
Lara, 2003; Bartczak & Breiter, 2003). Several mathematical models have been
proposed to overcome this difficulty.

Spherical harmonic expansion is often used to model planets, as these
celestial bodies have a sphere-like geometrical shape (Elipe & Lara, 2003).
However, when the object has an irregular shape, harmonic expansion is no
longer convenient and, often, convergence cannot be guaranteed (Riaguas et
al., 1999; Jiang & Baoyin, 2018; Lan et al., 2017). A model widely used to
represent the gravitational field of irregularly shaped objects is the
polyhedral model (Werner, 1994; Werner & Scheeres, 1996; Scheeres et al.,
2000) which is considered an accurate model and does not present convergence
problems (Scheeres et al., 2000). Due to its high accuracy, the polyhedral
model has been widely used in projects of real missions to analyze the
dynamics of the spacecraft near asteroids. As a drawback, the polyhedral model
generally contains thousands of parameters (vertices and facets) that are
necessary to guarantee its high accuracy. Therefore, this model requires a
high computational effort and time to perform the simulations. Besides, an
investigation regarding the general properties of the dynamics of the
spacecraft as a function of the parameters of the model becomes very
difficult. It happens because, in the polyhedral model, the parameters used
are mixed, producing a mixed influence on the effective gravitational field of
irregular bodies.

Simplified models can often be used to approximate the gravitational field of
irregular bodies, requiring less computational effort and producing
considerable results in a short amount of time. This topic has been addressed
in Chanut et al. (2015a) and Aljbaae et al. (2017) where the authors applied
the Mascon model using the polyhedral data to place the mass particles. As
shown by Aljbaae et al. (2017), the maximum difference between the classical
polyhedral approach (Werner, 1994; Werner & Scheeres, 1996) and the Mascon 8,
considering a uniform density, occurs at the location of the equilibrium point
represents 0.11 percent of the distance from the centre of the body, which, in
our opinion, may be consider satisfactory. That motivated us to use the Mascon
approach to validate our model later on in this work.

Besides that, a simplified model can often determine a simple analytical
expression for calculating the gravitational field of irregular bodies (Yang
et al., 2018). Another advantage of using a simplified model is that we can
easily analyse the consequences of a given parameter on the dynamics around
asteroids, such as, for example, the stability of the equilibrium points (Zeng
et al., 2015; Barbosa Torres dos Santos et al., 2017), the distribution of
stable periodic orbits (Lan et al., 2017), as well as the permissible hovering
regions Yang et al. (2015); Zeng et al. (2016a). Moreover, simplified models
can be used to assist orbit design (Wang et al., 2017) and feedback control
(Yang et al., 2017). Particular trajectories can be studied under more
accurate models when considering the final steps of a real mission, but a
general knowledge of the gravity field based on a small number of parameters
is very useful in the first stages of mission design.

Thus, several simplified models have been proposed to represent irregular
asteroids, such as the ones in Riaguas et al. (1999); Riaguas et al. (2001);
Bartczak & Breiter (2003); Elipe & Lara (2003), who investigated the dynamics
of a particle under the gravitational field of an asteroid modelled as a
straight segment. Other work that also investigated the motion of a particle
around small celestial bodies considering a simplified model used 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), a double and finite straight segments (Jain & Sinha, 2014), a rotating
mass dipole (Chermnykh, 1987; Kirpichnikov & Kokoriev, 1988; Kokoriev &
Kirpichnikov, 1988). In 2015, for the first time, Zeng et al. (2015) presented
a method for adjusting the rotating mass dipole model to a real asteroid. An
improvement of this model was carried out considering the oblateness of the
primary bodies (Zeng et al., 2016b, c, 2017), and the dipole segment model
(Zeng et al., 2018). Based on the elongated asteroid binary systems, Barbosa
Torres dos Santos et al. (2017) used the restricted synchronous three-body
problem model. These work, among others, showed that it is useful to apply
simplified models to identify the main parameters that dominate the motion of
a particle around certain asteroid systems.

Inspired by the work of Zeng et al. (2015), Lan et al. (2017) proposed the
rotating mass tripole model with symmetrical rotation. In the work developed
by Lan et al. (2017), it is argued that small convex bodies can be modeled
using the mass tripole model. Lan et al. (2017) showed that, from five
parameters (which are determined with the help of the polyhedral model), it is
possible to define the geometric shape and to obtain the physical
characteristics of a real asteroid and to obtain its gravitational field. In
order to generalize the work of Lan et al. (2017), Barbosa Torres dos Santos
et al. (2020) investigate the qualitative dynamics in the vicinity of an
asteroid with an convex shape using a rotating tripole model from a semi-
analytical study.

The axisymmetric triple-particle-linkage model, proposed by Lan et al. (2017),
considers that the points of mass are located in the $xy$ plane of the
asteroid (two-dimensional configuration). But we know that asteroids have a
spatial mass distribution ($xyz$ axes). Given this, the purpose of this work
is to improve the two-dimensional axisymmetric triple-particle-linkage model
(2D tripole) maintaining a simplified format, but considering a three-
dimensional mathematical modeling, taking another step towards a more
realistic scenario.

We rely on non-linear optimization to find the parameters of our proposed
model, as we will see later in this paper. The advantages of this model are
described next.

Because the proposed model considers the distribution of mass in three
dimensions, this model reproduces the dynamics of the spacecraft more
realistically when compared to the two-dimensional axisymmetric rotating mass
tripole (Lan et al., 2017), and, consequently, the rotating mass dipole model
(Zeng et al., 2015).

Although this is a simplified model, it is beneficial to carry out a
qualitative investigation to analyse the effects of the model’s parameters on
the orbital dynamics. Qualitative analyses can be used to investigate the
effects of some parameters, such as relating the characteristics of the motion
around elongated asteroids and the dynamical properties of the body (Yang et
al., 2015), or to analyse the distribution of stable periodic orbits close to
the equatorial plane (Lan et al., 2017), etc. Furthermore, qualitative
analysis has the potential to conduct preliminary real mission design (Yang et
al., 2018), leaving detailed studies for a later stage, based on more accurate
models that require more computer time, but will study a small number of
trajectories.

This paper is organized as follows. Section 2 provides the normalized equation
of motion of a particle around elongated asteroids that can be modeled as a
two/three-dimensional axisymmetric rotating mass tripole. A general equation
is provided, from where we can derive the equation of motion for the
situations of the planar or spatial mass of the asteroid, in the inertial and
in the rotating frame. Subsequently, the methodology adopted to determine the
parameters of the simplified model using a nonlinear optimization method is
presented in Section 3. In Section 4, the proposed methodology and model are
applied to three irregular asteroids: 1620 Geographos, 433 Eros and 243 Ida. A
comparison of the three-dimensional axisymmetric rotating mass tripole and
two-dimensional axisymmetric rotating mass tripole with respect to the Mascon
gravity framework using a shaped polyhedral source is performed in Section 5,
whose objective is to show the advantage of this new model. In Section 6, we
perform numerical simulations where we construct grids of initial conditions,
relating semi-major axis and eccentricity, to characterize bounded and
unbounded regions around the asteroids 1620 Geographos, 433 Eros and 243 Ida.
Finally, Section 7 provides our conclusions.

##  2 Equation of Motion

In this section, we present the mathematical model of the restricted four-body
problem, which is a basis for the three-dimensional rotating mass tripole
model.

Figure 1 shows three points of mass ($M_{1}$, $M_{2}$ and $M_{3}$) arranged
inside a body (asteroid or comet) that has an irregular shape. The equations
of motion developed in the rotating frame consider the asteroid-spacecraft
system, which means that perturbations from other celestial bodies and other
external forces are not considered. The rods that connect $M_{1}$ to $M_{3}$
and $M_{2}$ to $M_{3}$, shown in Fig. 1, are assumed to have negligible mass
and the same length $L$. This length $L$ (in canonical unit) can vary
depending on the asteroid (or comet) under study. The distance between $M_{1}$
and $M_{2}$ is denoted by $d$ and it defines one canonical unit (c.u.).

Figure 1: The rotating mass tripole in $xy$ plane.

Note, from Fig. 1, that it is possible to see the arrangement of the points of
mass in the $xy$ plane. Since this study uses a three-dimensional model, we
have to analyse the arrangement of these points of mass in space. Figure 2
provides a representative image of the position of $M_{1}$, $M_{2}$, and
$M_{3}$ in the $xyz$ space. To determine the locations of the primary bodies,
we use two angles, called azimuth ($\Phi$) and elevation angle ($\Psi$). Note
that the geometric configuration of the asteroid (or comet) depends on $\Phi$
and $\Psi$. Figure 2 shows the same geometric configuration of the primary
bodies. The segment green represents the position of the bodies $M_{1}$,
$M_{2}$, and $M_{3}$ with respect to the centre of mass of the asteroid. Note
that when $M_{1}$ and $M_{2}$ have a positive location on the $y$ axis,
$M_{3}$ has a negative position on the same axis. Likewise, when bodies
$M_{1}$ and $M_{2}$ have a positive position on the $z$ axis, $M_{3}$ has a
negative position on this same axis, such that the centre of mass of the
system remains at the origin of the Cartesian coordinate system (red vectors).
From Fig. 2, we can see that the position of $M_{1}$ depends on the angles
$\Phi$ and $\Psi$. $M_{2}$ is spatially positioned symmetrically to $M_{1}$.
Although it does not appear in Figure 2, the position of all primary bodies
($M_{1}$, $M_{2}$, and $M_{3}$) depends on the angles $\Phi$ and $\Psi$, as we
will see in the equations of motion. We chose not to place all the spatial
geometry of the primary bodies with respect to $\Phi$ and $\Psi$ to avoid
visual pollution in the figure.

When $\Psi=90^{\circ}$ and $\Phi\neq 0^{\circ}$, we have the particular case
of the rotating mass tripole model (two-dimensional). On the other hand, if
$\Psi=90^{\circ}$ and $\Phi=0^{\circ}$, which means that $M_{1}$, $M_{2}$, and
$M_{3}$ are aligned on the $x$ axis, there is symmetry with the equatorial
plane. We consider that the asteroid’s centre of mass is the origin of the
reference system ($xyz$).

Figure 2: Perspective of the positions of bodies $M_{1}$, $M_{2}$ and $M_{3}$
in the three-dimensional space.

In this work, we consider that both rods have length $L$, make the same angle
$\Phi$, in magnitude, with the $x$-axis and the same angle $\Psi$ with the
$z$-axis.

Next, we suppose that a body with negligible mass (spacecraft) is located at
P($x$, $y$, $z$). It is convenient to write the equations of motion on a unit
scale. For this, let’s assume that the distance between $M_{1}$ and $M_{2}$ is
the unit of length. The time unit is selected such that, the rotational period
of the tripole 3D is given by 2 _$\uppi$_. It is also assumed that the sum of
the masses ($M=m_{1}+m_{2}+m_{3}$) of the primaries is the unit of mass, where
$m_{1}$, $m_{2}$, and $m_{3}$ are the masses of $M_{1}$, $M_{2}$, and $M_{3}$,
in kilograms, respectively. The coordinates of the bodies $M_{1}$, $M_{2}$ and
$M_{3}$, in canonical units, are, respectively, given

$x_{1}=-L\cos{\Phi}\sin{\Psi},$ $y_{1}=(1-2\mu^{*})L\sin{\Phi},$ (1)
$z_{1}=L\cos{\Phi}\cos{\Psi},$

$x_{2}=L\cos{\Phi}\sin{\Psi},$ $y_{2}=(1-2\mu^{*})L\sin{\Phi},$ (2)
$z_{2}=L\cos{\Phi}\cos{\Psi},$

$x_{3}=0,$ $y_{3}=-2\mu^{*}L\sin{\Phi},$ (3)
$z_{3}=-2\mu^{*}L\cos{\Phi}\cos{\Psi}/(1-2\mu^{*}),$

where $\mu^{*}$ is the mass ratio defined by

$\mu^{*}=\frac{m_{2}}{m_{1}+m_{2}+m_{3}}.$ (4)

Note that the distance $d$, in canonical unit, between $M_{1}$ and $M_{2}$
depends on $L$, $\Phi$, and $\Psi$, which can be expressed by

$x_{2}-x_{1}=d=2L\cos{\Phi}\cos{\Psi}=1.$ (5)

In this work, we will use the notation $d$ to represent the distance between
$M_{1}$ and $M_{2}$ in canonical unit and $d^{*}$ to represent the distance
between $M_{1}$ and $M_{2}$ in meters.

Using the canonical units mentioned before, the Lagrangian function of the
system is given by:

$\mathcal{L}_{rotating}(x,y,z,\dot{x},\dot{y},\dot{z})=(\dot{x}^{2}+\dot{y}^{2}+\dot{z}^{2})/2+(x\dot{y}-y\dot{x})+(x^{2}+y^{2})/2\\\
+k(\mu^{*}/r_{1}+\mu^{*}/r_{2}+(1-2\mu^{*})/r_{3}),$ (6)

where

$r_{1}=\sqrt{(x-x_{1})^{2}+(y-y_{1})^{2}+(z-z_{1})^{2}},$ (7)
$r_{2}=\sqrt{(x-x_{2})^{2}+(y-y_{2})^{2}+(z-z_{2})^{2}},$ (8)
$r_{3}=\sqrt{(x-x_{3})^{2}+(y-y_{3})^{2}+(z-z_{3})^{2}}.$ (9)

$k$ is called the force ratio and is given by Eq. 10 (Zeng et al., 2015)

$k=\frac{GM}{\omega^{2}d^{*^{3}}}.$ (10)

Note that $k$ is a dimensionless quantity and depends on the angular velocity
vector of the asteroid, denoted by $\omega=[0,0,\omega]^{T}$, in radians per
second, $M$ is the mass of the body under study in kilogrammes, $d^{*}$ is the
distance between $M_{1}$ and $M_{2}$, in metres and $G$ is the universal
gravitation constant, given in the international system of units (Lan et al.,
2017; Zeng et al., 2015).

The Hamiltonian formulation can be obtained from the Lagrangian formulation.
The conjugate momenta ($p$) are given by

$p_{x}=\frac{\partial\mathcal{L}_{rotating}}{\partial\dot{x}}=\dot{x}-y,$ (11)
$p_{y}=\frac{\partial\mathcal{L}_{rotating}}{\partial\dot{y}}=\dot{y}+x,$ (12)
$p_{z}=\frac{\partial\mathcal{L}_{rotating}}{\partial\dot{z}}=\dot{z}.$ (13)

Using Legendre transformation
($q,\dot{q})\overset{\Lambda}{\rightarrow}(\frac{\partial\mathcal{L}_{rotating}}{\partial\dot{q}}$)
where $q,\dot{q}$ are coordinates and velocities of the particle,
respectively, we get

$\mathcal{H}=\dot{x}p_{x}+\dot{y}p_{y}+\dot{z}pz-\mathcal{L}_{rotating}.$ (14)

After some simplifications, we obtain the corresponding Hamiltonian

$\mathcal{H}(x,y,z,p_{x},p_{y},p_{z})=\frac{p_{x}^{2}+p_{y}^{2}+p_{z}^{2}}{2}-xp_{y}+yp_{x}-\\\
k\left(\frac{\mu^{*}}{r_{1}}+\frac{\mu^{*}}{r_{2}}+\frac{1-2\mu^{*}}{r_{3}}\right).$
(15)

Note that the Hamiltonian $\mathcal{H}$ is a function of the coordinates and
the conjugate momenta $p$ of the spacecraft. Due to the fact that the
Hamiltonian is independent of time, it is possible to say that its value is
conserved (Benacchio, 2007; Broucke, 1968; Worthington, 2012). $p_{x}$,
$p_{y}$ and $p_{z}$ are associated with the spacecraft’s conjugate momenta
with respect to the $x$-axis, $y$-axis and $z$-axis, respectively.

We define the effective potential by

$\Omega=\frac{x^{2}+y^{2}}{2}+k\left(\frac{\mu^{*}}{r_{1}}+\frac{\mu^{*}}{r_{2}}+\frac{1-2\mu^{*}}{r_{3}}\right).$
(16)

From the Hamiltonian function, given by Eq. 15, it is possible to determine
the spacecraft’s equations of motion in the rotating reference system, given
by:

$\ddot{x}-2\dot{y}=\Omega_{x},$ (17) $\ddot{y}+2\dot{x}=\Omega_{y},$ (18)
$\ddot{z}=\Omega_{z}.$ (19)

where $\Omega_{x}$ , $\Omega_{y}$ and $\Omega_{z}$ is the gradient of the
effective potential, as a function of $x$, $y$ and $z$, respectively.

The equations considered in this section, in particular the Hamiltonian
function (Eq.(15)) and the equations of motion given by Eqs.(17 \- 19), can be
used for the two-dimensional (tripole 2D) or three-dimensional (tripole 3D)
rotating mass tripole model.

In Section 6, we perform numerical simulations considering the gravitational
force of the body in which we want to place a spacecraft to orbit (asteroid or
comet), taking into account the solar radiation pressure (SRP) and the
gravitational attraction of the Sun. The equations of motion used in the
numerical simulation consider an inertial system whose origin of the reference
system coincides with the centre of mass of the asteroid. The equations of
motion used in the numerical simulation are shown next.

$\ddot{x}=-\frac{k\mu^{*}(x^{\prime}-x_{1}^{\prime})}{r^{\prime}{{}_{1}}^{3}}-\frac{k\mu^{*}(x^{\prime}-x_{2}^{\prime})}{r^{\prime}{{}_{2}}^{3}}-\frac{k(1-2\mu^{*})(x^{\prime}-x_{3}^{\prime})}{r^{\prime}{{}_{3}}^{3}}-P_{rad_{x}}\\\
+k_{sun}\frac{M_{sun}}{M_{ast}}\left(\frac{(x_{sun}-x^{\prime})}{r^{\prime
3}_{sun}}-\frac{x_{sun}}{r^{3}_{sun}}\right),$ (20)
$\ddot{y}=-\frac{k\mu^{*}(y^{\prime}-y_{1}^{\prime})}{r^{\prime}{{}_{1}}^{3}}-\frac{k\mu^{*}(y^{\prime}-y_{2}^{\prime})}{r^{\prime}{{}_{2}}^{3}}-\frac{k(1-2\mu^{*})(y^{\prime}-y_{3}^{\prime})}{r^{\prime}{{}_{3}}^{3}}-P_{rad_{y}}\\\
+k_{sun}\frac{M_{sun}}{M_{ast}}\left(\frac{(y_{sun}-y^{\prime})}{r^{\prime
3}_{sun}}-\frac{y_{sun}}{r^{3}_{sun}}\right),$ (21)
$\ddot{z}=-\frac{k\mu^{*}(z^{\prime}-z_{1}^{\prime})}{r^{\prime}{{}_{1}}^{3}}-\frac{k\mu^{*}(z^{\prime}-z_{2}^{\prime})}{r^{\prime}{{}_{2}}^{3}}-\frac{k(1-2\mu^{*})(z^{\prime}-z_{3}^{\prime})}{r^{\prime}{{}_{3}}^{3}}-P_{rad_{z}}\\\
+k_{sun}\frac{M_{sun}}{M_{ast}}\left(\frac{(z_{sun}-z^{\prime})}{r^{\prime
3}_{sun}}-\frac{z_{sun}}{r^{3}_{sun}}\right),$ (22)

where $k_{sun}$ is the force ratio of the Sun, $M_{sun}/M_{ast}$ is the mass
of the Sun in canonical unit, $x_{sun}$, $y_{sun}$ and $z_{sun}$ define the
position of the Sun on the axis $x$, $y$ and $z$, respectively. In the
simulations, we consider that the Sun is in the plane $xy$, consequently,
$z_{sun}$ = 0. The positions of the primary bodies $M_{1}$ , $M_{2}$ and
$M_{3}$, in the inertial frame, are given by the next equations, respectively.

$x^{\prime}_{1}=-L\cos{\Phi}\sin{\Psi}\cos{T}-(1-2\mu^{*})L\sin{\Phi}\sin{T},$
$y^{\prime}_{1}=-L\cos{\Phi}\sin{\Psi}\sin{T}+(1-2\mu^{*})L\sin{\Phi}\cos{T},$
(23) $z^{\prime}_{1}=L\cos{\Phi}\cos{\Psi},$

$x^{\prime}_{2}=L\cos{\Phi}\sin{\Psi}\cos{T}-(1-2\mu^{*})L\sin{\Phi}\sin{T},$
$y^{\prime}_{2}=L\cos{\Phi}\sin{\Psi}\sin{T}+(1-2\mu^{*})L\sin{\Phi}\cos{T},$
(24) $z^{\prime}_{2}=L\cos{\Phi}\cos{\Psi},$

$x^{\prime}_{3}=2\mu^{*}L\sin{\Phi}\sin{T},$
$y^{\prime}_{3}=-2\mu^{*}L\sin{\Phi}\cos{T},$ (25)
$z^{\prime}_{3}=-2\mu^{*}L\cos{\Phi}\cos{\Psi}/(1-2\mu^{*}).$

where $T$ is the time in canonical unity. The distances from the point mass
particle to the bodies $M_{1}$ , $M_{2}$ and $M_{3}$, in the inertial system
are, respectively, given by Eqs, 26 to 28.

$r_{1}=\sqrt{(x^{\prime}-x^{\prime}_{1})^{2}+(y^{\prime}-y^{\prime}_{1})^{2}+(z^{\prime}-z^{\prime}_{1})^{2}},$
(26)
$r_{2}=\sqrt{(x^{\prime}-x^{\prime}_{2})^{2}+(y^{\prime}-y^{\prime}_{2})^{2}+(z^{\prime}-z^{\prime}_{2})^{2}},$
(27)
$r_{3}=\sqrt{(x^{\prime}-x^{\prime}_{3})^{2}+(y^{\prime}-y^{\prime}_{3})^{2}+(z^{\prime}-z^{\prime}_{3})^{2}}.$
(28)

On the other hand, $r^{\prime}_{sun}$ and $r_{sun}$ is the Sun-spacecraft
distance and Sun-asteroid distance, respectively, given by Eqs. 29 and 30.

$r^{\prime}_{sun}=\sqrt{(x^{\prime}-x_{sun})^{2}+(y^{\prime}-y_{sun})^{2}+(z^{\prime}-z_{sun})^{2}},$
(29) $r_{sun}=\sqrt{x_{sun}^{2}+y_{sun}^{2}+z_{sun}^{2}},$ (30)

where $x_{sun}$, $y_{sun}$ and $z_{sun}$ define the position of the Sun in the
inertial system on the $x$, $y$ and $z$ axis, respectively. $P_{rad_{x}}$,
$P_{rad_{y}}$ and $P_{rad_{z}}$ represent the components $x$, $y$ and $z$,
respectively, of the acceleration due to the SRP (Pinheiro Lopes Masago et
al., 2016; Montenbruck et al., 2002; Beutler, 2005), given by Eq. 31.

$P_{rad}=C_{r}\frac{A}{m}P_{s}\frac{r_{0}^{2}}{r^{\prime 2}_{sun}}\hat{r},$
(31)

where $C_{r}$ is the radiation pressure coefficient. In the simulations we use
$C_{r}$ = 1.5. $m$ is the mass of the spacecraft and $A$ is the area of the
straight section of the spacecraft illuminated by the Sun. $P_{s}$ is the SRP
in the vicinity of the Earth-Sun distance and its numerical value is
approximately $4.55\times 10^{-6}$ N/m2; $r_{0}$ is the distance between the
Sun and the Earth and $\hat{r}$ is the radial unit vector of the Sun with
respect to the spacecraft. The value we adopted for the mass of the spacecraft
in the simulations was $m$ = 100 kilogrammes (kg) and the area of the solar
panel was $A=1.0$ m2, so $A/m$ = 0.01 m2/kg. Recalling that Equation 31 is
written in the asteroid-centred inertial frame.

##  3 Determination of the parameters for the symmetrical 2D and 3D triple
particle linkage model

Two different particle-linkage models were presented in the previous section
(2D and 3D tripole). Each particle-linkage model has unknown parameters, and
these parameters need to be determined to give a complete set of equations for
the model. In this study, considering the simplified model, the total mass and
angular velocity of the asteroids are equivalent to the actual values of those
same asteroids. Therefore, the 2D triple-particle linkage model only has four
unknown parameters, which are $\Phi$ [rad], $L$ [canonical unit], $k$
[dimensionless], and $\mu^{*}$ [canonical unit]. The three-dimensional triple-
particle linkage model is more complicated. It has five unknown parameters,
which are $\Phi$ [rad], $\Psi$ [rad], $L$ [canonical unit], $k$
[dimensionless] and $\mu^{*}$ [canonical unit].

The unknown parameters are determined using the locations of the external
equilibrium points obtained from the Mascon approach, similar to the approach
used by Zeng et al. (2015, 2016a); Yang et al. (2018). The main idea is to
find the parameters that generate equilibrium points as close as possible to
the ones obtained by using the Mascon model. The locations of the equilibrium
points [$\hat{x}_{E},\hat{y}_{E},\hat{z}_{E}$] can be determined as shown by
Eq. 32. This method can be used for asteroids that have external equilibrium
points. On the other hand, there are some asteroids do not have external
equilibrium points, for example, the asteroid 1998 KY26. For this type of
asteroid, the method that minimises the errors of the effective potential or
gradient of the effective potential can be used (Zeng et al., 2018).

$\Omega_{x}(\hat{x}_{E},\hat{y}_{E},\hat{z}_{E})=\Omega_{y}(\hat{x}_{E},\hat{y}_{E},\hat{z}_{E})=\Omega_{z}(\hat{x}_{E},\hat{y}_{E},\hat{z}_{E})=0.$
(32)

The idea is to find the parameters of the tripole model that minimise the
differences of all external equilibrium points between the simplified and the
polyhedron model. The method used to determine the parameters for the
two/three-dimensional triple-particle-linkage model is described below. The
optimisation variables for the simplified model are
$X=[\Phi,\Psi,L,k,\mu^{*}$]. Before the optimisation, we define the
constraints for each variable using the lower limits
[$\Phi_{min},\Psi_{min},L_{min},k_{min},\mu^{*}_{min}$] and the upper limits
[$\Phi_{max},\Psi_{max},L_{max},k_{max},\mu^{*}_{max}$] for those parameters.
For the nonlinear optimisation problem, the performance index is defined as

$\mathbf{J}(\Phi,\Psi,L,k,\mu^{*})=\\\
\sum_{i=1}^{n}\sqrt{(\hat{x}_{T_{i}}d^{*}-x_{P_{i}})^{2}+(\hat{y}_{T_{i}}d^{*}-y_{P_{i}})^{2}+(\hat{z}_{T_{i}}d^{*}-z_{P_{i}})^{2}},$
(33)

where [$\hat{x}_{T_{i}},\hat{y}_{T_{i}},\hat{z}_{T_{i}}$] are the positions,
in canonical units, of the equilibrium points obtained by the simplified
models (tripole 2D or 3D). Note that $d$ can be obtained from Eq. 10. On the
other hand, [$x_{P_{i}},y_{P_{i}},z_{P_{i}}$] represent the locations of the
equilibrium points (in metres) obtained by the Mascon model. The index i
corresponds to the i-th equilibrium point and, finally, $n$ is the total
number of external equilibrium points. The equality ($c_{eq}=0$) constraints
for tripole 2D and 3D models are shown in Eqs. 34 and 35, respectively,

$c_{eq}=\begin{bmatrix}m_{1}+m_{2}+m_{3}-1&\\\ L\cos{\Phi}-1/2&\end{bmatrix},$
(34) $c_{eq}=\begin{bmatrix}m_{1}+m_{2}+m_{3}-1&\\\
L\cos{\Phi}\sin{\Psi}-1/2&\end{bmatrix}.$ (35)

Mathematically, we can write the performance index as a restricted
minimisation, as shown in Eq. 36.

$min~{}\mathbf{J}(\Phi,\Psi,L,k,\mu^{*})~{}such~{}that~{}c_{eq}=0,$ (36)

where $\mathbf{J}(\Phi,\Psi,L,k,\mu^{*})$ is a function that returns a scalar
value and $c_{eq}$ are functions that return vectors. Here we used nonlinear
optimisation routines developed in Matlab to find the optimal solutions that
minimise $\mathbf{J}$. The optimisation problem is defined and solved with a
nonlinear programming method (NLP).

##  4 Application to realistic elongated asteroids

In this section, we apply the particle-linkage models mentioned in the
previous sections to three realistic elongated asteroids. The asteroids are:
1620 Gegraphos, 433 Eros and 243 Ida. The optimisation method is used to
determine the parameters of the simplified model, as shown in this section.
After that, we demonstrate the improved performance of the symmetric three-
dimensional triple-particle-linkage model with respect to the two-dimensional
axisymmetric triple-particle-linkage model.

###  4.1 Parameters of the sample elongated asteroids

The physical and orbital parameters of the asteroids 1620 Geographos, 433 Eros
and 243 Ida were obtained from Rozitis & Green (2014); Carry (2012) and Baer
et al. (2011), respectively, and can be found in Table 2. In this study, we
also applied the Mascon gravity framework using a shaped polyhedral source to
precisely represent the gravitational field of these asteroids (Chanut et al.,
2015a). We recommend the readers to review the details of this method in
Venditti et al. (2013); Aljbaae et al. (2017).

Table 1: Physical and Mascon Parameters of 1620 Gegraphos, 433 Eros and 243
Ida.

Asteroid | bulk density | Rotation Period | M | Verticies $\&$ Facets
---|---|---|---|---
| ($g~{}cm^{-3}$) | (hours) | (kg) |
1620 Geographos | 2.15 $\pm$ 0.5 | 5.2233 | 1.65$\times 10^{13}$ | 1022 $\&$ 2040
433 Eros | 2.67 $\pm$ 0.03 | 5.2702 | 6.69$\times 10^{15}$ | 856 $\&$ 1708
243 Ida | 2.35 $\pm$ 0.29 | 4.63 | 3.82$\times 10^{16}$ | 1022 $\&$ 2040

Asteroid |  $E_{1}$ (km) |  $E_{2}$ (km) |  $E_{3}$ (km) |  $E_{4}$ (km)
---|---|---|---|---
1620 Geographos |  [2.6318, 0.1782, 0.0046] |  [- 0.0352, 1.9315, 0.0017] |  [- 2.6812, 0.2003, -0.0039] |  [- 0.0196, -1.9698, 0.0019]
433 Eros |  [19.1246, -2.5756, 0.1449] |  [0.4637,14.7338, - 0.06604] |  [-19.6711, - 3.2861, -0.1231] |  [-0.4469, -13.991, -0.0791]
243 Ida |  [29.7950, -2.5972, 0.6870] |  [-0.5688, 25.9245, -0.1113] |  [-30.3092, - 1.8763, 0.3958] |  [-0.4906, - 25.7145, -0.0993]

Table 2: Positions of Equilibrium Points for 1620 Gegraphos, 433 Eros and 243
Ida.

After adapting the shape of each asteroid in our study according to the
procedure presented in Chanut et al. (2015b), we used the Mascon approach to
precisely calculate the positions of the equilibrium points. Table 2 provides
the results.

We recall that we will only consider the external equilibrium points, due to
the fact that the internal equilibrium points have no physical meaning.

To perform the simulations, it is necessary to determine the boundary
restrictions of each asteroid under study using the particle-linkage model.
They are chosen as follows:

  1. 1.

3D triple-particle-linkage model:

For each asteroid, [$\Phi_{min}$, $\Phi_{max}$] are set to [-0.5, 0.5] rad,
[$\Psi_{min}$, $\Psi_{max}$] are set to [1.39626, 1.91986], [$L_{min}$,
$L_{max}$] are set to [0, 2], [$k_{min}$, $k_{max}$] are set to [0, 9] and
[$\mu^{*}_{min}$, $\mu^{*}_{max}$] are set to [0.001, 0.999]

  2. 2.

2D triple-particle-linkage model:

For each asteroid, [$\Phi_{min}$, $\Phi_{max}$] are set to [-0.5, 0.5] rad,
[$L_{min}$, $L_{max}$] are set to [0, 2], [$k_{min}$, $k_{max}$] are set to
[0, 9] and [$\mu^{*}_{min}$, $\mu^{*}_{max}$] are set to [0.001, 0.999]

The boundary constraints of the geometric parameters $\Phi$, $\Psi$, $L$, are
determined based on the actual shapes of the asteroids. Due to the geometric
shape of the body and the information we obtain from the equilibrium points
from the high fidelity model, we observed that the equilibrium points $E_{1}$
and $E_{3}$ do not have high values on the $y$ axis, making us to consider the
azimuth angle varying from -$30^{\circ}$ to + $30^{\circ}$. Observing the
values of the equilibrium points on the $z$ axis, it is possible to notice
that these asteroids do not have a very high mass distribution on the $z$
axis, therefore, we consider the elevation angle varying between -$20^{\circ}$
and + $20^{\circ}$. Finally, knowing that, by the theoretical definition of
our study, the distance between $M_{1}$ and $M_{2}$ is one canonical unit, we
concluded that the value of $L$ is not much above the numerical value.

The boundary constraints of the physical parameters, $k$ and $\mu^{*}$, were
determined as follows. Through the articles in the literature that use the
particle-linkage model (Lan et al., 2017; Yang et al., 2018; Zeng et al.,
2015), we observed that the value of $k$ found in asteroids varies between 0.4
and 9, so we consider an initial guess between 0 and 9. On the other hand, the
contour restrictions of the masses were defined based on their theoretical
intervals.

After the optimisation has been carried out, if the value of one of the
parameters is close to the contour limit, we modify its range and carry out
the optimisation again. This procedure is carried out until the values of the
optimised parameters are located between the specified limits. The initial
guesses for the 2D and 3D model, were chosen as follows:

  1. 1.

3D triple-particle-linkage model:

For 1620 Geographos is [$\Phi$, $\Psi$, $L$, $k$ , $\mu^{*}$] = [0.36, 1.5, 1,
0.2, 0.24], 433 Eros is [$\Phi$, $\Psi$, $L$, $k$ , $\mu^{*}$] = [0.3, 1.5,
0.5, 0.2, 0.28] and for 243 Ida is [$\Phi$, $\Psi$, $L$, $k$ , $\mu^{*}$] =
[0.3, 1.5, 1, 0.4, 0.2].

  2. 2.

2D triple-particle-linkage model:

For 1620 Geographos is [$\Phi$, $\Psi$, $L$, $k$ , $\mu^{*}$] = [0.36,
$\pi/2$, 1, 0.2, 0.24], 433 Eros is [$\Phi$, $\Psi$, $L$, $k$ , $\mu^{*}$] =
[0.3, $\pi/2$, 0.5, 0.2, 0.28] and 243 Ida is[$\Phi$, $\Psi$, $L$, $k$ ,
$\mu^{*}$] = [0.3, $\pi/2$, 1, 0.4, 0.2].

The optimal values of $\Phi$, $\Psi$, $L$, $k$ and $\mu^{*}$ parameters
obtained using the optimisation method mentioned above are shown in Tables 4
and 4, where we see the performance index $J_{1}$ and $J_{2}$. These
performance indices are calculated after the optimal parameters $J_{0}$ are
obtained. These performance indexes, $J_{1}$ and $J_{2}$, provide the relative
errors, maximum and minimum, respectively, for a single equilibrium point and
are defined as follows

$\mathbf{J_{max}}(\Phi,\Psi,L,k,\mu^{*})=\\\
max\left(\frac{\sqrt{(\hat{x}_{T_{i}}d^{*}-x_{P_{i}})^{2}+(\hat{y}_{T_{i}}d^{*}-y_{P_{i}})^{2}+(\hat{z}_{T_{i}}d^{*}-z_{P_{i}})^{2}}}{L}\right)\\\
\times 100\%,~{}~{}~{}~{}i=1,2,3,4,$ (37)

$\mathbf{J_{min}}(\Phi,\Psi,L,k,\mu^{*})=\\\
min\left(\frac{\sqrt{(\hat{x}_{T_{i}}d^{*}-x_{P_{i}})^{2}+(\hat{y}_{T_{i}}d^{*}-y_{P_{i}})^{2}+(\hat{z}_{T_{i}}d^{*}-z_{P_{i}})^{2}}}{L}\right)\\\
\times 100\%,~{}~{}~{}~{}i=1,2,3,4.$ (38)

Table 3: Optimisation Results for the tripole 3D model.

Asteroid | $\Phi$ | $\Psi$ | $L$ | $k$ | $\mu^{*}$ | J | $J_{max}$ | $J_{min}$
---|---|---|---|---|---|---|---|---
| [deg] | [deg] | [c. u.] | [adim] | [c. u.] | [km] | % | %
1620 Geographos | 7.1780 | 89.773 | 0.5039 | 0.2991 | 0.2024 | 0.120 | 1.5342 | 0
433 Eros | -19.892 | 88.7891 | 0.5318 | 0.5195 | 0.2815 | 1.743 | 3.313 | 0
243 Ida | -6.3105 | 87.2621 | 0.5036 | 0.6529 | 0.1636 | 2.005 | 2.975 | 0

Table 4: Optimisation Results for the tripole 2D model.

Asteroid | $\Phi$ | $\Psi$ | $L$ | $k$ | $\mu^{*}$ | J |  $J$max |  $J$min
---|---|---|---|---|---|---|---|---
| [deg] | [deg] | [c. u.] | [adim] | [c. u.] | [km] | % | %
1620 Geographos | 7.0884 | 90 | 0.5038 | 0.2843 | 0.1971 | 0.147 | 1.348 | 0.373
433 Eros | -18.935 | 90 | 0.5286 | 0.4454 | 0.2619 | 2.144 | 3.529 | 1.359
243 Ida | -7.3614 | 90 | 0.5041 | 0.7169 | 0.1672 | 5.020 | 7.438 | 3.125

##  5 Comparison between the simplified and Mascon model

We have verified that the positions of the equilibrium points in the
simplified model compared to the Mascon model are closer, giving the first
evidence of the validation of the simplified model. Next, we analyse the
classifications of the equilibrium points of the three asteroids studied in
this paper and we found that our analysis is consistent with the analysis made
using the Mascon model. After that, we compared the relative errors between
the potential estimated by the 3D tripole model and the Mascon approach
(Chanut et al., 2015a), to justify the simplified model.

Figures 3, 4 and 5 show the relative errors of the gravitational potentials
between the two models (simplified and Mascon model) for the asteroid 1620
Geographos, 433 Eros and 243 Ida, respectively. The black vertical line refers
to the surface of the asteroid. The left side of the black vertical line is
the relative error of the potential inside the asteroid, and has no meaning
for the purpose of this work. On the other hand, outside the asteroid, we
observe the relative error of the potential when we consider the asteroid as a
mass point (red), 2D tripole model (green) and 3D tripole model (blue).

We see in Fig. 3 that, when the spacecraft is close to the asteroid 1620
Geographos (less than 10 km from centre), we cannot consider the asteroid as a
point of mass, due to the high relative error between the gravitational
potentials. On the other hand, when we model the asteroid as a rotating mass
tripole (2D or 3D model), the result shows a good agreement between the
models. The 3D tripole model has almost the same accuracy as the 2D tripole
model, due to the fact that the mass distribution of this asteroid is
predominantly on the $xy$ axis.

Figure 3: Relative error of the gravitational potential for the asteroid 1620
Geographos.

Figure 4 shows the relative error between the potentials for the asteroid 433
Eros. In contrast to 1620 Geographos, the 3D tripole model is more accurate
than the 2D version, when the spacecraft is less than 60 km from the centre of
the asteroid.

Figure 4: Relative error of the gravitational potential for the asteroid 433
Eros.

Finally, Fig. 5 shows similar behaviour as in the case of 433 Eros. The 3D
tripole model for the asteroid 243 Ida is more accurate than the 2D one.

Figure 5: Relative error of the gravitational potential for the 243 Ida
asteroid.

These results indicate that the 3D tripole model approximates with good
accuracy the external potential of some asteroids with elongated bodies.
Besides that, it is worth to mention that the tripole model can approximate
external potentials for 1030225 points uniformly distributed in the $xy$ plane
around an asteroid in about 3 seconds, using an Intel processor of 3.6 GHz. It
is about 200 times faster than the Mascon approach shown in Chanut et al.
(2015a), which takes about 10 min. Thus, for a first analysis of a space
mission whose objective is to visit an asteroid, it is possible to use the 3D
tripole model to obtain preliminary results, since this model is close to the
high precision model (Mascon model) and demands well less computational time
to perform the simulations.

##  6 Numerical investigations

A spacecraft that is close to the surface of an asteroid undergoes complex
perturbation, which may lead it to collide with the asteroid or even escape
from the domain of its gravitational field. Understanding these perturbations
is extremely important to determine regions in space where it is possible to
find natural orbits around asteroids. These regions allow a spacecraft to
orbit an asteroid for a long time. We will call these regions “bounded
regions”. In contrast, the regions around the asteroids where no natural orbit
is possible, will be called “unbounded regions”.

The present section was developed with the objective of finding possible
regions in the vicinity of the asteroids 1620 Geographos, 433 Eros and 243 Ida
capable of maintaining a spacecraft for a long time without the need for
orbital manoeuvres.

In the simulations performed here, we considered the gravitational field of
the asteroid (modelled as a three-dimensional rotating mass tripole) in which
it is desired to orbit, the gravitational field of the Sun (modelled as a mass
point), and the SRP. The results are shown in the inertial frame with an
origin coinciding with the centre of mass of the asteroid under analysis. In
the simulations where we take into account the SRP we assume that the area-to-
mass ratio is 0.01 $m^{2}/kg$. The methodology adopted and the results
obtained are described next.

###  6.1 Numerical simulations methodology

To determine the initial conditions of the spacecraft to start the
simulations, we initially consider the location of the spacecraft in terms of
the Keplerian elements $\mathbf{\kappa}=[a,e,i,\Omega,\omega,f]$. To perform
the numerical simulations and build the initial condition grids, we start the
orbit with eccentricity 0 and vary it until eccentricity = 0.99, in steps of
0.01. The spacecraft’s semi-major axis varies from an initial value $a_{i}$ to
a final value $a_{f}$, whose values depend on the asteroid that we are
analysing. In this study, we investigated direct orbits ($i_{0}=0^{\circ}$)
and retrograde orbits ($i_{0}=180^{\circ}$). The maximum integration time is
365 days. We use the Runge-Kutta 7/8 method with a time interval, in canonical
unit, of 0.01. We consider the initial Keplerian elements as shown in Table 5.
The acronym RAAN in the table means Right Ascension of the Ascending Node.

Table 5: Initial Keplerian elements. Initial semi-major axis [m] | $a_{i}\in[a_{initial},~{}a_{final}]$
---|---
Initial eccentricity [adim] |  $e_{0}\in$ [0, 0.99]
Initial inclination [deg] |  $i_{0}\in$ [0 and 180]
Initial RAAN [deg] | $\Omega_{0}=0$
Initial periselene anomaly [deg] | $\omega_{0}=0$
Initial true anomaly [deg] | $f_{0}=0$

In this initial consideration, we assume that there is only the asteroid that
we wish to orbit and the spacecraft. The asteroid, initially, is modelled as a
mass point. Based on the initial conditions (the osculating orbit of the
spacecraft), we obtain the initial conditions of position ($\mathbf{\rho_{0}}$
= ($\mathbf{[x^{\prime}_{0},~{}y^{\prime}_{0},~{}z^{\prime}_{0}}]^{T}$) and
velocity ($\mathbf{\eta_{0}}$ =
($\mathbf{[\dot{x}^{\prime}_{0},~{}\dot{y}^{\prime}_{0},~{}\dot{z}^{\prime}_{0}}]^{T}$)
of the spacecraft in the Cartesian system, given by $\mathbf{X_{0}}$ =
[$\mathbf{\rho_{0}},~{}\mathbf{\eta_{0}}$], to be used in the equations of
motion. $\mathbf{\rho_{0}}$ and $\mathbf{\eta_{0}}$ define the state vector of
the spacecraft at $t$ = 0 in the restricted problem, where the vector
$\mathbf{X_{0}}$ $\in$ $\mathbb{R}^{6}$, $\mathbf{\rho_{0}}$ $\in$
$\mathbb{R}^{3}$ (initial position) and $\mathbf{\eta_{0}}$ (initial velocity)
$\in$ $\mathbb{R}^{3}$. The initial conditions obtained, in the Cartesian
plane, can be written as $\mathbf{X_{0}}$ = [$x^{\prime}_{0}$ 0, 0, 0,
$\pm~{}\dot{y}^{\prime}_{0}$, 0]T. We assign + $\dot{y}^{\prime}_{0}$ when the
orbit is direct ($i_{0}$= 0∘) and - $\dot{y}^{\prime}_{0}$ for retrograde
orbits ($i_{0}$ = 180∘).

Once the initial conditions of the spacecraft are found with respect to the
asteroid’s centre of mass, in the cartesian system, the numerical integration
is carried out including the spacecraft, the shape of the asteroid (3D
tripole) and the perturbation of the Sun. Regarding the perturbation of the
Sun, in some simulations, we consider only the solar gravity perturbation,
while in others we consider both, the solar gravity perturbation and the SRP.
This is done to represent the situation where the area-to-mass ratio is very
small. Initially, ($t$ = 0) the three bodies (Asteroid, spacecraft and Sun)
are aligned, as shown in Fig. 6 (out of scale). Every simulation is performed
in the asteroid-centred inertial frame. To determine the initial conditions of
the Sun, we assign the value of the semi-major axis and eccentricity according
to references Rozitis & Green (2014); Carry (2012); Baer et al. (2011) for the
1620 Geographos, 433 Eros and 243 Ida, respectively. The other Keplerian
elements are considered equal to zero. It means that, initially, the asteroid
is located at the Pericentre of its orbit around the Sun. We chose this
starting position because the asteroid is in a region close to Earth. After
that, we convert the initial Keplerian elements of the asteroid to Cartesian
elements and use them in Equations 20, 21 and 22.

Figure 6: Graphic representation of the system studied.

During the numerical integrations, we monitor the particles that survived
during the integration time, the particles that collided with the asteroid and
the particles that were ejected from the system. In the dynamics of the
spacecraft, a collision with the asteroid is considered to occur when the
spacecraft crosses the disk with a mean radius of the asteroid. We consider an
ejection when the spacecraft’s distance from the asteroid reaches a value
greater than the Hill radius of the asteroid (Murray & Dermott, 2000).

The results are presented in graphs that show the lifetime of the spacecraft
that survives, for each set of initial conditions, with respect to the semi-
major axis and eccentricity of the initial orbit. The regions in which the
spacecraft survives for 1 year are called “bounded regions”. The regions in
which the spacecraft survives for a few days are called “unstable regions”.

Figure 7(a) shows direct orbits and Fig. 7(b) retrograde orbits taking into
account the shape of the asteroid 433 Eros and the solar gravity perturbation.
They provide the duration of the orbits as a function of the initial semi-
major axis and eccentricity of the particle around 433 Eros. The colour code
provides the time the particle remains in orbit for each initial condition.
The white area in Fig. 7 presents regions where the initial conditions of the
spacecraft are inside the asteroid, so it has no physical significance. We
investigated the orbital dynamics of a spacecraft with
$40~{}\leq~{}a_{0}~{}\leq~{}1000$ km around 433 Eros asteroid. These values of
$a_{0}$ were chosen because it is the region where the relative error of the
gravitational potential of the rotating triple-particle-linkage three-
dimensional model is similar to the potential using the Mascon model, as shown
in Fig. 4. Note that the value of $a_{f}$ is below the Hill radius ($\approx$
3270 km) of this asteroid. That allows us to safely neglect the gravitational
effects of other celestial bodies. The same reasoning was used to find the
range of the semi-major axis of the asteroids studied in this work (1620
Gegraphos and 243 Ida).

It is observed, in Fig. 7, that the solar gravity perturbation has little
effect on the dynamics of a spacecraft in the vicinity of asteroid 433 Eros,
while the non-spherical gravitational field near the asteroid significantly
affects the motion of a particle.

A similar result was found in Ni et al. (2014), in which the author showed
that solar gravity perturbation is not strong enough to modify families of
periodic orbits around the asteroid 433 Eros. In addition, the disturbance of
the solar gravity also does not change the stability of the periodic orbits
obtained here.

Similar results were found for asteroid 1620 Geographos and 243 Ida. The next
results take into account the shape of the asteroid (3D tripole), the solar
gravity perturbation and the SRP.

(a) Lifetimes in the regions close to 433 Eros without SRP (direct orbits).

(b) Lifetimes in the regions close to 433 Eros without SRP (retrograde
orbits).

Figure 7: Plots of $a_{0}$ versus $e_{0}$ showing the evolution of the
lifetimes in the regions close to 433 Eros asteroid (direct orbits 7(a) and
retrograde orbit 7(b)) without SRP.

Figure 8 provides the initial condition grids in the vicinity of the asteroid
1620 Geographos, where Fig. 8(a) shows direct orbits and Fig. 8(b) retrograde
orbits. Note that the spacecraft cannot naturally orbit the asteroid 1620
Geographos for a period of 365 days. The maximum duration is around 10 days,
but this time is long enough to make observations around the asteroid, so the
orbits found here can have practical applications. These results occur due to
the small mass of the asteroid 1620 Geographos, which makes the SRP to become
significant, making the spacecraft to collide with the asteroid or escape the
sphere of its influence.

Note that there is a narrow strip (NS) in Fig. 8, which initially appears at
$a$ $\approx$ 11 km and $e$ = 0 and extends to $a$ = 30 km and $e$ $\approx$
0.9. These regions exist both for direct and retrograde orbits. The maximum
survival time of these orbits around the asteroid 1620 Geographos is
approximately 10 days. For times longer than this, the spacecraft will escape
the system or collide with the asteroid. It is possible to observe that the
orbits survive from 4 to 8 days in the regions on the right side of the NS,
while the orbits in the region at the left of the NS survive from 0 to 3 days,
approximately. The radius of the Pericentre ($r_{p}=a(1-e))$ coinciding with
the position of the narrow strip (NS) are called here $r_{p_{NS}}$.

The regions where $r_{p}~{}<~{}r_{p_{NS}}$ is close to the asteroid. In this
region, the asteroid’s gravitational field is strong, making the spacecraft to
be attracted by the asteroid or perform a swing-by and quickly escape from the
system. On the other hand, in regions where $r_{p}~{}>~{}r_{p_{NS}}$, the
perturbation of the Sun, especially the SRP, dominates the dynamics making the
spacecraft to escape from the system in a few days.

(a) Lifetimes in the regions close to Geographos (direct orbits).

(b) Lifetimes in the regions close to Geographos (retrograde orbits).

Figure 8: Plots of $a_{0}$ versus $e_{0}$ showing the evolution of the
lifetimes in the regions close to 1620 Geographos asteroid (direct orbits 8(a)
and retrograde orbits 8(b)).

Figures 9(a) and 9(b) provide information about the final destination of the
orbits after the integration time shown in Figs. 8(a) and 8(b), respectively.
The yellow regions indicate that the spacecraft has escaped the sphere of
influence of the asteroid. The blue regions show the initial conditions that
lead the spacecraft to collide with the asteroid. Proving our argument, note
that, for orbits that start close to the Geographos surface, the spacecraft is
captured, ending in collision with the asteroid. Orbits with
$r_{p}~{}>~{}r_{p_{NS}}$ makes the particles to eject from the system. Observe
that the boundary between the blue and yellow regions is the narrow strip of
1620 Geographos. It is this boundary that defines whether the orbit will be
captured or ejected from the system. Overall, we see that a spacecraft can not
orbit the asteroid 1620 Geographos for a long period.

(a) Collision region for direct orbits around Geographos.

(b) Collision region for retrograde orbits around Geographos.

Figure 9: Direct orbits 9(a) and retrograde orbit 9(b) around Geographos.
Collision regions with the asteroid (blue) and regions of ejection from the
system (yellow). The white region represents the locations with initial
conditions inside the asteroid.

We also make an analysis considering 433 Eros. Some previous studies have
analysed the orbital dynamics around this asteroid, neglecting the Sun’s
disturbances or considering only the solar gravity perturbation, due to the
gravitational domain of the asteroid (Scheeres et al., 2000; Chanut et al.,
2014; Ni et al., 2014). To carry out an investigation complementary to the
results available in the literature, we investigated the orbital dynamics of a
spacecraft with $40$ km $\leq~{}a_{0}~{}\leq~{}1000$ km, taking into account
the effect of the solar gravity and SRP. Our results are presented in Fig.
10(a) for direct orbits and Fig. 10(b) for retrograde orbits. We see that
there is a narrow strip (starting near $a_{0}$ $\approx$ 300 km and
$e_{0}~{}=~{}0$) where the spacecraft survives around 30 days. This region
corresponds to the $r_{p_{NS}}$ of the 433 Eros.

Note that, for direct orbits, when $r_{p}$ $<$ $r_{p_{NS}}$, the movement of
the particle around the asteroid 433 Eros is unstable due to the shape of 433
Eros and the perturbations from the Sun. Observe that, making a small
variation in the initial conditions, it generates completely different
trajectories for the particle, reflecting on the lifetime of the orbits and,
consequently, in the final destination of the particle. On the other hand,
there are regions with solutions that remain bounded close to the asteroid 433
Eros when we consider retrograde orbits (see Fig. 10(b)).

Space dust most commonly exists around asteroids and planets in direct orbits.
On the other hand, retrograde orbits are unlikely to arise naturally. Thus,
regions where retrograde orbits survive and direct orbits do not survive are
great options for placing a spacecraft, due to the fact that, in these
regions, the probability of having space dust is low, decreasing the risk of a
spacecraft colliding with some dust particle (Araujo et al., 2015; Araujo et
al., 2017).

(a) Lifetimes in the regions close to 433 Eros (direct orbits).

(b) Lifetimes in the region closes to 433 Eros (retrograde orbits).

Figure 10: Plots of $a_{0}$ versus $e_{0}$ showing the evolution of the
lifetimes in the regions close to 433 Eros asteroid (direct orbits 10(a) and
retrograde orbits 10(b)).

Note that, when we consider the SRP, as presented in this article, the bounded
orbits that once existed around 433 Eros cease to exist. This fact shows that
the effect of the perturbation of SRP is much larger than the solar gravity
perturbation. It happens because the solar gravity perturbation involves
direct and indirect terms, that is, the difference between the action of the
Sun in the spacecraft and in the asteroid, as shown by Eqs. 20 \- 22. Those
differences are small when the spacecraft gets closer to the asteroid ($\sim$
$10^{-7}$). In contrast, the SRP acts directly on the vehicle, with no
indirect term to compensate it, and this makes the effect of SRP to be much
larger ($\sim$ $10^{-4}$). There is no doubt that, if we decrease the area
mass ratio $A/m$, or if we use the lowest reflectivity coefficient ($C_{r}$),
we can find an increasingly larger regions where the spacecraft remains in
orbit around the asteroid for longer times.

Figures 11(a) and 11(b) provide information about the final destination of the
orbits after the integration time shown in Figs. 10(a) and 10(b). Note that,
for direct orbits, when $r_{p}~{}<~{}r_{p_{NS}}$ km, the spacecraft collides
with the asteroid or escapes from the sphere of influence of the asteroid due
to the complex dynamics of the movement. For retrograde orbits, we see that
there are orbits close to the surface of the asteroid (in black) that survives
through numerical integration, that is, that remain bounded.

(a) Collision region for direct orbits for 433 Eros.

(b) Collision region for retrograde orbits for 433 Eros.

Figure 11: Direct orbits 11(a) and retrograde orbits 11(b) around 433 Eros.
Collision region with the asteroid (blue), orbits that survived for 365 days
(black) and regions of ejection from the system (yellow). The white region
represents the locations with initial conditions inside the asteroid.

Through Figs. 10(a) and 10(b), we see that, when $r_{p}~{}>~{}r_{p_{NS}}$, for
different initial conditions, the spacecraft’s lifetime is similar. Note that
for small variations of the initial conditions, the spacecraft has
approximately the same lifetimes, and probably the same final destination.
This statement can be seen in Figs. 11(a) and 11(b). Note that, for
$r_{p}~{}>~{}r_{p_{NS}}$, all orbits are ejected from the system (yellow).
This happens because, as we move away from the central body, the disturbing
effect of the Sun becomes larger and, consequently, making the spacecraft (or
dust particle) to escape the asteroid’s sphere of influence.

Finally, grids of initial conditions were used to investigate the motion of a
spacecraft around 243 Ida, as shown in Figs. 12(a) and 12(b). The orbital
elements used to build the grids of initial conditions shown in Figure 12 were
identical to those used for asteroid 433 Eros. Note that, for direct orbits
close to the asteroid ($50$ km $<~{}a_{0}~{}<$ 108 km) and low eccentricity,
it is possible to observe solutions that remain bounded around asteroid 243
Ida. It was precisely in this region that a moon was observed orbiting the
asteroid 243 Ida, called Dactyl, which is a proof that the technique used here
works very well.

(a) Lifetimes in the regions close to 243 Ida (direct orbits).

(b) Lifetimes in the regions close to 243 Ida (retrograde orbits).

Figure 12: Plots of $a_{0}$ versus $e_{0}$ showing the evolution of the
lifetimes in the regions close to 243 Ida asteroid (direct orbits 12(a) and
retrograde orbits 12(b)).

On the other hand, there is a large region with solutions that remain bounded
around 243 Ida when we consider retrograde orbits and $a_{0}~{}<~{}200$ km, as
shown in Fig. 12(b) (Sanchez & Prado, 2019). Note that this region exist even
when the eccentricity is high.

(a) Collision region for direct orbits for 243 Ida.

(b) Collision region for retrograde orbits for 243 Ida.

Figure 13: Direct orbits 13(a) and retrograde orbit 13(b) for 243 Ida
asteroid. Collision region with the asteroid (blue), orbits that survived for
365 days (black) and regions of ejection from the system (yellow). The white
region represents the locations with initial conditions inside the asteroid.

Figure 13 provides information about the final destination of the orbits after
the integration time shown in Fig. 12. From the information shown in Fig.
13(a), it is possible to conclude that the direct orbits around asteroid 243
Ida are unstable when $a_{0}~{}>~{}120$ km. Note that small variations in the
initial conditions make the final destination of the spacecraft to change
completely, causing the spacecraft to collide or move away from the asteroid.

For retrograde orbits, note that, when $200$ km $<~{}a_{0}~{}<$ 850 km, most
of the orbits collide with asteroid 243 Ida. There are few orbits that eject
from the system when $a_{0}~{}\sim~{}700$ km. When $a_{0}~{}>~{}850$ km, due
to the low gravitational force of 243 Ida, the orbits tend to escape from the
sphere of influence of the asteroid.

Similarly to the study made for the asteroid 433 Eros, it is possible to find
regions around asteroid 243 Ida in which for direct orbits the regions are
unstable, while for retrograde orbits the regions are bounded, those are
interesting locations to place a spacecraft, as already explained.

##  7 Conclusion

In this work, we developed a simplified model to represent the gravity field
of elongated asteroids (convex or not) with spatial mass distribution,
inspired by the existing particle-linkage model. The purpose of this
simplified model is to use a three-dimensional axisymmetric triple-particle-
linkage model that consists, essentially, of three-point particles distributed
in space and two rigid rods with negligible mass. Compared to the mass point
model, it was showed that there are ranges of particle-asteroid distances
where the accuracy of the proposed model is better, but it still keeps a small
processing time compared to high fidelity models. The simplified model was
then been applied to the asteroids 1620 Geographos, 433 Eros and 243 Ida. The
results indicate that the three-dimensional rotating mass tripole has
advantages over the two-dimensional rotating tripole model in terms of
accuracy, and get results close to the high fidelity models, using much less
computer time.

The results focused on the dynamical study of a spacecraft around the
asteroids 1620 Geographos, 433 Eros and 243 Ida, for direct and retrograde
orbits, considering an area/mass ratio $A/m~{}=~{}0.01$ m2/kg. We found
solutions that remain bounded in the vicinity of these asteroids. Orbits that
remain bounded are important, because they are adequate locations to place a
spacecraft. Furthermore, they indicate good places to search for a new member
of the asteroid system, as exemplified by the asteroid 243 Ida.

It was observed that no orbit survives for 365 days around the asteroid
Geographos when we consider the SRP. On the other hand, we found orbits that
remain bounded around the 433 Eros and 243 Ida. Some solutions exist for
direct orbits, while other solutions exist for retrograde orbits. Considering
the asteroid 433 Eros, the orbits that remain bounded exist only for
retrograde orbits. On the other hand, when we investigated the asteroid 243
Ida, we found orbits that remain bounded to both direct and retrograde senses.
Determining these regions is important when it comes to astronautical
applications, because regions in space where the orbits remain bounded for
retrograde orbits and that escape from asteroid for direct orbits are adequate
to place a spacecraft, due to the fact that there is a low risk of collisions
with dust particles and they provide stability for the orbit, reducing the
orbital manoeuvres required to keep the spacecraft orbiting the asteroid. So,
simplified models can be used to assist in the pre-designing of real missions.

Finally, It is worth to mention that the small size of our targets allow us to
safely neglect the shadowing phenomenon in the SRP. However, Future
applications of this model could involve the study of polar orbits considering
the effects of the shadowing phenomenon in the SRP.

##  8 Acknowledgements

The authors wish to express their appreciation for the support provided by:
grants 140501/2017-7, 150678/2019-3, 422282/2018-9 and 301338/2016-7 from the
National Council for Scientific and Technological Development (CNPq); grants
2016/24561-0 and 2016/18418-0, from São Paulo Research Foundation (FAPESP);
grant 88887.374148/2019-00 from the National Council for the Improvement of
Higher Education (CAPES) and to the National Institute for Space Research
(INPE).

## Data Availability

The data underlying this article will be shared on reasonable request to the
corresponding author

## References

  * Aljbaae et al. (2017) Aljbaae S., Chanut T. G. G., Carruba V., Souchay J., Prado A. F. B. A., Amarante A., 2017, mnras, 464, 3552
  * Araujo et al. (2015) Araujo R. A. N., Winter O. C., Prado A. F. B. A., 2015, mnras, 449, 4404
  * Araujo et al. (2017) Araujo R. A. N., Moraes R. V., Prado A. F. B. A., Winter O. C., 2017, mnras, 472, 3999
  * Baer et al. (2011) Baer J., Chesley S. R., Matson R. D., 2011, aj, 141, 143
  * Barbosa Torres dos Santos et al. (2017) Barbosa Torres dos Santos L., Bertachini de Almeida Prado A. F., Merguizo Sanchez D., 2017, Astrophysics and Space Science, 362, 61
  * Barbosa Torres dos Santos et al. (2020) Barbosa Torres dos Santos L., Marchi L., Sousa-Silva P. A., Sanchez D. M., Aljbaae S., Prado A. F. B. A., 2020, rmxaa, 56, 269
  * Bartczak & Breiter (2003) Bartczak P., Breiter S., 2003, Celestial Mechanics and Dynamical Astronomy, 86, 131
  * Benacchio (2007) Benacchio T., 2007, Master’s thesis, Universit‘a degli studi di Padova
  * Beutler (2005) Beutler G., 2005, Methods of celestial mechanics. Vol. II: Application to planetary system geodynamics and satellite geodesy
  * Blesa (2006) Blesa F., 2006, Monografías del Seminario Matemático García de Galdeano, 33, 67
  * Broucke (1968) Broucke R. A., 1968, Periodic orbits in the restricted three-body problem with earth-moon masses
  * Carry (2012) Carry B., 2012, planss, 73, 98
  * Chanut et al. (2014) Chanut T. G. G., Winter O. C., Tsuchida M., 2014, mnras, 438, 2672
  * Chanut et al. (2015a) Chanut T. G. G., Aljbaae S., Carruba V., 2015a, mnras, 450, 3742
  * Chanut et al. (2015b) Chanut T. G. G., Winter O. C., Amarante A., Araújo N. C. S., 2015b, mnras, 452, 1316
  * Chermnykh (1987) Chermnykh S. V., 1987, Vest. Leningrad Univ, 2, 10 pp.
  * Elipe & Lara (2003) Elipe A., Lara M. A., 2003, Journal of the Astronautical Sciences, 51, 391
  * Gabern et al. (2006) Gabern F., Koon W. S., Marsden J. E., Scheeres D. J., 2006, SIAM Journal on Applied Dynamical Systems, 5, 252
  * Jain & Sinha (2014) Jain R., Sinha D., 2014, Astrophysics and Space Science, 351, 87
  * Jiang & Baoyin (2018) Jiang Y., Baoyin H., 2018, Advances in Space Research, 62, 3199
  * Kirpichnikov & Kokoriev (1988) Kirpichnikov S. N., Kokoriev A. A., 1988, Vest. Leningrad Univ, 3, 73
  * Kokoriev & Kirpichnikov (1988) Kokoriev A. A., Kirpichnikov S. N., 1988, Vest. Leningrad Univ, 1, 75
  * Lan et al. (2017) Lan L., Yang H., Baoyin H., Li J., 2017, Astrophysics and Space Science, 362, 169
  * Liu et al. (2011a) Liu X., Baoyin H., Ma X., 2011a, Astrophysics and Space Science, 333, 409
  * Liu et al. (2011b) Liu X., Baoyin H., Ma X., 2011b, Astrophysics and Space Science, 334, 357
  * Montenbruck et al. (2002) Montenbruck O., Gill E., Lutze F., 2002, Applied Mechanics Reviews, 55, B27
  * Murray & Dermott (2000) Murray C. D., Dermott S. F., 2000, Solar System Dynamics
  * Ni et al. (2014) Ni Y., Baoyin H., Junfeng L., 2014, Proceedings of the International Astronautical Congress, IAC, 7, 4610
  * Pinheiro Lopes Masago et al. (2016) Pinheiro Lopes Masago B. Y., Prado A. F. B. A., Marins A. P., Chiaradia Martins Gomes V., 2016, Advances in Space Research, 57, 962
  * Riaguas et al. (1999) Riaguas A., Elipe A., Lara M., 1999, in Henrard J., Ferraz-Mello S., eds, Impact of Modern Dynamics in Astronomy. p. 169
  * Riaguas et al. (2001) Riaguas A., Elipe A., López-Moratalla T., 2001, Celestial Mechanics and Dynamical Astronomy, 81, 235
  * Rozitis & Green (2014) Rozitis B., Green S. F., 2014, aap, 568, A43
  * Sanchez & Prado (2019) Sanchez D. M., Prado A. F. B. A., 2019, Journal of Spacecraft and Rockets, 56, 1775
  * Scheeres et al. (2000) Scheeres D. J., Williams B. G., Miller J. K., 2000, Journal of Guidance Control Dynamics, 23, 466
  * Venditti et al. (2013) Venditti F. C. F., Rocco E. M., Prado A. F. B. A., 2013, in Journal of Physics Conference Series. p. 012008, doi:10.1088/1742-6596/465/1/012008
  * Wang et al. (2017) Wang W., Yang H., Zhang W., Ma G., 2017, Astrophysics and Space Science, 362, 229
  * Werner (1994) Werner R. A., 1994, Celestial Mechanics and Dynamical Astronomy, 59, 253
  * Werner & Scheeres (1996) Werner R. A., Scheeres D. J., 1996, Celestial Mechanics and Dynamical Astronomy, 65, 313
  * Worthington (2012) Worthington J., 2012, Master’s thesis, University of Sydney
  * Yang et al. (2015) Yang H.-W., Zeng X.-Y., Baoyin H., 2015, Research in Astronomy and Astrophysics, 15, 1571
  * Yang et al. (2017) Yang H., Baoyin H., Bai X., Li J., 2017, Astrophysics and Space Science, 362, 27
  * Yang et al. (2018) Yang H.-W., Li S., Xu C., 2018, Research in Astronomy and Astrophysics, 18, 084
  * Zeng et al. (2015) Zeng X., Jiang F., Li J., Baoyin H., 2015, Astrophysics and Space Science, 356, 29
  * Zeng et al. (2016a) Zeng X., Gong S., Li J., Alfriend K. T., 2016a, Journal of Guidance Control Dynamics, 39, 1223
  * Zeng et al. (2016b) Zeng X., Baoyin H., Li J., 2016b, apss, 361, 14
  * Zeng et al. (2016c) Zeng X., Baoyin H., Li J., 2016c, apss, 361, 15
  * Zeng et al. (2017) Zeng X.-Y., Liu X.-D., Li J.-F., 2017, Research in Astronomy and Astrophysics, 17, 2
  * Zeng et al. (2018) Zeng X., Zhang Y., Yu Y., Liu X., 2018, aj, 155, 85