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

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.


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 . 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 arXiv:2207.03527v1 [astro-ph.EP] 7 Jul 2022 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 . In this paper we study the orbital history of the targets of the exploration mission Janus : 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.

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.

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 ), 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.

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 (CNEOS) 2021). 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.

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. (2003Valsecchi et al. ( , 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.

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: 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: The separation is convenient to expand in the ratio of semi-major axes α jk as well as sines and cosines of { j , Ω j , λ j , k , Ω k , λ k , }. The ratio of semi-major axes is α jk = a j /a k if the perturber is external, or α jk = a j /a k if the perturber is internal. All the terms that depend on the longitudes {λ, λ j } 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: Whereᾱ jk = a j /a k if the perturber is external, orᾱ jk = 1 if the perturber is internal. The coefficients A jj , A jk , B jj , B jk are: 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, the potential becomes: Our complete set of states includes the mean anomaly at epoch σ j and L j = GM a j . The equations of motion become:ṗ 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 σ 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: 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 β i , γ i given a set of initial conditions. In order to solve for e ji , I ji , β i , γ i we proceed as follows. From the matrices of normalized eigenvectorsē ji ,Ī ji and the initial conditions h, k, p, q we form: These are four linear systems of equations, where S i , T i are the scaling factors of each eigenvector. Solving for the combined factors we can reconstruct the vectors e ji , I ji and the phase angles β i , γ i . 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. Table 1. Initial conditions of the Solar System propagation in Figure 3. Note-From ephemeris DE431 at Epoch: JD0 = 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.
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 σ j . Short term applications benefit from the improved characterization of the position of the bodies in their orbits. Solving for σ 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 σ j becomes: and the solution depends on whether the perturber is external or internal: wherec jk = a j /a 2 k in the case of an external perturber andc jk = −1/a j if the perturber is internal. The solution of the equation is simply a constant drift given by the rateσ. This element completes the set of elements of the secular model.

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: 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: where the constant κ is found as the ratio between Laplace coefficients: 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: with constants of integration: 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 (t), Ω(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) = h 2 1 (t) + k 2 1 (t) and inclination i 1 (t) = p 2 1 (t) + q 2 1 (t). The maximum and minimum values of eccentricity and inclination are: 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 Ω, ω 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 Ω, . In addition, given the relationship between the frequencies g and f , the relationship between the arguments rates isω = −2Ω. . 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 Ω,ω. 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: 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ω and˙ . 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.

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.

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. (2003Valsecchi et al. ( , 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 γ at the time of closest approach: 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.

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:Ḋ where D is the set of elements of choice and L(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 K = [a, e, i, Ω, ω, σ] and Cartesian coordinates. From the orbital elements representations available to choose, the current implementation uses the Delaunay elements: 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: 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: Figure 5. Logarithm of the errors in the computation of the final Keplerian elements given by pseudo-Öpik theory (PÖpikfirst 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.
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: 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 ∞ , 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.

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.8au, 0.2, 10 deg and with standard deviations of 0.05au, 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.

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.

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.. 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.

Short-term propagation of near-Earth objects using different models
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.

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. (2018Hedo et al. ( , 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.
The distributions of closest approach distances are shown in figure 8 and the unperturbed relative velocity V ∞ 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.

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 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 Ω − ω 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 . 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.
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 σ ∝ √ 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 ω, Ω. In the duration of our simulations of 500,000 most cases approach this uniformization as we showed in figure 11.
We can compute more rigorously whether the distributions can be considered uniform by conducting the chi-squared test on the longitude of perihelion = ω+Ω. 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 is true. Cases 1, 3, 4, (175706) 1996 FG3 and (35107)   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 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.

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 Ω and ω 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.

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 . 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. 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.  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.

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 semimajor 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. 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.

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;, 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 . 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  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.
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 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∞ (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. 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∞ (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. possible relative velocities is larger in all the planets. In addition, many encounters are found with quite slower V ∞ , 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.
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.

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.
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, Ω, ω, σ}. 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.

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: A jk e j e k cos ( j − k ) + B jk I j I k cos (Ω j − Ω k )    The coefficients A jk ,B jk ,B jj and A jj are: 3/2 (α jk ) (B1) 3/2 (α jk ) (B2) The definition of Laplace coefficient is: 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 d dα . b