# Semi-analytical near-Earth objects propagation:
the orbit history of (35107) 1991 VH and (175706) 1996 FG3

Oscar Fuentes-Muñoz Research Assistant
Ann and H. J. Smead Department of Aerospace Engineering Sciences, University
of Colorado Boulder, 3775 Discovery Dr.,
Boulder, CO 80303, USA  Alex J. Meyer Research Assistant
Ann and H. J. Smead Department of Aerospace Engineering Sciences, University
of Colorado Boulder, 3775 Discovery Dr.,
Boulder, CO 80303, USA  Daniel J. Scheeres University of Colorado
Distinguished Professor
Ann and H. J. Smead Department of Aerospace Engineering Sciences, University
of Colorado Boulder, 3775 Discovery Dr.,
Boulder, CO 80303, USA

(Accepted June 30, 2022)

###### Abstract

The propagation of small bodies in the Solar system is driven by the
combination of planetary encounters that cause abrupt changes in their orbits
and secular long-term perturbations. We propose a propagation strategy that
combines both of these effects into a single framework for long-term, rapid
propagation of small bodies in the inner Solar System. The analytical secular
perturbation of Jupiter is interrupted to numerically solve planetary
encounters, which last a small fraction of the simulation time. The proposed
propagation method is compared to numerical integrations in the Solar system,
effectively capturing properties of the numerical solutions in a fraction of
the computational time. We study the orbital history of the Janus mission
targets: (35107) 1991 VH and (175706) 1996 FG3, obtaining a stochastic
representation of their long-term dynamics and frequencies of very close
encounters. Over the last million years the probability of a strongly
perturbing flyby is found to be small.

Asteroids (72), Near-Earth objects (1092), Astrodynamics (76), Asteroid
dynamics (2210), Close encounters (255), Opik theory (1162), Flyby missions
(545)

##  1 Introduction

As remnants of primordial planetary formation, Near Earth Objects (NEOs) are
relevant targets for scientific exploration. The number of discovered NEOs is
expected to continue increasing with future surveys, offering new
opportunities to the scientific community (Jones et al., 2018). The interest
and availability of specific objects is assessed with long term predictions of
their orbits. Trajectories of Near Earth Objects are dominated by close
encounters with the inner Solar System planets and secular perturbations
(Michel et al., 1996). Planetary encounters cause a high dependence on the
initial conditions as the flybys cause neighboring trajectories to diverge
(Tancredi, 1998). This divergence of the dynamics also causes the uncertainty
to grow very rapidly. Thus, an analytical model that captures the main
dynamical effects can avoid the computational cost of using high-fidelity
models while capturing the overall statistical evolution of the orbits
accurately. The semi-analytical propagation of asteroids allows the rapid
propagation of NEO orbits, in the interest of the analysis of large databases
of NEOs.

The long-term study of asteroid orbits has been achieved in the past using a
wide variety of analytical, semi-analytical and numerical methods. Analytical
methods are based on the study of the gravity potential to obtain secular and
resonant perturbations (Milani & Knežević, 1990). Semi-analytical methods are
used to map orbital elements to the locations of linear secular resonances,
which are resonances involving one planetary and one asteroid frequency
(Michel et al., 1997; Michel & Froeschlé, 1997). Both types of solutions
represent the dynamics of asteroids in the absence of planetary encounters by
averaging the perturbing potential.

On the other hand, previous studies focus on the accumulation of planetary
encounters in contrast to numerical integration (Dones et al., 1999). The
effect of close encounters on the orbit of asteroids can be computed using
analytical (Öpik, 1976), semi-analytical or numerical methods. Semi-analytical
solutions (Alessi & Sánchez, 2015) allow the computation of flybys treating
the planet as a perturbing force in the Lagrange Planetary Equations. Specific
numerical integrators are convenient to propagate orbits of asteroids in the
long-term, in which symplecticity is desired along with the capacity to
accurately solve close encounters (Wisdom & Holman, 1991; Chambers, 1999).
Under multiple resonances asteroids start to encounter planets while their
eccentricity increases. This increase often causes the asteroids to eventually
collide with the Sun, planets or to be ejected from the Solar System on a
hyperbolic orbit (Farinella et al., 1994; Gladman et al., 1997; Milani et al.,
1989; Dones et al., 1999; Michel et al., 2005).

In this paper we aim to provide a simulation framework for the propagation of
particles in the Solar System. Our approach consists in the analytical
propagation of the particle until a close encounter is found. The propagation
is stopped when the trajectory is close to a planet, then the close encounter
is evaluated numerically. The evaluation of the encounters is based on a
quadrature of the Lagrange Planetary Equations (LPE) around the closest
approach date. After the encounter the analytical propagation of the orbit is
resumed. The propagation under secular perturbations provides a realistic
prediction of when the next encounter can occur as the orbit of the asteroid
drifts between different regions of the inner Solar System. This approach
reduces substantially the computational time of solutions obtained entirely by
numerical integration while providing deeper insight into the dynamics.

The use of the analytical secular model allows the prediction of long-term
properties of the asteroid dynamics. Eccentricities, inclinations and angles
of asteroid and planets drift secularly. Thus, we can propagate the minimum
orbit intersection distance (MOID). The MOID constrains the minimum closest
approach distance between the asteroid and the planets and defines if
asteroids are potentially hazardous (PHAs). The long-term dynamics of the
orbits of NEOs and the MOID are studied by sampling a large number of virtual
asteroids from their uncertainty distributions. We use the semi-analytical
propagation of these asteroids to show the stochastic nature of the orbital
evolution of NEOs.

The semi-analytical propagation allows us to track the encounters experienced
by asteroids in the inner Solar System, which can perturb the physical
properties of asteroids. The orbits of binary asteroids can be disrupted by a
very close encounter (Meyer & Scheeres, 2021). In this paper we study the
orbital history of the targets of the exploration mission Janus (Scheeres et
al., 2020): the two binaries (35107) 1991 VH and (175706) 1996 FG3. The
stochastic long-term dynamics in the last million years are modelled by
sampling a large number of particles from their current orbit uncertainties.
We model the evolution of these statistical distributions by a random walk in
semi-major axis, eccentricity and inclination and a uniform distribution in
longitude of perihelion. Then, we compute the probability that (35107) 1991 VH
and (175706) 1996 FG3 could have been potentially disrupted by a close
encounter in this period of a million years.

This section presents the scope of the paper and is followed by the background
of this work in section 2. Next, section 3 describes the methodology including
a detailed study of flybys evaluation and the derivation of an analytical
N-body secular problem solution. Section 4 shows examples of the long-term
propagation of asteroids and how the long-term dynamics can be characterized
stochastically. Section 5 discusses the limitations of the semi-analytical
propagation tool. Section 6 studies the orbital history of the Janus targets
and the frequency of close encounters. Last, section 7 concludes by evaluating
the aspects in which this methodology proves beneficial, questions that
remained unanswered, and future work.

##  2 Background

The long-term dynamics of NEOs are governed by their gravitational
interactions with the other bodies of the Solar System. The effects of the
most massive and external planets have timescales of millennia. However,
planetary close encounters can abruptly change an orbit over a timescale of
days. The accumulation of such planetary encounters cause the orbits of NEOs
to be chaotic (Tancredi, 1998). This section describes this phenomenon in more
detail. The evaluation of close encounters is necessary for the propagation of
NEOs, hence the variety of possible flybys is demonstrated later for the
validation of the method.

Many asteroids experience long periods of time without flybys. The dominant
dynamics in those periods of time are the secular perturbations from massive
planets in the Solar System. Likewise, the orbits of the planets evolve
secularly over similar timescales. The Laplace-Lagrange secular theory
qualitatively describes the evolution of the elements of the planets at any
distant time in the future or past. As for the asteroid, the secular solution
from external perturbers represents the orbital dynamics of asteroids between
encounters.

The presence of repeated encounters is one of the main characteristics of the
long-term propagation of asteroids in the inner Solar System. Repeated close
encounters cause a random walk in the elements of the asteroids. Very close
encounters occur less frequently but change substantially the orbits of NEOs,
modifying predictions on the long-term evolution of their orbits. Thus, we
propose an informed analytical propagation of the orbits while characterizing
planetary close encounters. The proposed methodology is born from the
combination of these two dynamical regimes: the long-term effects of secular
dynamics and the frequent changes in elements experienced in planetary
encounters. Considering the secular drift of the asteroid we model the
seasonal variation of the possible encounters with planets.

###  2.1 Chaotic dynamics in the inner Solar System

An accurate description of the evolution of orbits of near-Earth asteroids
beyond a few centuries is challenging. This is because the succession of
planetary encounters disperse neighboring trajectories to become chaotic
(Tancredi, 1998). Small deviations in the orbital period change the timing of
the flybys, spreading the uncertainty along the Line of Variation (Milani et
al., 2005). After successive flybys the resulting imaginary stream of
particles is spread in highly non-linear distributions. For this reason the
study of long-term dynamics is often left to a statistical analysis requiring
a large number of particles and computational efforts. In this context we
propose the use of this semi-analytical tool to obtain long-term simulations
in short computational times.

Figure 1: Chaotic dynamics of (35107) 1991 VH as obtained from numerical
integration. Each axis represents the variation from the initial value of
elements pairs: (left) semi-major axis-eccentricity, (center) semi-major axis-
inclination, (right) argument of the node-argument of perihelion in degrees.
The orbital evolution is shown at four instants of time: initial (first row),
after 500 years (second row), after 5000 years (third row) and after 1 million
years (bottom).

We exemplify the sensitivity to initial conditions in a numerical integration
of asteroid (35107) 1991 VH, which is one of the two targets of mission Janus
(Scheeres et al., 2020), a NASA SIMPLEx mission. Figure 1 shows 1440 particles
generated from the uncertainty in the orbit solution of (35107) 1991 VH, which
is included in appendix A. These particles are propagated in the N-body
integrator IAS15 (Rein & Spiegel, 2014) including the Solar System planets
from Venus to Neptune. The particles are propagated for a million years,
although in this section we study in more detail the distributions after
shorter periods of time.

After 500 years the initial normal distribution already becomes a stream of
particles. While the variation in the elements from the nominal is similar for
all the particles, there is a dispersion orders of magnitude smaller that
represents the stream of particles. After 5000 years, the distribution becomes
completely different: the presence of planetary encounters disperses the
particles around the initial orbit. The variation in eccentricities and
inclinations has a secular component. However, the variation on the argument
of the node and argument of perihelion is dominantly secular after a few
millennia. After a million years, the particles are spread along a large
region of near-Earth space. In argument of perihelion and ascending node we
observe that the distribution becomes almost uniform in the whole 2D angular
space.

The secular drift in the arguments defines the possibility of encounters over
time. For this reason, it is important to characterize this drift and the
secular cycles under the perturbation of the large bodies of the Solar System.
When encounters are possible with the inner Solar System planets, these need
to be accounted as perturbers of the orbit evolution.

The stochastic nature of the long-term dynamics of NEOs under close encounters
implies that the precise determination of their position after hundreds of
thousands of years is unachievable. However, we can still collect statistics
that give us insight on their orbital history. Another implication is that the
inclusion of higher order dynamics is shadowed under the stochastic dispersion
caused by the main gravitational perturbations. For example, the magnitude of
the Yarkovsky effect is typically $~{}10^{-4}$ au/Myr (Vokrouhlický et al.,
2000; Nesvorný & Bottke, 2004), which is still two orders of magnitude smaller
than a typical dispersion after 10,000yr under repeated close encounters, as
observed in the example of Figure 1. In section 6.1 we show that (35107) 1991
VH is not under a particularly high frequency of close encounters compared to
other NEOs.

Similarly, relativistic effects can have a non-negligible effect in the
secular rates of the argument of perihelion. These are usually measured in
arcseconds per year or century, and typical values are 1-2 orders of magnitude
smaller than the typical secular periods of the order of 100,000 years
(Benitez & Gallardo, 2008). Even if the secular rate has an error, the
presence of encounters already causes the distributions to become uniform in
argument of the node and perihelion after a few secular periods.

###  2.2 NEO close encounters in the inner Solar System

Flybys can occur with multiple planets over short periods of time. Even if the
encounters are with the same planet, the closest approach distance and
relative velocities can change depending on the timing of the flyby. The
geometry of the flyby is constrained by the heliocentric elements of the
asteroid. If shallow encounters are considered, the position in the asteroid
orbit in which the planet is encountered can significantly change the relative
velocity. These variations are not well captured by analytical theories, but
the proposed propagation tool aims to accurately model these variations. These
are different regimes of flybys in which the evaluation tool needs to be
accurate.

Figure 2: Relative velocity at closest approach of flybys generated from the
propagation of NEOs with semi-major axis smaller than 2 au for 50 years. The
symbol indicates the planet that the asteroid is encountering and the
horizontal axis shows respectively closest approach distance (au), inclination
(deg) and eccentricity.

In order to broadly show the diversity in flybys that different NEOs
experience, we generate a list of flybys that will be used to validate the
evaluation of close encounters. From the database of NEOs we select the ones
with semi-major axis smaller than 2 au (JPL Solar System Dynamics & Center for
NEO Studies (2021), CNEOS). Then, we propagate their positions using the
secular model for 50 years. For such a brief period of time the change in the
elements is insignificant for our purposes. Figure 2 shows more than 30,000
flybys generated with the described method. Shallow encounters are much more
frequent than the very close encounters that cause large variations in orbit
elements. Thus, we want to consider them even if their individual contribution
is not as significant.

The range of possible relative velocities in figure 2 depends on the planet in
question, with increasing maximum relative velocity for the planet closest to
the Sun. The relative velocity is defined by the heliocentric orbit of the
asteroid, with an increasing range of possible values depending on the
inclination and eccentricity of the orbit. Overall, after millions of years
asteroids experience a variety of encounters that can be computed with
different methods. With this purpose the list of generated flybys is used to
decide the method to compute the post-encounter elements of flybys. In section
3 we compute the error of different close encounter evaluation methods
referenced to numerical integration of the trajectories.

##  3 Methodology

This semi-analytical propagation tool consists in the following process.
First, the orbit of the asteroid is propagated by an analytical secular
solution. This perturbed motion is interrupted when an encounter is found with
a nearby planet. Then, the trajectory during the planetary encounter is
modeled using a numerical method. Next, the secular propagation is continued
until the subsequent encounter.

The simplest way to find encounters is to track the distance between the
asteroid and planets at all times. While the regions in which encounters are
possible are determined by the geometry of the asteroid around the central
body, searching at all times is the most generic approach. The state of the
planets is obtained from the secular solution of the 8 main planets
interacting with each other. The state of the small body is corrected given
the secular dynamics model. Once we determine the initial conditions of the
encounter, the change in orbit elements is computed through the proposed
numerical procedure.

There are many methods to compute planetary encounters available in the
literature. Analytical solutions for Keplerian elements before and after close
encounters in Öpik’s Theory (Öpik, 1976) were extended for multiple
applications by Valsecchi et al. (2003, 2015). However, these analytical
expressions are constrained to encounters that are very close and small bodies
that are not co-moving with the planet. Asteroid and planet are co-moving when
they have a small inclination and at least one of the node crossings close to
the planet orbit.

We name shallow encounters those with large close approach distance but non-
negligible effects. Shallow encounters are very frequent and influence the
long-term evolution of small bodies in the Solar System. In order to account
for shallow encounters, semi-analytical methodologies can be used to map
before and after encounter conditions (Alessi & Sánchez, 2015). These methods
are based on the quadrature of Lagrange Planetary Equations around the
encounter. In this work we derive a quadrature of Lagrange Planetary Equations
in Delaunay elements which is solved using a numerical integration scheme. In
the case of extremely close or slow encounters we solve the encounters
numerically.

Once the solution of most encounters is obtained satisfactorily we focus
efforts in the computation of the perturbed motion of the asteroid in absence
of encounters. We evaluate the solution of N-bodies interacting secularly to
generate the orbits of the planets. Then, we obtain the perturbed motion of
the asteroid including only the planets relevant to its secular influence.
Taking into account the influence of only Jupiter is a valid generic approach
to estimate the secular dynamics of NEOs (Vokrouhlický et al., 2012; Pokorný &
Vokrouhlický, 2013; Fuentes-Munoz & Scheeres, 2020). In this work we use the
Laplace-Lagrange secular model. The secular rates as obtained by the
analytical theory are compared to numerical integration to validate the range
of validity of the solution. This defines a range of applicability of the
tool, as we discuss later.

In this section we compare the individual pieces of the semi-analytical
propagation tool to numerical methods. Last, we compare the combined semi-
analytical propagation tool with trajectories obtained through numerical
integration and evaluate the computational efficiency of the method.

###  3.1 Analytical secular dynamics of multibody systems

The dynamical landscape of the Solar System is complex with gravitational
interactions between all planets. This landscape leads to resonances and
secular motion in asteroids in the system. Well inside the inner Solar System,
the dynamics are dominantly secular. The secular solution of a planetary
system formed by N-planets can be obtained analytically to first order in
inclinations and eccentricities and in the absence of resonances. This section
derives an implementation of the solution following the procedure in Chapter 7
of Murray & Dermott (2000). The perturbing potential is written for the $N$
bodies considered. Then Lagrange Planetary Equations are used to compute the
equations of motion of the elements of each particle, leading to a system of
differential equations solved together.

The secular model is obtained as follows: (1) The perturbing potential is
split in a direct part and an indirect part based on the dependency on fast
angles, (2) then the perturbing potential is expanded in Keplerian Elements.
(3) The important terms of the expansion are selected based on the averaging
principle. (4) The terms are rewritten in semi-equinoctial elements to ease
the solution of the global system of equations. (5) Take the necessary
partials to solve the set of Lagrange Planetary Equations. The perturbing
potential experienced by a mass $j$ by a second mass $k$ is:

$R_{jk}=\frac{Gm_{k}}{a_{k}}\left({R_{jk}}_{D}+{R_{jk}}_{I}\right)$ (1)

Where $a_{k}$ is the semi-major axis of body the external body. The perturbing
potential is separated in the direct ${R_{jk}}_{D}$ and indirect
${R_{jk}}_{I}$ parts:

$\begin{aligned}
{R_{jk}}_{D}=\frac{a_{k}}{\left|{r}_{j}-{r}_{k}\right|}\end{aligned}\qquad\qquad\begin{aligned}
{R_{jk}}_{I}=-\frac{a_{k}^{2}}{a_{j}}\frac{{r}_{j}\cdot{r}_{k}}{\left|{r}_{k}\right|^{3}}\end{aligned}$
(2)

The separation is convenient to expand in the ratio of semi-major axes
$\alpha_{jk}$ as well as sines and cosines of
$\left\\{\varpi_{j},\Omega_{j},\lambda_{j},\varpi_{k},\Omega_{k},\lambda_{k},\right\\}$.
The ratio of semi-major axes is $\alpha_{jk}=a_{j}/a_{k}$ if the perturber is
external, or $\alpha_{jk}=a_{j}/a_{k}$ if the perturber is internal. All the
terms that depend on the longitudes $\left\\{\lambda,\lambda_{j}\right\\}$ are
of short-period, so it can be argued that they do not contribute to the
averaged potential $R_{j}$. The secular potential lowest order in
eccentricities and inclinations is:

$R_{j}=R_{0,j}+R_{1,j}=\sum^{N}_{k=1,k\neq
j}Gm_{k}\frac{1}{2a_{k}}b^{(0)}_{1/2}\left(\alpha_{jk}\right)+R_{1,j}$ (3)

$R_{1,j}=n_{j}a^{2}_{j}\left[\frac{1}{2}A_{jj}e^{2}_{j}+\frac{1}{2}B_{jj}I^{2}_{j}+\sum^{N}_{\begin{subarray}{c}k=1\\\
k\neq
j\end{subarray}}A_{jk}e_{j}e_{k}\cos{\left(\varpi_{j}-\varpi_{k}\right)}+B_{jk}I_{j}I_{k}\cos{\left(\Omega_{j}-\Omega_{k}\right)}\right]$
(4)

Where $\bar{\alpha}_{jk}=a_{j}/a_{k}$ if the perturber is external, or
$\bar{\alpha}_{jk}=1$ if the perturber is internal. The coefficients
$A_{jj},A_{jk},B_{jj},B_{jk}$ are:

$A_{jk}=-n_{j}\frac{1}{4}\frac{m_{k}}{m_{c}+m_{j}}\alpha_{jk}\bar{\alpha}_{jk}b^{(2)}_{3/2}\left(\alpha_{jk}\right)$
(5)
$B_{jk}=+n_{j}\frac{1}{4}\frac{m_{k}}{m_{c}+m_{j}}\alpha_{jk}\bar{\alpha}_{jk}b^{(1)}_{3/2}\left(\alpha_{jk}\right)$
(6) $A_{jj}=+n_{j}\frac{1}{4}\sum^{N}_{k=1,k\neq
j}\frac{m_{k}}{m_{c}+m_{j}}\alpha_{jk}\bar{\alpha}_{jk}b^{(1)}_{3/2}\left(\alpha_{jk}\right)=\sum^{N}_{k=1,k\neq
j}B_{jk}$ (7) $B_{jj}=-n_{j}\frac{1}{4}\sum^{N}_{k=1,k\neq
j}\frac{m_{k}}{m_{c}+m_{j}}\alpha_{jk}\bar{\alpha}_{jk}b^{(1)}_{3/2}\left(\alpha_{jk}\right)=-\sum^{N}_{k=1,k\neq
j}B_{jk}$ (8)

where the coefficients $b^{(k)}_{s}$ are Laplace Coefficients. More details on
their computation can be found in Appendix B. The coefficients
$A_{jj},A_{jk},B_{jj},B_{jk}$ form the matrices ${A}$ and ${B}$. We can
rewrite the potential in semi-equinoctial elements,

$\begin{aligned} h_{j}&=e_{j}\sin\varpi_{j}\\\ k_{j}&=e_{j}\cos\varpi_{j}\\\
\end{aligned}\qquad\qquad\begin{aligned} p_{j}&=I_{j}\sin\Omega_{j}\\\
q_{j}&=I_{j}\sin\Omega_{j}\\\ \end{aligned}$ (9)

the potential becomes:

$R_{1,j}=n_{j}a^{2}_{j}\left[\frac{1}{2}A_{jj}(h_{j}^{2}+k_{j}^{2})+\frac{1}{2}B_{jj}(p_{j}^{2}+q_{j}^{2})+\sum^{N}_{k=1,k\neq
j}A_{jk}(h_{j}h_{k}+k_{j}k_{k})+B_{jk}(p_{j}p_{k}+q_{j}q_{k})\right]$ (10)

Our complete set of states includes the mean anomaly at epoch $\sigma_{j}$ and
$L_{j}=\sqrt{GMa_{j}}$. The equations of motion become:

$\begin{aligned} \dot{p}_{j}&=\frac{1}{n_{j}a^{2}_{j}}\frac{\partial
R_{j}}{\partial q_{j}}\\\ \dot{q}_{j}&=-\frac{1}{n_{j}a^{2}_{j}}\frac{\partial
R_{j}}{\partial p_{j}}\\\ \end{aligned}\qquad\qquad\begin{aligned}
\dot{h}_{j}&=\frac{1}{n_{j}a^{2}_{j}}\frac{\partial R_{j}}{\partial k_{j}}\\\
\dot{k}_{j}&=-\frac{1}{n_{j}a^{2}_{j}}\frac{\partial R_{j}}{\partial h_{j}}\\\
\end{aligned}\qquad\qquad\begin{aligned} \dot{L}_{j}&=\frac{\partial
R_{j}}{\partial\sigma_{j}}\\\ \dot{\sigma}_{j}&=-\frac{\partial
R_{j}}{\partial L_{j}}\\\ \end{aligned}$ (11)

The solution of $h_{j},k_{j},p_{j},q_{j}$ only depends on $R_{1,j}$. For this
reason the perturbing potential is often only expressed with those components.
However, if we want the solution of the mean anomaly at epoch $\sigma_{j}$ it
is necessary to take into account $R_{0,j}$. In the process of averaging the
terms that would effect the semi-major axis are removed, meaning that under
this assumption that element remains constant. The solution of
$h_{j},k_{j},p_{j},q_{j}$ is:

$\begin{aligned}
h_{j}(t)&=\sum^{N}_{i=1}e_{ji}\sin\left(g_{i}t+\beta_{i}\right)\\\
k_{j}(t)&=\sum^{N}_{i=1}e_{ji}\cos\left(g_{i}t+\beta_{i}\right)\\\
\end{aligned}\qquad\qquad\begin{aligned}
p_{j}(t)&=\sum^{N}_{i=1}I_{ji}\sin\left(f_{i}t+\gamma_{i}\right)\\\
q_{j}(t)&=\sum^{N}_{i=1}I_{ji}\cos\left(f_{i}t+\gamma_{i}\right)\\\
\end{aligned}$ (12)

where two sets of eigenvalue problems are solved for
$e_{ji},I_{ji},f_{i},g_{i}$. The frequencies $g_{i}$ are the eigenvalues of
${A}$, and the frequencies $f_{i}$ are the eigenvalues of ${B}$. $e_{ji}$ and
$I_{ji}$ are related to the eigenvectors of ${A}$ and ${B}$, but need to be
solved with $\beta_{i},\gamma_{i}$ given a set of initial conditions. In order
to solve for $e_{ji},I_{ji},\beta_{i},\gamma_{i}$ we proceed as follows. From
the matrices of normalized eigenvectors $\bar{e}_{ji},\bar{I}_{ji}$ and the
initial conditions ${h},{k},{p},{q}$ we form:

$\begin{aligned} {h}&=\bar{e}_{ji}\left[S_{i}\sin\beta_{i}\right]\\\
{k}&=\bar{e}_{ji}\left[S_{i}\cos\beta_{i}\right]\\\
\end{aligned}\qquad\qquad\begin{aligned}
{p}&=\bar{I}_{ji}\left[T_{i}\sin\gamma_{i}\right]\\\
{q}&=\bar{I}_{ji}\left[T_{i}\cos\gamma_{i}\right]\\\ \end{aligned}$ (13)

These are four linear systems of equations, where $S_{i},T_{i}$ are the
scaling factors of each eigenvector. Solving for the combined factors
$\left[S_{i}\sin\beta_{i}\right],\left[S_{i}\cos\beta_{i}\right],\left[T_{i}\sin\gamma_{i}\right]$
and $\left[T_{i}\cos\gamma_{i}\right]$ we can reconstruct the vectors
$e_{ji},I_{ji}$ and the phase angles $\beta_{i},\gamma_{i}$.

Table 1: Initial conditions of the Solar System propagation in Figure 3.

Planet |  $a$ (au) | $e$ |  $i$ (deg) |  $\Omega$ (deg) |  $\omega$ (deg) |  $M_{0}$ (deg)
---|---|---|---|---|---|---
Mercury | 0.39703 | 0.21337 | 6.936 | 48.264 | 31.991 | 52.745
Venus | 0.73096 | 0.012687 | 3.378 | 76.799 | 45.020 | 16.566
Earth | 1.0030 | 0.018402 | 0.001 | 154.979 | 296.322 | 8.654
Mars | 1.5177 | 0.093083 | 1.852 | 49.461 | 288.507 | 322.879
Jupiter | 5.1904 | 0.047388 | 1.305 | 100.514 | 273.897 | 353.761
Saturn | 9.5499 | 0.05412 | 2.487 | 113.612 | 339.598 | 91.261
Uranus | 19.207 | 0.04628 | 0.772 | 73.997 | 96.864 | 189.506
Neptune | 30.109 | 0.0091006 | 1.770 | 131.780 | 265.440 | 291.693

Note. — From ephemeris DE431 at Epoch: JD${}_{0}=2455562.5$ (2011 January 1)
TDB

Figure 3: Semi-equinoctial elements of the inner Solar System planets obtained
using three models. The analytical secular model derived is shown in blue, the
integration of the 9BP of the main planets and the Sun in red and the
ephemeris file DE431 is shown in grey. The initial conditions are obtained
from the time average of the ephemeris file over the first orbital periods of
the planets.

Figure 3 shows the solution of eq. 12 for Mercury, Venus, Earth and Mars as
perturbed mutually and from the rest of planets of the Solar System. This
model is compared to two other models for 15,000 years into the past. The
first one is a numerical integration of the N-body problem taking into account
the main 8 planets of the Solar System and the Sun. Then, we also compare to
the planetary ephemerides DE431 (Folkner et al., 2014). While the complete
ephemerides models show the short period effects, the secular component is
modeled by the two simplified models.

The initial conditions of the 9BP integration and the secular theory are
obtained by averaging the full ephemeris model for two orbit periods. As a
result, the initial conditions visually appear to be off from the mean of the
full ephemeris solution, but they are the actual time average.

The short-term components have a significant effect in the evolution of
$h_{j},k_{j}$. In the case of $p_{j},q_{j}$, the secular component is the
dominant effect of the evolution. Because of the assumptions of small
eccentricities and inclinations, the predicted frequencies are not perfectly
accurate, as observed in the drift between the 9BP solution and the analytical
theory of figure 3, especially in $p_{j},q_{j}$.

A similar agreement in these elements is found for the gas giants. However,
only the inner Solar System planets are shown as they are the bodies that are
encountered by near-Earth objects. Thus, these are the planets for which we
want to guarantee an accurate model of their secular dynamics. The analytical
propagation of Mercury drifts the most from the full ephemeris solution,
although it is the least relevant inner planet. Close encounters with Mercury
are unfrequent and have a small effect, as Mercury is the least massive planet
and it is encountered with very high relative velocity.

As a result of the averaging of the perturbing potential, the semi-major axis
of the bodies remains constant. The complete set of secular solutions includes
the mean anomaly at epoch $\sigma_{j}$. Short term applications benefit from
the improved characterization of the position of the bodies in their orbits.
Solving for $\sigma_{j}$ is straightforward if we ignore the contribution of
$R_{1,j}$, which has a small effect compared to $R_{0,j}$. The equation of
motion for $\sigma_{j}$ becomes:

$\dot{\sigma}_{j}=-\frac{\partial R_{j}}{\partial
L_{j}}=-\frac{2}{n_{j}a_{j}}\frac{\partial R_{0,j}}{\partial a_{j}}$ (14)

and the solution depends on whether the perturber is external or internal:

$\dot{\sigma}_{j}=-\sum^{N}_{k=1,k\neq
j}\frac{Gm_{k}}{n_{j}a^{2}_{j}}\bar{c}_{jk}Db^{(0)}_{1/2}$ (15)

where $\bar{c}_{jk}=a_{j}/a^{2}_{k}$ in the case of an external perturber and
$\bar{c}_{jk}=-1/a_{j}$ if the perturber is internal. The solution of the
equation is simply a constant drift given by the rate $\dot{\sigma}$. This
element completes the set of elements of the secular model.

###  3.2 Analytical secular dynamics of near-Earth asteroids

The secular dynamics of asteroids can be modelled as a particular case of the
secular dynamics of multibody systems described above. In the present work we
apply this solution to the evolution of the asteroid under the external
perturbation of Jupiter. The solutions of equation 13 simplify in the case of
a system of 2 bodies with a massless internal body. We follow the same process
to obtain the solution. Matrices $A_{jk}$ and $B_{jk}$ simplify to:

$\begin{aligned} A_{jk}=\left(\begin{matrix}B_{12}&A_{12}\\\
0&0\end{matrix}\right)\end{aligned}\qquad\qquad\begin{aligned}
B_{jk}=\left(\begin{matrix}-B_{12}&B_{12}\\\
0&0\end{matrix}\right)\end{aligned}$ (16)

where the subindexes $1,2$ correspond respectively to the massless particle
and the external perturber. The coefficients of the matrices are found as in
equations 5-8 above. The solution to the eigenvalue problem yields the secular
frequencies of the secular propagation $g_{1}=B_{12}$, $g_{2}=0$,
$f_{1}=-B_{12}$, $f_{2}=0$. As expected, the elements of the perturber
$h_{2},k_{2},p_{2},k_{2}$ remain constant. The eigenvectors are the columns of
the matrices:

$\begin{aligned} \bar{e}_{jk}=\left(\begin{matrix}1&\kappa\\\
0&1\end{matrix}\right)\end{aligned}\qquad\qquad\begin{aligned}
\bar{i}_{jk}\left(\begin{matrix}1&\frac{\sqrt{2}}{2}\\\
0&\frac{\sqrt{2}}{2}\end{matrix}\right)\end{aligned}$ (17)

where the constant $\kappa$ is found as the ratio between Laplace
coefficients:

$\kappa=\frac{A_{12}}{-B_{12}}=\frac{b^{(2)}_{3/2}}{b^{(1)}_{3/2}}$ (18)

Note that the vector ($e_{12},e_{22}$) is not normalized. This is not
necessary because in the process of obtaining the integration constants from
the initial conditions the scaling of the eigenvectors is found. The solution
of the elements of the massless particle becomes:

$\begin{aligned} h_{1}{(t)}&=S_{1}\sin{(g_{1}t+\beta_{1})}+\kappa h_{2}\\\
k_{1}{(t)}&=S_{1}\cos{(g_{1}t+\beta_{1})}+\kappa k_{2}\\\
\end{aligned}\qquad\qquad\begin{aligned}
p_{1}{(t)}&=T_{1}\sin{(f_{1}t+\gamma_{1})}+p_{2}\\\
q_{1}{(t)}&=T_{1}\cos{(f_{1}t+\gamma_{1})}+q_{2}\\\ \end{aligned}$ (19)

with constants of integration:

$\begin{aligned} S^{2}_{1}&=e^{2}_{1,0}+\kappa^{2}e^{2}_{2}-2\kappa
e_{1,0}e_{2}\cos{(\varpi_{1,0}-\varpi_{2})}\\\
T^{2}_{1}&=i^{2}_{1,0}+i^{2}_{2}-2i_{1,0}i_{2}\cos{(\Omega_{1,0}-\Omega_{2})}\\\
\end{aligned}\qquad\qquad\begin{aligned} \tan\beta_{1}&=\frac{h_{1,0}-\kappa
h_{2}}{k_{1,0}-\kappa k_{2}}\\\
\tan\gamma_{1}&=\frac{p_{1,0}-p_{2}}{q_{1,0}-q_{2}}\\\ \end{aligned}$ (20)

The time evolution of the Keplerian elements set can be obtained from the
relationships with the semi-equinoctial set in equation 9. The solutions of
$\varpi(t),\Omega(t)$ are the secular drift with frequencies $g_{1},f_{1}$
that are equal with opposite signs. The solutions of $e(t),i(t)$ are
oscillations with frequencies $g_{1},f_{1}$ as obtained from the development
of eccentricity $e_{1}(t)=\sqrt{h^{2}_{1}(t)+k^{2}_{1}(t)}$ and inclination
$i_{1}(t)=\sqrt{p^{2}_{1}(t)+q^{2}_{1}(t)}$. The maximum and minimum values of
eccentricity and inclination are:

$\begin{aligned} e^{2}_{1,min}&=S^{2}_{1}+\kappa^{2}e^{2}_{2}-2S_{1}\kappa
e_{2}\\\ e^{2}_{1,max}&=S^{2}_{1}+\kappa^{2}e^{2}_{2}+2S_{1}\kappa e_{2}\\\
\end{aligned}\qquad\qquad\begin{aligned}
i^{2}_{1,min}&=T^{2}_{1}+i^{2}_{2}-2T_{1}i_{2}\\\
i^{2}_{1,max}&=T^{2}_{1}+i^{2}_{2}-2T_{1}i_{2}\\\ \end{aligned}$ (21)

The secular model is computed for the fictitious asteroid of Case 1 of table 2
with the perturbation of Jupiter given by the elements of table 1. These cases
are used later to demonstrate the propagation tool. For a nominal eccentricity
of 0.15 the minimum eccentricity is 0.14946 and maximum is 0.17466. For a
nominal inclination of 10 degrees, the minimum inclination is 7.41823 degrees
and maximum is 10.02508 degrees. The characteristic period of the secular
motion $T_{sec}$ is 154,116 years.

This model assumes small eccentricities and inclinations. While these
conditions are usually not fulfilled, it is important to remark that
eccentricity and inclination are under frequent disturbance due to close
encounters. Most importantly, the secular drift in $\Omega,\omega$ controls
the evolution of the possible planetary encounters.

The assumptions on the heliocentric orbit of the asteroid for the analytical
secular perturbation solution are not always fulfilled among the NEO
population. In this section we show that the analytical theory represents the
dynamics of the perturbation by Jupiter. For this reason, we integrated the
orbits of 4462 NEOs with $e<0.7$ and $i<0.5$ rad for 50,000 years. Note that
in the solution of equation 19 if the terms of the external perturber are
small the solution tends to a linear drift of the angles $\Omega,\varpi$. In
addition, given the relationship between the frequencies $g$ and $f$, the
relationship between the arguments rates is $\dot{\omega}=-2\dot{\Omega}$.

Figure 4: Error in the secular rates of near-Earth objects as obtained by
linear regression of the propagation using Laplace-Lagrange secular theory and
compared to numerical integration of the three-body problem in percent.
Secular rates in ascending node (left) and argument of perihelion (center) as
function of the initial conditions semi-major axis and inclination. Secular
rates in the angles as obtained with the two methods and function of the
initial semi-major axis (right). The dashed lines indicate the region in which
we compute the average errors. This region includes the initial conditions
used throughout the paper, indicated as cross marks.

Figure 4 shows the secular rates computed from linear regression of the time
histories of $\Omega$,$\omega$. Note that from an initially larger list of
NEOs, a significant fraction (12%) was discarded because either eccentricity
or inclination were larger than 0.5. An additional 11% of the solutions were
discarded because the error in the regression was too large or during the
propagation close encounters with Jupiter moved the orbit of the NEO to a
completely different location than the initial conditions. While the linear
regression secular rates are not equivalent to the frequencies of the
analytical solution, they serve as a comparison between the two dynamical
models. The error is computed in percent relative to the rate measured from
the regression of the numerically integrated trajectories, given by:

$E(\dot{\omega})=100\left|\frac{\dot{\omega}_{sec}-\dot{\omega}_{3BP}}{\dot{\omega}_{3BP}}\right|$
(22)

As observed in Figure 4 the rates obtained with the two methods agree towards
the smaller end of semi-major axis. Since these are near-Earth objects, the
condition of being in the vicinity of Earth means that eccentricity increases
with semi-major axis. We can see that past 1.5-2 au the difference between the
two models is increased, as well as the secular rates values themselves also
increase. This difference is also appreciated in the rates as function of
semi-major axis, in which we show the agreement in the NEO region. Using the
current model we find secular rates for 60% of the population with an error
less than 30% in both $\dot{\omega}$ and $\dot{\varpi}$. If we limit the
application of the secular theory to semi-major axes between 0.8-1.4 au
(dashed region in Figure 4) we find that this agreement improves to 88%. It is
important to note that the examples chosen to demonstrate the semi-analytical
propagation tool fall within this region. Outside of this region we can verify
if the secular rates found are reliable by using numerical integration. This
test integration must be long enough to observe the secular rates, but still
orders of magnitude shorter than the time-scales that we can more efficiently
study using the semi-analytical propagation.

At larger semi-major axis the effect of mean-motion resonances becomes
important, and that Lidov-Kozai dynamics may better represent the dynamics for
large eccentricities and inclinations (Michel et al., 1996; Morbidelli et al.,
2009). The implementation of additional analytical long-term dynamics models
to model any generic asteroid is left as future work.

###  3.3 Finding the subsequent encounter

The analytical propagation of particles is interrupted when an encounter with
a planet is detected. In principle it is not necessary to track the distance
to planets at all times, since the regions in which encounters are possible
are determined by the geometry of the asteroid around the central body. If the
inclination relative to the planet is high, then the encounters are only
possible in the vicinity of the ascending and descending node. However, the
most generic approach is to track the distance between the asteroid and the
crossing planets at all times. Thus, the results in this work follow the
latter approach to find encounters within a closest approach distance of 0.1
au. When the two bodies are close, the unperturbed closest approach distance
is found using a bisection method where the function is the derivative of the
distance as obtained by finite differences. This process results in less
evaluations of the relative distance function based on the heliocentric
elements of the bodies.

The elements of the planets and the asteroid are propagated using the secular
solution at the date of start of the encounter, which is defined below. The
transition between models consists in the conversion between the sets of
elements, obtaining the necessary Keplerian elements in the process. These are
semi-equinoctial elements for the analytical perturbed propagation as in
equation 9 and Delaunay elements for the quadrature of the Lagrange planetary
equations.

###  3.4 Evaluation of planetary encounters

Close encounters are commonly solved using the analytical Öpik theory (Öpik,
1976). While this theory requires the least computational resources, its
accuracy is limited to specific circumstances. The quadrature of Lagrange
Planetary Equations can be used to solve close encounters (Alessi & Sánchez,
2015). In this work we derive a solution using this method for generic close
encounters using Delaunay elements. The two methods are compared to the
integration of the three body problem from the same date and during the same
period of time.

####  3.4.1 Öpik theory of close encounters

An analytical solution to the planetary close encounter problem was derived by
Öpik (Öpik, 1976). This solution was extended and studied in detail by
Valsecchi et al. (2003, 2015). The encounter solution uses a linearized
mapping from orbital elements to a planetocentric Cartesian frame, that is
later expressed in B-plane coordinates. Then, the encounter is assumed
instantaneous and the incoming asymptote and B-plane both rotate. The new
B-plane coordinates are mapped back to the orbit elements space.

The analytical solution is derived for a hyperbolic flyby around a point
secondary mass. This mapping between B-plane coordinates and orbit elements is
linearized in the impact parameter. Thus, the encounter must be close for the
method to be reliable. Additionally, if the inclination is small the relative
velocity coordinates become undefined. A possible way to avoid this is by
using a method sometimes referred as pseudo-Öpik (Greenberg et al., 1988). In
this case the relative velocity vector is computed directly and defines the
turn angle $\gamma$ at the time of closest approach:

$\tan{\frac{\gamma}{2}}=\frac{m}{bU^{2}}$ (23)

where $m$ is the mass of the planet in units of the mass of the Sun, $b$ is
the impact parameter and $U^{2}$ is the relative velocity in units of the
circular velocity of the planet. Here we use the unperturbed trajectory of the
planet and asteroid to find these quantities. That is, the impact parameter
and relative velocity are found as the planetocentric distance and velocity at
closest approach.

Figure 5: Logarithm of the errors in the computation of the final Keplerian
elements given by pseudo-Öpik theory (PÖpik - first and second rows) and the
quadrature of LPE (QLPE - third and fourth rows). The list of flybys used is
shown in figure 2 and generated as described in the text. The flybys are
represented in the plane of relative velocity at closest approach $V_{inf}$
(km s-1) and distance of closest approach (au). The dashed area in the QLPE
error plots separates the region in which the encounters are computed using an
integration of the 3BP during the semi-analytical propagation.

####  3.4.2 Lagrange Planetary Equations

The proposed computation of the encounter effect is computed as follows. The
variation in elements over the encounter event is obtained from a quadrature
of the Lagrange Planetary Equations assuming the geometry of the unperturbed
flyby. The elements used are obtained from the secular propagation of the
asteroid. The Lagrange Planetary Equations describe the evolution of orbit
elements due to a perturbing potential. The derivation of Lagrange Planetary
Equations can be found in some references for different sets of orbit elements
(Brouwer & Clemence, 1961; Roy, 2004). In general, they have the form of:

$\bm{\dot{D}}=\left[L(\bm{D})\right]\frac{\partial
R}{\partial\bm{D}}(\bm{D},t)$ (24)

where $\bm{D}$ is the set of elements of choice and $L(\bm{D})$ is a function
of the elements that depends on the chosen elements. In this case the
perturbing potential $R$ is the gravitational potential of the encountered
planet. Note that without further simplifying assumptions the partials of the
perturbing potential are function of the elements and function of time. For
the elements of choice we take the partial derivatives that relate the set of
elements to Keplerian elements $\bm{K}=[a,e,i,\Omega,\omega,\sigma]$ and
Cartesian coordinates. From the orbital elements representations available to
choose, the current implementation uses the Delaunay elements:

$\begin{aligned} L&=\sqrt{\mu a}\\\ G&=L\sqrt{1-e^{2}}\\\
H&=G\cos{i}\end{aligned}\qquad\qquad\begin{aligned} l&=\sigma\\\ g&=\omega\\\
h&=\Omega\end{aligned}$ (25)

The Lagrange Planetary Equations with the perturbing potential of Equation 1
with $j$ being the asteroid and $k$ the encountered planet can be written as:

$\begin{aligned} \frac{\mathrm{d}L}{\mathrm{d}t}&=\frac{\partial
R}{\partial{r}}\frac{\partial{r}}{\partial\bm{K}}\frac{\partial\bm{K}}{\partial
l}\\\ \frac{\mathrm{d}G}{\mathrm{d}t}&=\frac{\partial
R}{\partial{r}}\frac{\partial{r}}{\partial\bm{K}}\frac{\partial\bm{K}}{\partial
g}\\\ \frac{\mathrm{d}H}{\mathrm{d}t}&=\frac{\partial
R}{\partial{r}}\frac{\partial{r}}{\partial\bm{K}}\frac{\partial\bm{K}}{\partial
h}\end{aligned}\qquad\qquad\begin{aligned}
\frac{\mathrm{d}l}{\mathrm{d}t}&=-\frac{\partial
R}{\partial{r}}\frac{\partial{r}}{\partial\bm{K}}\frac{\partial\bm{K}}{\partial
L}\\\ \frac{\mathrm{d}g}{\mathrm{d}t}&=-\frac{\partial
R}{\partial{r}}\frac{\partial{r}}{\partial\bm{K}}\frac{\partial\bm{K}}{\partial
G}\\\ \frac{\mathrm{d}h}{\mathrm{d}t}&=-\frac{\partial
R}{\partial{r}}\frac{\partial{r}}{\partial\bm{K}}\frac{\partial\bm{K}}{\partial
H}\end{aligned}$ (26)

The proposed solution is the integration of these differential equations
around the encounter date $t_{e}$ and assuming the unperturbed geometry of the
flyby. Hence, the asteroid coordinates are obtained from the heliocentric
elements secularly propagated until the start of integration date $D_{0}$ and
the quadrature is only a function of time:

$\Delta{\bm{D}}=\int_{t_{e}-\delta t}^{t_{e}+\delta
t}\bm{\dot{D}}(\bm{D}_{0},t)dt$ (27)

The integration is conducted for a fraction of the orbit period around the
closest approach distance. This fraction is a constant set large enough such
that the effect of the encounter is captured completely. In this work we use a
self-coded fast quadrature function based on the midpoint rule and a total
integration time of 20% of the orbital period. This method avoids the frame
transformation to the center of the planet, since it considers the planet as
an external perturber of the asteroid motion around the Sun. For this reason,
it is not possible to obtain a closed form solution of the integral.
Nonetheless, this approach does not imply further assumptions that limit its
range of applicability. Future work will be done in finding the optimal set
for this application. This approach is accurate for the vast majority of
encounters, but it is less accurate for the closest ones, as we explore in the
following section.

Using the list of flybys generated in figure 2 we computed the errors of Öpik
theory and the quadrature of LPE compared to the solution of the encounter
using the three-body problem. The error $E(K)$ is relative to the variation
and in percent, given by:

$E(K)=100\frac{\Delta K_{QLPE}-\Delta K_{3BP}}{\Delta K_{3BP}}$ (28)

The results of this evaluation are described in figure 5. Using pseudo-Öpik
theory there is a region in the space of relative velocity and closest
approach distance in which flybys can be computed accurately. However, this
region is not constant for all Keplerian elements. In addition, most flybys
are not computed correctly using this method. Slow flybys break the assumption
in Öpik theory that the behavior during the flyby is modeled by the two-body
hyperbolic interaction. Many of the faster flybys occur with Venus and
Mercury. The two-body hyperbolic flyby model fails to characterize the effect
of flybys with these less massive planets even if they are faster.

Using the quadrature of the Lagrange Planetary Equations 99% of the flybys
list are computed with less than 3% of error, and more than 88% with less than
0.1% error. The flybys that are not computed accurately with this method are
very close and with a slow relative velocity. These flybys can break the
assumption of the unperturbed geometry of the flyby. Given that these
encounters also cause significant variations in the elements, these infrequent
encounters are solved using a three-body problem integration in Cartesian
coordinates. The criteria to solve these encounters using the alternative
method is by defining three threshold regions in the encounter parameters:
with very small $V_{\infty}$, very small closest approach distance, and a
combination of both close to zero. This process simplifies the detection of
collisions with the planets during the numerical integration in Cartesian
space in the heliocentric frame.

###  3.5 Semi-analytical propagation vs. numerical integration

In the previous sections we validate the individual pieces of the semi-
analytical propagation tool. Once combined, we want to compare the resulting
trajectories to trajectories obtained using numerical integration. With this
purpose we generate a fictitious NEO population and propagate their orbits
using both methods.

The fictitious NEO population we define consists in normal distributions for
the perihelion distance, eccentricity, inclination. The distributions are
centered respectively around $0.8$au, $0.2,10\deg$ and with standard
deviations of $0.05$au, $0.05,3\deg$. Arguments of the node, perihelion and
initial mean anomalies are defined as uniform distributions in the $[0,360]$
degrees range. We sample 1000 particles from these distributions as our test
set for the comparison between the two methodologies. We setup the numerical
integration of the asteroid orbits considering the planets as third body
perturbers. The model we use for the orbits of the planets is the secular
theory developed in section 3.1.

Figure 6 shows the results of the propagation using both methods as well as
the distributions that defined the initial conditions. In addition, we find
that the cumulative distribution of the mean number of encounters versus
closest approach distance is matched very closely. After 200,000 years, the
presence of planetary encounters causes a dispersion of the initial
distributions. The distribution obtained through semi-analytical propagation
is able to track this drift.

The main difference between the results using the two methods is that the
numerically integrated distribution shows a small drift in the center of the
distribution, but very similar dispersions. In terms of the longitude of the
perihelion, the resulting distributions after 200,000 years remain uniform
using either of the two propagation methods. The other significant difference
between the two methods is in the required the computational time, which is
discussed next.

Figure 6: Comparison between the 1000 trajectories obtained using numerical
integration (blue) and semi-analytical propagation (orange) of a fictitious
population of NEOs after 200,000 years. The probability density function is
shown for the perihelion distance, inclination and argument of perihelion. The
analytical probability density function of the initial conditions is shown as
a continuous red line. On the right, the mean number of encounters found as
function of the closest approach distance for both methods.

###  3.6 Computational time of the semi-analytical propagation tool

The semi-analytical propagation of near-Earth asteroids reduces the
computational time required to obtain long-term trajectories. The use of
numerical techniques is limited almost exclusively to the computation of
planetary encounters, that represent a fraction of the simulation time.

In order to estimate the speedup of the propagation we generate a fictitious
population of NEOs and propagate them for 100,000 years with the semi-
analytical propagation tool and with numerical integration. The semi-
analytical propagation tool uses the secular solution of the Solar System to
compute the orbits of the planets over long periods of time. Then, planetary
encounters can occur with any planet. The secular propagation of the asteroid
is derived as described in the text accounting only for Jupiter, which
accurately represents the asteroid motion in the absence of resonances. The
semi-analytical propagation of 100,000 years with this setup is computed in
around 5 seconds. In order to account for these effects in numerical
integrations we use the N-Body problem integrator IAS15 (Rein & Spiegel, 2014)
with all the planets and the asteroid.

The simulation is setup using high-level programming code that runs libraries
in more efficient low-level code. This is the case for both numerical
integrations and semi-analytical propagation, running in the same 2.5GHz Intel
Core i7 processor. The result is a speed-up of x500-x1000 of the semi-
analytical propagation tool as compared to the numerical integration. The
current implementation allows room for significant speed-up that is left for
future work. The perturbed long-term propagation could be extended to use
other suitable models of interest. The computational cost is not expected to
increase while we use analytical solutions of these long-term perturbations.

##  4 Semi-analytical propagation results

In this section we demonstrate the semi-analytical propagation tool in a
variety of scenarios. First, we want to compare the semi-analytical model with
trajectories obtained through numerical integration. Matching very accurately
trajectories obtained with more complex models is outside the scope of the
comparison. Even if the models were identical, trajectories under encounters
are very sensitive to the initial conditions and under small perturbations
they diverge into different paths. This effect was visualized in figure 1
using only numerical integration. For this reason, long term simulations may
focus on the statistical analysis of the dynamical evolution rather than
individual trajectories. Throughout the section, the simulations include
Jupiter as the only planet secularly perturbing the asteroids. All the inner
Solar System bodies are considered to evaluate planetary encounters. These are
secularly evolving due to mutual perturbations and the perturbations of the
outer Solar System planets, as described in section 3.a..

###  4.1 Short-term propagation of near-Earth objects using different models

Figure 7: Trajectories of a particle that mirrors the binary asteroid (35107)
1991 VH for 500 years. Comparison of three methods. First, numerical
integration of the heliocentric orbit of the asteroid with two external
perturbers in constant sets of elements: Earth and Jupiter (4BP, in red). Full
Ephemeris as accessed in HORIZONS (Giorgini & JPL Solar System Dynamics, 2021)
(green), and the semi-analytical propagation method (SA-prop, blue).

The asteroid chosen for the tool demonstrations is the binary (35107) 1991 VH.
The first reference trajectory is obtained from the HORIZONS system of JPL
(Giorgini & JPL Solar System Dynamics, 2021). The second model is the
numerical integration of the asteroid motion under the influence of the Sun,
Earth and Jupiter. Jupiter is the main driver of the secular motion, which is
observed as a linear drift in the argument of perihelion and argument of the
ascending node. Given the current orbit of (35107) 1991 VH, it only
experiences planetary encounters with Earth in the next few centuries. In
Figure 7, we compare these two models from numerical integration with the
present semi-analytical propagation tool.

In the three trajectories we observe similar behavior, although it manifests
differently in every element. First, there is a close agreement in the
encounter dates of the Sun-Earth-Jupiter integration and the encounters found
by the propagation tool that shows in all elements. The encounter dates can be
distinguished as the discontinuities in the trajectories, especially in the
semi-analytical propagation trajectory. The variation of the semi-major axis
is characterized by a random walk from the planetary encounters. In both
eccentricity and inclination there is a relevant role of the encounters with
an additional secular component that the secular model is able to model.

The dynamics of the argument of the ascending node are dominated by the
secular drift. There is not a significant effect that can be perceived by the
planetary encounters and this is expected when the encounters occur with a
unique planet and close to the node. What we observe is that the secular
dynamics including the complete effects of Earth and Jupiter are very similar,
and the analytical secular drift is off by about a degree after the 500 years
of the propagation. The secular drift rate in the argument of perihelion is
not as trivial to compare since it must be observed between encounters,
although good agreement is found too. Last, the mean anomaly at epoch evolves
over time with an increasing amplitude present in all three models.

Figure 8: Close encounters in the semi-analytical propagation of the
fictitious NEOs Cases 1-3 with initial conditions given in table 2. Closest
approach distance $d_{CA}$ (au) of the encounters and MOID with the planets.
Dots indicate close encounters, solid lines indicate the MOID. Color code
indicates the planets: (Green) Venus, (Blue) Earth, (Red) Mars. Table 2: Near-
Earth objects used as example cases for the demonstration of the semi-
analytical propagation tool.

Asteroid | $a(au)$ | $e$ |  $i$ (deg) |  $\Omega$ (deg) |  $\omega$ (deg) |  $M_{0}$ (deg) | JD (TBD)
---|---|---|---|---|---|---|---
Case 1 | 1.1 | 0.15 | 10 | 90 | 90 | 90 | 2451545.0
Case 2 | 1.2 | 0.35 | 40 | 90 | 90 | 90 | 2451545.0
Case 3 | 1.3 | 0.5 | 10 | 90 | 270 | 90 | 2451545.0
Case 4 | 0.95 | 0.07 | 20 | 90 | 90 | 90 | 2451545.0
Case 5 | 0.9 | 0.25 | 15 | 90 | 90 | 90 | 2451545.0
1991 VH | 1.1373 | 0.14426 | 13.912 | 139.37 | 206.88 | 302.39 | 2456902.5
1996 FG3 | 1.0543 | 0.34987 | 1.9903 | 299.88 | 23.930 | 147.277 | 2454796.5

Note. — The elements of asteroids 1991 VH and 1996 FG3 were retrieved from
HORIZONS (Giorgini & JPL Solar System Dynamics, 2021). Using DE431 and
SB431-N16. 1991 VH: Orbit solution date 2021 April 15, 1996 FG3: Orbit
solution date of 2021 April 26.

###  4.2 Long-term propagation and the MOID

The minimum orbit intersection distance (MOID) indicates what the minimum
distance between any two points of the two heliocentric Keplerian orbits is.
In this case we focus on the orbit of Earth and the orbit of the asteroid. The
MOID is also used as one of the criteria to define an asteroid as a
potentially hazardous asteroid. There are many algorithms available in the
literature to compute the MOID (Gronchi, 2005; Wiśniowski & Rickman, 2013;
Armellin et al., 2010). In this paper we use the tool derived in Hedo et al.
(2018, 2020) based on an asymptotic approach.

The MOID constrains the minimum distance of a possible close encounter. In
other words, the periods with a large MOID are absent of close encounters.
Three examples are used to visualize time histories of close encounters and
the evolution of the MOID for 100,000 years. These are obtained for high and
low eccentricities and inclinations, and the initial conditions are found in
table 2.

Figure 9: Close encounters in the semi-analytical propagation of the
fictitious NEO Cases 1-3 with initial conditions given in table 2.
$V_{\infty}$ at closest approach of the encounter (km s-1). Color code
indicates the planets: (Green) Venus, (Blue) Earth, (Red) Mars.

The distributions of closest approach distances are shown in figure 8 and the
unperturbed relative velocity $V_{\infty}$ at those encounters is found in
figure 9. Case 1 is an example of a NEO with relatively low eccentricity and
inclination. In these conditions, close encounters are only possible with
Earth and at a low relative velocity. The MOID oscillates secularly with long
periods of low MOID. Case 2 is an example of an opposite scenario in which
both eccentricity and inclination are large. In the secular evolution of the
MOID this translates in short periods of low MOID and long periods absent from
encounters. This is scenario in which the semi-analytical propagation of the
asteroid allows a rapid propagation until the next period of frequent
encounters. Case 2 faces high velocity close encounters with Venus, Earth and
Mars.

Case 3 is an example of a NEO with high eccentricity and low inclination. This
combination of factors leads to a large frequency of close encounters with the
inner planets. In this case, encounters are very frequent with Venus, Earth
and Mars. The close encounters experienced by Case 3 are with a relative
velocity smaller than in Case 2 given the reduced inclination. Even under the
elevated frequency of close encounters, the secular signature of the MOID is
maintained. The structure shows until the event of an energetic close
encounter. The occurrence of such encounters is just a matter of probability
of having the right timing during the low-MOID intervals of the secular
propagation.

###  4.3 Statistics of long-term propagation

The chaotic nature of the dynamics implies that the study of the orbital
evolution over long timescales should be done statistically. Given the
uncertainty in the orbit solution of an asteroid, we can sample a large number
of particles and study the dynamical paths that the different particles take.
Because of the sensitivity to initial conditions in planetary encounters, very
well determined distributions diverge in a few centuries to widely different
paths. We demonstrate these effects by propagation of 500 particles that
sample uncertainty distributions around a nominal asteroid orbit for 500,000
years. We inspected 7 examples: first, in detail, orbit solution of (35107)
1991 VH; then, the orbit solution of (175706) 1996 FG3; last, the previous
Cases 1-5 of table 2 with artificial orbit uncertainties as described in
appendix A.

Figure 10: 500,000 year Monte Carlo semi-analytical propagation of asteroid
(35107) 1991 VH. Initial conditions are given in table 2 as obtained from
HORIZONS (Giorgini & JPL Solar System Dynamics, 2021) and uncertainties in the
distribution are obtained from JPL’s SSD/CNEOS Small-Body DataBase (SSD/CNEOS
API Service, 2021) as described in Appendix A. Elements shown are perihelion,
eccentricity, inclination, argument of the ascending node, argument of
perihelion and minimum orbit intersection distance (MOID). Grey lines show
individual simulations, black lines are the median of the 500 simulations of
each parameter shown.

Figure 11: Semi-analytical propagation of asteroid Cases 1-5, (35107) 1991 VH,
(175706) 1996 FG3 for 500,000 years. Histograms at the initial time, after
100,000 years and at the final time. Initial conditions are given in table 2
as obtained from HORIZONS (Giorgini & JPL Solar System Dynamics, 2021) and
uncertainties in the distribution are obtained from JPL’s SSD/CNEOS Small-Body
DataBase (SSD/CNEOS API Service, 2021) as described in Appendix A.

Figure 10 shows the time histories of the individual runs of the cloud of
points originally neighboring (35107) 1991 VH, showing that the cloud of
particles distributes over a wide region of near-Earth space. On the order of
hundreds of thousand years, the dispersion is accomplished by the less
frequent very close encounters. Eccentricity and inclination show the secular
component but are dispersed by the presence of encounters. The dynamics of
argument of ascending node and argument of perihelion remain mainly secular
with a degree of dispersion because of the presence of encounters. Note that
the initial uncertainty on the orbit of (35107) 1991 VH is very small as shown
while demonstrating the chaotic dynamics nature in Figure 1 of the background
section. As it was observed in the detailed analysis of a shorter simulation
in figure 7, the mean anomaly at epoch changes completely with small changes
in the semi-major axis. This fact reflects in the long-term simulations as a
complete uniformization after just a few centuries.

The binary (35107) 1991 VH currently presents a MOID that is decreasing. This
means that after a few millenia the probability of experiencing very close
encounters increases. In the statistical analysis, this probability shows in
that a fraction of the fictitious asteroids experience such encounters. We
observe that towards the end of the simulation there is a large dispersion in
inclination and perihelion distance. By the end of the simulation, the angles
$\Omega-\omega$ are dispersed along a linear drift as caused by numerous close
encounters that not only change these angles, but modify the secular dynamic
frequencies.

The long-term dynamics of 6 more examples are integrated for 500,000 years.
These are the Cases that we used to illustrate the long-term dynamics with
initial conditions in Table 2. The case of the binary (175706) 1996 FG3 is
also added to the discussion, as it is the other target of exploration of the
Janus mission (Scheeres et al., 2020). The uncertainties of (175706) 1996 FG3
and (35107) 1991 VH are sampled based on their publicly available orbit
solutions. In the case of the fictitious asteroids we used an arbitrary
distribution. Both approaches are explained in detail in the appendix A. The
statistical distributions after 100,000 and 500,000 years are shown in Figure
11. The combined final and initial distributions are shown in Figure 12. The
recorded number of encounters are shown in Figure 13 and more details on the
statistical distributions over time are shown in Figure 14. Next, we describe
the dynamical evolution of these test cases.

In Figure 12, the 500 virtual asteroids that are generated at the initial time
per case are in the same bar of the histogram. The effect of repeated
encounters causes the distributions to spread along near-Earth space. This
dispersion is clearly shown in the perihelion distance in all cases, with a
general trend of a decrease in the distance. Figure 12 additionally shows this
spread of all the cases together.

The presence of mean motion resonances in the inner Solar System can protect
asteroids from close encounters Milani et al. (1989). In the semi-analytical
propagation of near-Earth asteroids the orbits may drift to these regions, as
observed in Figure 12. The resonance regions are found for semi-major axes
larger than $a=1$ au. In these cases encounters with the Earth stop occurring
for a period of time. However, this clustering of particles is not found in
numerically integrated populations or in the discovered population of near-
Earth asteroids. Thus, we investigate the trajectories obtained through NBP
dynamics.

The process of obtaining the secular long-term perturbation eliminates the
short-period perturbations. The latter perturbations cause an oscillation in
the orbit elements of non-negligible amplitude. In order to measure the
influence of this effect, we included an analytical oscillation in the semi-
major axis with the frequency of the orbital period and an amplitude of the
order of 0.01 au. This extension completely eliminates the artificial mean
motion resonance regions. The analytical characterization of the short-period
perturbation is left as future work, as its contribution must be considered
for the complete set of orbital elements and is not expected to significantly
modify the obtained distributions.

Figure 12: Dispersion in the distributions of the 5 asteroid in table 2,
(175706) 1996 FG3 and (35107) 1991 VH after 100,000 years (top) and 500,000
years (bottom). Semi-major axis vs. eccentricity (left) and semi-major axis
vs. inclination (right). Initial uncertainty distributions detailed in
appendix A and shown in red. See text for the explanation of the resonances
found in the semi-analytical propagation.

Figure 13: Number of encounters experienced by the 5 test cases given in table
2, (175706) 1996 FG3 and (35107) 1991 VH during 500,000 years. Mean of the
distributions of 500 particles (solid lines) and 1-$\sigma$ bounds (shadowed
area).

In general, we observe that asteroids encountering the planets more frequently
disperse their distributions faster. This is the case for the Janus targets
and Case 1, that present the largest standard deviations increase in the
simulation time (Fig. 14). This dispersion is not only shown in the elements,
but also in the number of encounters (Fig. 13).

The evolution of the distribution as caused by encounters could be modelled as
a random walk. If this hypothesis is true, then the standard deviation in the
population increases linearly with the square root of time. Figure 14 tests
graphically this hypothesis for semi-major axis, eccentricity and inclination.
Initially in all cases there is a fast increase in the standard deviations.
After the few first millennia, some of the distributions follow the hypothesis
of the linear relationship $\sigma\propto\sqrt{t}$, especially in the semi-
major axis.

In eccentricity and inclination, we observe that there is a secular component
in the evolution of the distribution. The secular component is observed in
both the evolution of the standard deviation and the mean of the variations
(Fig. 14). These are shown with respect to the initial values to have a common
reference in the comparison of cases.

The dispersion in semi-major axis modifies the secular rates of the drift in
argument if perihelion and argument of the node. This effect combined with the
direct variation of the angles during planetary encounters leads to the
uniformization of the distribution in $\omega,\Omega$. In the duration of our
simulations of 500,000 most cases approach this uniformization as we showed in
figure 11.

Figure 14: Statistical evolution of the distributions of the 5 test cases
given in table 2, (175706) 1996 FG3 and (35107) 1991 VH during 500,000 years.
Standard deviation of semi-major axis, eccentricity, inclination as function
of the number of encounters (top row), square root of time (center row).
Variation of the mean value over time with respect of the initial value of
semi-major axis, eccentricity and inclination (bottom row).

We can compute more rigorously whether the distributions can be considered
uniform by conducting the chi-squared test on the longitude of perihelion
$\varpi=\omega+\Omega$. Figure 15 shows the result of the chi-squared test of
a uniform distribution over the simulation time for all cases. If the p-value
is larger than our threshold of $p=0.05$, we can consider that our null
hypothesis of the uniform distribution of $\varpi$ is true. Cases 1, 3, 4,
(175706) 1996 FG3 and (35107) 1991 VH reach this threshold, while Cases 2 and
5 do not approach the significance by the end of the simulation time. It is
interesting to show the evolution of the p-value compared to the mean number
of encounters. In figure 15 we show how the distribution of Case 3 tends to
the uniformization in $\varpi$ with less encounters than other cases that
achieve this distribution earlier in the simulation time. This is expected
since this case experiences more frequent encounters with the most massive
planets and with a slower relative velocity, which means that the impact of
these encounters in the dispersion of their distributions is larger.

From the general trends that we observed, the only case that is very different
is Case 2, that experiences much fewer encounters. As we illustrate in Figure
8, Case 2 experiences close encounters much less frequently than the other
cases. The relative velocity is also larger in this case, which means that the
effects of the encounters are not as strong. The case of (175706) 1996 FG3 is
more difficult to fit in the general description of the dynamics, as the close
encounters do not cause such a fast dispersion of the distribution. This
binary asteroid is studied in more detail in section 6.

Figure 15: P-value of the chi-squared test of the uniform distribution of the
longitude of perihelion $\varpi=\omega+\Omega$ for the 5 test cases given in
table 2, (175706) 1996 FG3 and (35107) 1991 VH during 500,000 years. P-value
is shown as function of the mean number of encounters (top) and function of
time (bottom).

##  5 Discussion

The semi-analytical propagation tool shows the main dynamical effects observed
in long-term numerical integration of the inner Solar System. The secular
drifts caused by Jupiter move the asteroids between the vicinities of the
different planets of the inner Solar System. This effect causes a seasonal
presence of strong close encounters that can disturb both the orbit of the
asteroid and its physical properties. While the time-scales of these events is
of millions of years (Fang & Margot, 2012), if we sample a large enough number
of particles we can measure the probabilities of collisions or the potential
disruption of other physical properties of asteroids. The measurement of the
collision probabilities was outside the scope of this paper, but in this work
a few collisions were detected in the uncertainty sampling of the asteroids.

The analytical modelling of the dynamics far from the planets was done using
the Laplace-Lagrange theory, which works well in a large fraction of the NEO
population. For this reason we defined a region in near-Earth space in which
the secular model works best, as shown in Figure 4. However, we could extend
the modelling in the regions of large eccentricity and inclination. In the
previous section we describe the low frequency of encounters that is
characteristic of asteroids with high inclination, specifically with Case 2.
An asteroid of these characteristics would be likely to be dominated by the
Lidov-Kozai effect, in which there is an exchange between high inclination-low
eccentricity periods and low inclination-high excentricity periods. This would
mean that Case 2 evolves to become a case closer to Case 3, in which
encounters are more frequent. The use of analytical models of the Lidov-Kozai
model (Kinoshita & Nakai, 2007) for the perturbed propagation is left for
future work.

Using the semi-analytical propagation tool we observe the stochastic nature of
the dynamics. However, the effect is different on each of the elements. While
in semi-major axis we observe what could be described as a random-walk
process, the angles $\Omega$ and $\omega$ become uniformly distributed.
Eccentricity and inclination show a mixed effect between a random-walk that
adds dispersion to the distribution and the oscillations driven by the secular
theory.

##  6 Orbit histories of the Janus mission targets

The binary asteroids (175706) 1996 FG3 and (35107) 1991 VH are the targets of
the exploration mission Janus (Scheeres et al., 2020). The long-term orbital
dynamics of the asteroids are studied in this section with two goals. First,
characterizing the stochastic nature of their long-term dynamics. Then,
estimating the probability that they have been potentially disrupted by a very
close encounter. With these purposes, we study their orbital origins by
sampling 1000 particles and propagating them for 1Myr into the past using the
semi-analytical propagation tool.

Figure 16: Orbit history of (35107) 1991 VH in the last million years. Initial
conditions are given in table 2 as obtained from HORIZONS (Giorgini & JPL
Solar System Dynamics, 2021) and uncertainties in the distribution are
obtained from JPL’s SSD/CNEOS Small-Body DataBase (SSD/CNEOS API Service,
2021) as described in Appendix A. Elements shown are perihelion distance,
eccentricity, inclination, argument of the ascending node, argument of
perihelion and minimum orbit intersection distance (MOID) with Venus, Earth
and Mars. Grey lines show individual simulations, black lines are the median
of the 1000 simulations of each parameter shown. Only 50 example runs and the
median of the full distribution are shown in orbit elements (Top 2 rows).

The orbit time histories of (35107) 1991 VH and (175706) 1996 FG3 are shown
respectively in figure 16 and figure 17. For clarity, we show only a subset of
the runs and the median of the full distribution of 1000 runs. The minimum
orbit intersection distance (MOID) is shown for the inner Solar System planets
Venus, Earth and Mars. These metrics show when close encounters with these
planets are possible. The presence of frequent close encounters causes the
dispersion of the orbit histories. This feature manifests in the orbit history
of (175706) 1996 FG3, in which the period of very low Venus MOID corresponds
with a dispersion in the overall statistical representation of the orbit.

Figure 17: Orbit history of (175706) 1996 FG3 in the last million years.
Initial conditions are given in table 2 as obtained from HORIZONS (Giorgini &
JPL Solar System Dynamics, 2021) and uncertainties in the distribution are
obtained from JPL’s SSD/CNEOS Small-Body DataBase (SSD/CNEOS API Service,
2021) as described in Appendix A. Elements shown are perihelion distance,
eccentricity, inclination, argument of the ascending node, argument of
perihelion and minimum orbit intersection distance (MOID) with Venus, Earth
and Mars. Grey lines show individual simulations, black lines are the median
of the 1000 simulations of each parameter shown. Only 100 example runs and the
median of the full distribution are shown in orbit elements (Top 2 rows).

Figure 18: Histograms of the orbit history of (35107) 1991 VH and (175706)
1996 FG3 at the initial time, 100,000 years ago and 1Myr ago. Initial
conditions are given in table 2 as obtained from HORIZONS (Giorgini & JPL
Solar System Dynamics, 2021) and uncertainties in the distribution are
obtained from JPL’s SSD/CNEOS Small-Body DataBase (SSD/CNEOS API Service,
2021) as described in Appendix A.

Similarly to the long-term dynamics into the future studied in section 4, the
initially very close distribution becomes a wide statistical distribution when
propagated far into the past. In Figure 18 we show histograms of the orbit
elements and the number of encounters recorded below a closest approach
distance threshold of 0.1 au. The orbit evolution of (35107) 1991 VH is mostly
a spread around the initial conditions. However, (175706) 1996 FG3 is in a
particular initial orbit with low inclination. On average, the very low
inclination and high eccentricities drift toward a smaller eccentricity and
higher inclination that are more frequent in the secular cycle. In both cases,
the longitude of the perihelion becomes uniformly distributed. In the next
section we characterize this uniformization process.

###  6.1 Stochastic modelling of the long-term dynamics

In section 4 we show that we can model the long-term dynamics with a random
walk in semi-major axis, eccentricity and inclination. In addition, the latter
two present also the influence of the oscillations of the secular theory. We
also want to study the uniformization in the longitude of the perihelion, as
this process occurs with time but also with a repeated number of close
encounters.

The random-walk model can be characterized by a linear increase in standard
deviation with the square root of time. Figure 19 shows the standard deviation
of the 1000 Monte Carlo experiment that we conducted into the past of the two
Janus targets (35107) 1991 VH and (175706) 1996 FG3. We fit a linear model to
the standard deviation evolution, and measure the slopes to compare the
evolution of the two targets. In the case of (35107) 1991 VH we avoid the use
of the full simulation time, as the standard deviation bends from the initial
linear increase. The slower increase after this period occurs when the
distribution migrates from a configuration with slower encounters. The
opposite case occurs with inclinations, in which the rapid increase of
(175706) 1996 FG3 from the low-inclination initial regime slows down after the
initial growth.

The measured slopes are reported in table 3, showing that the random walk of
(175706) 1996 FG3 is faster in semi-major axis and inclination. However,
because of the bends in the progression after a few hundred thousand years,
the final standard deviations are not substantially larger than the ones of
(35107) 1991 VH after the million years into the past in eccentricity and
inclination.

Figure 19: Random walk statistical modelling of the evolution of semi-major axis, eccentricity, inclination of (35107) 1991 VH (top) and (175706) 1996 FG3 (bottom). Standard deviation of the 1000 Monte Carlo runs as function of the square root of time and random walk model fit. Table 3: Stochastic modelling of the long-term dynamics of the Janus targets Target | Random walk $k_{x}$ constant ($s_{x}=k_{x}\sqrt{t}$) |  $\varpi$ Uniformization
---|---|---
Name | $k_{a}$ | $k_{e}$ | $k_{i}$ | Time | Encounters $<0.1$ au
| (au/$\sqrt{yr}$) $\cdot 10^{-3}$ | (1/$\sqrt{yr}$) $\cdot 10^{-3}$ | (deg/$\sqrt{yr}$) $\cdot 10^{-3}$ | $yr\cdot 10^{3}$ | Number
(35107) 1991 VH | 0.2661 | 0.1799 | 6.0055 | -434 | 11800
(175706) 1996 FG3 | 0.3688 | 0.0944 | 11.598 | -797 | 29800

The process of uniformization of the longitude of perihelion is shown in
figure 20. We conduct the chi-squared test of the uniform distribution over
the 1Myr simulation, to find when the hypothesis of the uniform distribution
is significant. In table 3 we show the first time in which this criterion is
satisfied, both in time and mean number of encounters: -434,000 years and a
mean of 11800 encounters for (35107) 1991 VH, and -797,000 years and a mean of
29800 encounters for (175706) 1996 FG3.

The uniformization of (35107) 1991 VH is faster than the uniformization of
(35107) 1991 VH in both time and mean number of encounters. It is remarkable
that (35107) 1991 VH takes a much lower mean number of encounters. This is
explained by the faster relative velocities of the encounters of (175706) 1996
FG3 and a larger fraction occuring with Mars, a less massive planet. The
relative velocities of a few of the recorded flybys are shown in figures 21
and 22 in the context of studying the probability that a close encounter could
potentially disrupt the binaries.

The comparison between the two binary systems highlights how the effect of the
encounters depends on the relative velocities and the mass of the planets. In
general, slow encounters and with larger planets are more efficient at causing
the distributions to become uniform. However, depending on the heliocentric
orbit these encounters may be more or less frequent. Thus, leveraging both
effects is required to obtain a stochastic representation of the long-term
dynamics of NEOs under frequent encounters.

Figure 20: P-value of the chi-squared test of the uniform distribution of the
longitude of perihelion $\varpi=\omega+\Omega$ of (35107) 1991 VH (top) and
(175706) 1996 FG3 (bottom). P-value is shown over the past 1Myr (left) and as
function of the mean number of encounters (center). The mean number of
encounters of the 1000 Monte Carlo runs is shown over time the corresponding
with 1-$\sigma$ bounds (right).

###  6.2 Potentially disruptive planetary encounters

The two binary targets of the Janus mission present different relative orbits
as observed by radar and photometry (Pravec et al., 2016; Meyer et al., 2021),
showing that (175706) 1996 FG3 is in a stable orbital state and (35107) 1991
VH is in a chaotic state. The perturbed state of (35107) 1991 VH could be
explained by a recent very close encounter with the inner Solar System planets
(Heggie & Rasio, 1996). Thus, it is of interest to characterize the frequency
of such encounters in the orbital history of asteroid binaries.

Using the semi-analytical propagation tool we obtain the history of flybys
over a long time period of time. The perturbation in the orbit of a binary
system during a planetary close encounter is studied in detail as described in
(Meyer & Scheeres, 2021). In this section we combine both results to predict
the frequency of a disrupting flyby.

The effect of the close encounter on the binary can be modelled as an
impulsive variation in the binary Keplerian elements. In Meyer & Scheeres
(2021) the effect of close encounters to singly synchronous binary asteroids
is studied. The variation in semi-major axis, eccentricity, and inclination
obtained with numerical methods was compared to analytical expressions for the
impulsive variation in binary Keplerian elements (Fang & Margot, 2011). We
used these analytical expressions as they provide an estimate of the variation
as function of the relative velocity and distance of closest approach.

Figure 21: Potentially disruptive encounters recorded by the 1000 virtual
(35107) 1991 VH binaries during a million years into the past. The background
contour represents the logarithm of the mean variation in binary eccentricity
during a close encounter with a planet: Earth (blue, top) or Mars (red,
bottom). The encounters found during the propagation in this threshold are
shown by their closest approach distance (au) and $V_{\infty}$ (km s-1). The
radius of the planet is shown as a dashed line. The initial uncertainty
distribution of (35107) 1991 VH is detailed in appendix A.

For every binary and encountered planet we can generate contours of the
variation of the eccentricity. In Figures 21 and 22 we show these contours
respectively for (35107) 1991 VH and (175706) 1996 FG3. In every figure we
highlight two levels, a variation of 0.1 in eccentricity, and a variation of
1, which would mean that the binary is completely separated. The probability
of disruption in the binary orbit increases as the relative velocity and
closest approach distance are reduced.

Using the semi-analytical propagation tool we track the close encounters below
the threshold of 0.003 au, above which the variation in the binary Keplerian
elements becomes negligible. For each binary we generate 1000 trajectories for
a million years into the past. All the encounters that are found in this
threshold are plotted in figures 21 and 22 and separated by encountered
planet.

The encounters potentially disruptive recorded for (35107) 1991 VH with Earth
and Mars are shown in figure 21. Less than 1% of the recorded encounters are
with Venus so the map with this planet is not included. The relative
velocities of the flybys are mostly found between 5 and 15 km s-1. In the case
of (175706) 1996 FG3 this range of possible relative velocities is larger in
all the planets. In addition, many encounters are found with quite slower
$V_{\infty}$, which means that even if they are not as close, they can still
potentially cause a disruption.

As we observed in figure 13, (175706) 1996 FG3 experiences more frequent
encounters. However, the regions in which the encounters are potentially
disruptive depend on the current orbital configurations of the binaries. In
this case, (175706) 1996 FG3 requires closer and slower encounters to obtain
the same mean variation in binary elements.

Figure 22: Potentially disruptive encounters recorded by the 1000 virtual
(175706) 1996 FG3 binaries during a million years into the past. The
background contour represents the logarithm of the mean variation in binary
eccentricity during a close encounter with a planet: Venus (green, top), Earth
(blue, center) or Mars (red, bottom). The encounters found during the
propagation in this threshold are shown by their closest approach distance
(au) and $V_{\infty}$ (km s-1). The radius of the planet is shown as a dashed
line. The initial uncertainty distribution of (175706) 1996 FG3 is detailed in
appendix A.

Considering the orbit history in a million years, a non-negligible probability
exists that both binaries have been potentially disrupted at some point.
However, it is possible that the signature of these potential disruptions
vanish if there is energy dissipation in the system. Thus it is relevant to
study the probability that the binary orbits have potentially been disrupted
recently in the orbit histories. We can study the first fraction of the long-
term secular periods, in which the particles still have not mixed. Figure 23
shows the history of potentially disruptive encounters in the recent periods
of frequent encounters.

The last period of possible very close encounters that we find for (35107)
1991 VH starts beyond the last 10,000 years. As determined by a mean variation
in binary eccentricity of 0.1, we find that 61 of the 1000 test runs
experience a potential disruption in the period of time up until the last
30,000 years. When we incorporate the next period of close encounters, the
probability of a potential disruption increases to 131 out of 1000 runs in the
last 60,000 years.

Given the current orbit state of (175706) 1996 FG3, we find very close
encounters with Venus in the last millennia. Because of the more restrictive
closest approach distance for a potential disruption to happen, only 10 out of
the 1000 test runs experienced this potential disruption in this period of
frequent encounters. If we consider the last 10,000 years, then the
probability increases to 54 cases in which at least one potentially disruptive
encounter was found.

If we keep increasing the time in which we consider all the potentially
disruptive encounters, we can estimate the probability that a potentially
disruptive encounter occurs. These probabilities are shown in figure 23 with
the corresponding 95% confidence intervals. The probability of suffering a
disruptive encounter of (175706) 1996 FG3 increases faster than the
probability of (35107) 1991 VH. This is explained by the significantly higher
number of recorded close encounters. Thus, it is not possible to explain the
chaotic state of (35107) 1991 VH only from the long-term probability of
experiencing such encounters. However, a low probability in recent times
combined with the incapability to dissipate the perturbation in a long time
could explain the chaotic state of (35107) 1991 VH. Thus, future work will be
done in the lines of characterizing the timescales of the dissipation of
perturbations due to close encounters.

Figure 23: (Top) Potentially disruptive encounters time history in the most
recent periods of close encounters of (35107) 1991 VH and (175706) 1996 FG3.
In the most recent 50,000 years, (35107) 1991 VH experiences Earth and Mars
encounters. In the latest 10,000 years, (175706) 1996 FG3 experiences Venus
and Earth encounters. The close encounters below the line of mean variation in
binary eccentricity of 0.1 are highlighted with a black circle and the dashed
line marks the radius of Earth. (Bottom) Probability of disruption based on
the number of encounters found below thresholds of mean $\Delta e=0.1$ for low
disruption and $\Delta e=1$ for high disruption. Confidence intervals on the
probability of disruption prediction are shown in dashed lines.

##  7 Conclusions

In this paper we present a rapid semi-analytical propagation tool for
asteroids in the inner Solar System. The tool combines an analytical solution
for the secular dynamics and the evaluation of planetary encounters. The
derived solution of planetary encounters proves to accurately model the effect
of the majority of flybys that asteroids experience in the inner Solar System.

The long-term effect of the perturbation by Jupiter is captured by the
analytical secular solutions in a large fraction of the NEO population. The
combination with detection and evaluation of close encounters recreates the
full dynamics as we demonstrate for the case of (35107) 1991 VH.

The description of the orbits of NEOs in long-term timescales must be done
statistically. We showed how the different elements can be represented by
different distributions, and how the time it takes for the elements to become
statistical depends on the frequency of close encounters. Through the sampling
of different NEO cases we inspect the stochastic models that represent the
long-term evolution of the orbital elements.

The use of a fast semi-analytical propagation tool allows an efficient study
of the dynamics of asteroids. We studied in detail the orbital histories of
the Janus mission targets: (35107) 1991 VH and (175706) 1996 FG3. We
characterized the encounters that can cause a potentially disruption of the
binary orbits and computed the frequency of such encounters.

Additional modeling of the effects of close encounters to other physical
properties of asteroids will allow the study of the frequency of disruptive
events. These are just a few potential examples of the benefits of a fast
propagation tool for Solar System studies in the fashion of the presented
tool.

##  Appendix A Initial Conditions and uncertainties

The statistical representation of the uncertainty in the orbit solutions can
be done using the covariance matrix. This information is available for
multiple asteroids in JPL’s SSD/CNEOS Small-Body DataBase (SSD/CNEOS API
Service, 2021). The covariance matrices for (175706) 1996 FG3 and (35107) 1991
VH that are used in this work are:

Table 4: Initial covariance of the orbit of NEO binary (35107) 1991 VH

1991 VH | e | q (au) |  $t_{p}$ (TDB) |  $\Omega$ (deg) |  $\omega$ (deg) |  $i$ (deg)
---|---|---|---|---|---|---
e | 3.0691e-16 | -3.5095e-16 | -8.8918e-14 | -7.2651e-15 | -5.2217e-14 | 3.6255e-15
q (au) | -3.5095e-16 | 4.0175e-16 | 1.0484e-13 | 8.2268e-15 | 6.1722e-14 | -4.1006e-15
$t_{p}$ (TDB) | -8.8918e-14 | 1.0484e-13 | 7.5479e-11 | -7.8662e-12 | 6.4146e-11 | -3.8102e-12
$\Omega$ (deg) | -7.2651e-15 | 8.2268e-15 | -7.8662e-12 | 3.104e-11 | -3.4282e-11 | -4.876e-12
$\omega$ (deg) | -5.2217e-14 | 6.1722e-14 | 6.4146e-11 | -3.4282e-11 | 8.1951e-11 | 1.9766e-12
$i$ (deg) | 3.6255e-15 | -4.1006e-15 | -3.8102e-12 | -4.876e-12 | 1.9766e-12 | 8.4244e-12

Note. — As obtained from JPL’s SSD/CNEOS Small-Body DataBase (SSD/CNEOS API
Service, 2021). Using DE431 and SB431-N16, orbit solution date 2021 April 15
for epoch JD = 2456902.5.

Table 5: Initial covariance of the orbit of NEO binary (175706) 1996 FG3

1996 FG3 | e | q (au) |  $t_{p}$ (TDB) |  $\Omega$ (deg) |  $\omega$ (deg) |  $i$ (deg)
---|---|---|---|---|---|---
e | 1.2362e-16 | -1.3026e-16 | 7.9966e-15 | -4.1649e-14 | 3.9787e-14 | -2.3726e-14
q (au) | -1.3026e-16 | 1.3752e-16 | -9.0726e-15 | 4.7436e-14 | -4.5397e-14 | 2.5036e-14
$t_{p}$ (TDB) | 7.9966e-15 | -9.0726e-15 | 2.9877e-12 | 6.6638e-11 | -6.6225e-11 | -2.4698e-12
$\Omega$ (deg) | -4.1649e-14 | 4.7436e-14 | 6.6638e-11 | 7.0256e-09 | -6.9647e-09 | -6.8055e-11
$\omega$ (deg) | 3.9787e-14 | -4.5397e-14 | -6.6225e-11 | -6.9647e-09 | 6.9046e-09 | 6.7284e-11
$i$ (deg) | -2.3726e-14 | 2.5036e-14 | -2.4698e-12 | -6.8055e-11 | 6.7284e-11 | 6.7136e-12

Note. — As obtained from JPL’s SSD/CNEOS Small-Body DataBase (SSD/CNEOS API
Service, 2021). Using DE431 and SB431-N16, orbit solution date 2021 April 26
for epoch JD = 2454796.5.

In the case of the artificial cases used to illustrate the long-term dynamics,
we set the covariance matrix to be a diagonal matrix of 1e-6 in the Keplerian
set $\\{a,e,i,\Omega,\omega,\sigma\\}$. While this is orders of magnitude
larger than the uncertainties of (175706) 1996 FG3 and (35107) 1991 VH, the
uncertainty without further observations increases exponentially after only
tens of encounters. Thus, it is adequate for the studies in long-term
simulations of this work.

The individual particles are sampled considering a multidimensional normal
distribution centered around the nominal values shown in Table 2. Then, we use
the Cholesky factorization of the covariance matrices to add the corresponding
perturbation from the nominal.

##  Appendix B Computation of Laplace coefficients

The expansion of the potential requires the computation of Laplace
coefficients, as introduced by Laplace (1798). Brouwer & Clemence (1961);
Murray & Dermott (2000) detail both the expansion and computation of
coefficients. In the case of the expansion in equation 4:

$\left\langle
R_{j}\right\rangle=n_{j}a^{2}_{j}\left[\frac{1}{2}A_{jj}e^{2}_{j}+\frac{1}{2}B_{jj}I^{2}_{j}+\sum^{N}_{\begin{subarray}{c}k=1\\\
k\neq
j\end{subarray}}A_{jk}e_{j}e_{k}\cos{\left(\varpi_{j}-\varpi_{k}\right)}+B_{jk}I_{j}I_{k}\cos{\left(\Omega_{j}-\Omega_{k}\right)}\right]$

The coefficients $A_{jk}$,$B_{jk}$,$B_{jj}$ and $A_{jj}$ are:

$A_{jk}=-n_{j}\frac{1}{4}\frac{m_{k}}{m_{c}+m_{j}}\alpha_{jk}\bar{\alpha}_{jk}b^{(2)}_{3/2}\left(\alpha_{jk}\right)$
(B1)
$B_{jk}=+n_{j}\frac{1}{4}\frac{m_{k}}{m_{c}+m_{j}}\alpha_{jk}\bar{\alpha}_{jk}b^{(1)}_{3/2}\left(\alpha_{jk}\right)$
(B2) $A_{jj}=\sum^{N}_{k=1,k\neq j}B_{jk}$ (B3) $B_{jj}=-\sum^{N}_{k=1,k\neq
j}B_{jk}$ (B4)

The definition of Laplace coefficient is:

$\frac{1}{2}b_{s}^{(k)}(\alpha)=\frac{1}{2\pi}\int_{0}^{2\pi}\frac{\cos{(k\psi)}d\psi}{(1-2\alpha\cos\psi+\alpha^{2})^{s}}$
(B5)

This integral can be rewritten in a series expansion that simplifies the
computation of the Laplace coefficients numerically as function of the rising
factorial or Pochhammer symbol. However, it is found to be computationally
more efficient to compute the quadrature integral above. There are many
recursion and derivative rules that avoid computing the coefficients based on
the definition. These expressions use the nomenclature of $D$ being the
derivative operator $\frac{\mathrm{d}}{\mathrm{d}\alpha}$.

$b_{s+1}^{(j)}=\frac{s+j}{s}\frac{(1+\alpha^{2})}{(1-\alpha^{2})^{2}}b_{s}^{(j)}-\frac{2(j-s+1)}{s}\frac{\alpha}{(1-\alpha^{2})^{2}}b_{s}^{(j+1)}$
(B6)

$b_{s+1}^{(j+1)}=\frac{j}{j-s}\left(\alpha+\frac{1}{\alpha}\right)b_{s+1}^{(j)}-\frac{j+s}{j-s}b_{s+1}^{(j-1)}$
(B7)

$Db_{s}^{(j)}=s\left(b_{s+1}^{(j-1)}-2\alpha
b_{s+1}^{(j)}+b_{s+1}^{(j+1)}\right)$ (B8)

$D^{n}b_{s}^{(j)}=s\left(D^{n-1}b_{s+1}^{(j-1)}-2\alpha
D^{n-1}b_{s+1}^{(j)}+D^{n-1}b_{s+1}^{(j+1)}-2(n-1)D^{n-2}b_{s+1}^{(j)}\right)$
(B9)

## References

  * Alessi & Sánchez (2015) Alessi, E. M., & Sánchez, J. P. 2015, Journal of Guidance, Control, and Dynamics, 39, Number 2, doi: 10.2514/1.g001237
  * Armellin et al. (2010) Armellin, R., Di Lizia, P., Berz, M., & Makino, K. 2010, Celestial Mechanics and Dynamical Astronomy, doi: 10.1007/s10569-010-9281-7
  * Benitez & Gallardo (2008) Benitez, F., & Gallardo, T. 2008, Celest Mech Dyn Astr, 101, 289, doi: 10.1007/s10569-008-9146-5
  * Brouwer & Clemence (1961) Brouwer, D., & Clemence, G. M. 1961, in Methods of Celestial Mechanics (Academic Press), 610, doi: 10.1016/b978-1-4832-0075-0.50013-7
  * Chambers (1999) Chambers, J. E. 1999, Monthly Notices of the Royal Astronomical Society, 304, 793, doi: 10.1046/j.1365-8711.1999.02379.x
  * Dones et al. (1999) Dones, L., Gladman, B., Melosh, H. J., et al. 1999, Icarus, 142, 509, doi: 10.1006/icar.1999.6220
  * Fang & Margot (2011) Fang, J., & Margot, J.-L. 2011, The Astronomical Journal, 143, 25
  * Fang & Margot (2012) Fang, J., & Margot, J. L. 2012, Astronomical Journal, doi: 10.1088/0004-6256/143/1/25
  * Farinella et al. (1994) Farinella, P., Froeschlé, C., Froeschlé, C., et al. 1994, Nature, 371, 314, doi: 10.1038/371314a0
  * Folkner et al. (2014) Folkner, W. M., Williams, J. G., Boggs, D. H., Park, R. S., & Kuchynka, P. 2014, Interplanet. Netw. Prog. Rep
  * Fuentes-Munoz & Scheeres (2020) Fuentes-Munoz, O., & Scheeres, D. J. 2020, in 2020 AAS/AIAA Astrodynamics Specialist Conference - Lake Tahoe
  * Giorgini & JPL Solar System Dynamics (2021) Giorgini, J., & JPL Solar System Dynamics. 2021, NASA/JPL Horizons On-Line Ephemeris System.  http://ssd.jpl.nasa.gov/?horizons
  * Gladman et al. (1997) Gladman, B. J., Migliorini, F., Morbidelli, A., et al. 1997, Science, 277, 197, doi: 10.1126/science.277.5323.197
  * Greenberg et al. (1988) Greenberg, R., Carusi, A., & Valsecchi, G. B. 1988, Icarus, doi: 10.1016/0019-1035(88)90125-X
  * Gronchi (2005) Gronchi, G. F. 2005, Celestial Mechanics and Dynamical Astronomy, doi: 10.1007/s10569-005-1623-5
  * Hedo et al. (2020) Hedo, J. M., Fantino, E., Ruíz, M., & Peláez, J. 2020, Astronomy & Astrophysics, doi: 10.1051/0004-6361/201936502
  * Hedo et al. (2018) Hedo, J. M., Ruíz, M., & Peláez, J. 2018, Monthly Notices of the Royal Astronomical Society, doi: 10.1093/mnras/sty1598
  * Heggie & Rasio (1996) Heggie, D. C., & Rasio, F. A. 1996, Monthly Notices of the Royal Astronomical Society, 282, 1064, doi: 10.1093/mnras/282.3.1064
  * Jones et al. (2018) Jones, R. L., Slater, C. T., Moeyens, J., et al. 2018, The Large Synoptic Survey Telescope as a Near-Earth Object discovery machine, doi: 10.1016/j.icarus.2017.11.033
  * JPL Solar System Dynamics & Center for NEO Studies (2021) (CNEOS) JPL Solar System Dynamics, & Center for NEO Studies (CNEOS). 2021, JPL Small-Body Database Search Engine
  * Kinoshita & Nakai (2007) Kinoshita, H., & Nakai, H. 2007, Celestial Mechanics and Dynamical Astronomy, 98, 67
  * Meyer & Scheeres (2021) Meyer, A. J., & Scheeres, D. J. 2021, Icarus, 367, 114554, doi: https://doi.org/10.1016/j.icarus.2021.114554
  * Meyer et al. (2021) Meyer, A. J., Scheeres, D. J., Naidu, S., et al. 2021, in AAS/Division of Dynamical Astronomy Meeting, Vol. 53, AAS/Division of Dynamical Astronomy Meeting, 405.06
  * Michel & Froeschlé (1997) Michel, P., & Froeschlé, C. 1997, Icarus, doi: 10.1006/icar.1997.5727
  * Michel et al. (1996) Michel, P., Froeschlé, C., & Farinella, P. 1996, in Worlds in Interaction: Small Bodies and Planets of the Solar System, ed. H. Rickman & M. J. Valtonen (Dordrecht: Springer Netherlands), 151–164
  * Michel et al. (1997) Michel, P., Froeschlé, C., & Farinella, P. 1997, Celestial Mechanics and Dynamical Astronomy, 69, 133, doi: 10.1007/978-94-017-1321-4-11
  * Michel et al. (2005) Michel, P., Morbidelli, A., & Bottke, W. F. 2005, Comptes Rendus Physique, 6, 291, doi: 10.1016/j.crhy.2004.12.013
  * Milani et al. (1989) Milani, A., Carpino, M., Hahn, G., & Nobili, A. M. 1989, Icarus, 78, 212, doi: 10.1016/0019-1035(89)90174-7
  * Milani et al. (2005) Milani, A., Chesley, S. R., Sansaturio, M. E., Tommei, G., & Valsecchi, G. B. 2005, Icarus, doi: 10.1016/j.icarus.2004.09.002
  * Milani & Knežević (1990) Milani, A., & Knežević, Z. 1990, Celestial Mechanics and Dynamical AStronomy, 49, 347, doi: 10.1007/BF00049444
  * Morbidelli et al. (2009) Morbidelli, A., Brasser, R., Tsiganis, K., Gomes, R., & Levison, H. F. 2009, Astronomy and Astrophysics, 507, 1053, doi: 10.1051/0004-6361/200912876
  * Murray & Dermott (2000) Murray, C. D., & Dermott, S. F. 2000, Solar System Dynamics (Cambridge University Press), doi: 10.1017/cbo9781139174817
  * Nesvorný & Bottke (2004) Nesvorný, D., & Bottke, W. F. 2004, Icarus, 170, 324, doi: 10.1016/J.ICARUS.2004.04.012
  * Öpik (1976) Öpik, E. J. 1976, Interplanetary encounters: close-range gravitational interactions, Developments in solar system- and space science (Elsevier Scientific Pub. Co.).  https://books.google.com/books?id=-HXvAAAAMAAJ
  * Pokorný & Vokrouhlický (2013) Pokorný, P., & Vokrouhlický, D. 2013, Icarus, doi: 10.1016/j.icarus.2013.06.015
  * Pravec et al. (2016) Pravec, P., Scheirich, P., Kušnirák, P., et al. 2016, Icarus, 267, 267
  * Rein & Spiegel (2014) Rein, H., & Spiegel, D. S. 2014, Monthly Notices of the Royal Astronomical Society, doi: 10.1093/mnras/stu2164
  * Roy (2004) Roy, A. 2004, Orbital Motion, Fourth Edition (Routledge), doi: 10.1201/9781420056884
  * Scheeres et al. (2020) Scheeres, D., McMahon, J., Bierhaus, E., et al. 2020, Bulletin of the AAS, 52\.  https://baas.aas.org/pub/2020n6i217p06
  * SSD/CNEOS API Service (2021) SSD/CNEOS API Service. 2021, SSD/CNEOS Small-Body DataBase API.  https://ssd.jpl.nasa.gov/tools/sbdb_lookup.html#/
  * Tancredi (1998) Tancredi, G. 1998, Celestial Mechanics and Dynamical Astronomy, 70, 181, doi: 10.1023/A:1008331422678
  * Valsecchi et al. (2015) Valsecchi, G. B., Alessi, E. M., & Rossi, A. 2015, Celestial Mechanics and Dynamical Astronomy, 123, 151, doi: 10.1007/s10569-015-9631-6
  * Valsecchi et al. (2003) Valsecchi, G. B., Milani, A., Gronchi, G. F., & Chesley, S. R. 2003, Astronomy & Astrophysics, 408, 1179, doi: 10.1051/0004-6361:20031039
  * Vokrouhlický et al. (2000) Vokrouhlický, D., Milani, A., & Chesley, S. R. 2000, Icarus, 148, 118, doi: 10.1006/ICAR.2000.6469
  * Vokrouhlický et al. (2012) Vokrouhlický, D., Pokorný, P., & Nesvorný, D. 2012, Icarus, doi: 10.1016/j.icarus.2012.02.021
  * Wisdom & Holman (1991) Wisdom, J., & Holman, M. 1991, The Astronomical Journal, 102, 1528, doi: 10.1086/115978
  * Wiśniowski & Rickman (2013) Wiśniowski, T., & Rickman, H. 2013, Acta Astronomica