OSSOS. XXIII. 2013 VZ₇₀ and the Temporary Coorbitals of the Giant Planets

Coorbital objects are found in the 1:1 mean-motion resonance with a planet. Resonance membership is determined by inspecting the evolution of the resonant angle ϕ₁₁ = λ − λ_P , where λ = Ω + ω + M is the mean longitude, P denotes the planet, Ω is the longitude of the ascending node, ω is the argument of pericenter, and M is the mean anomaly. The resonant angle ϕ₁₁ must librate rather than circulate (i.e., ϕ₁₁ must occupy a bounded range) in order for an object to be considered to be in coorbital resonance. Like other n:1 resonances, the 1:1 mean-motion resonance includes multiple libration islands; objects in these islands are called leading Trojans (mean 〈ϕ₁₁〉 = +60°), trailing Trojans (〈ϕ₁₁〉 = 300° = −60°), quasi-satellites (〈ϕ₁₁〉 = 0°), or horseshoe coorbitals (〈ϕ₁₁〉 = 180°). The motions of Trojans librate around one of the L4 or L5 Lagrangian points, while the paths of horseshoe coorbitals encompass all of the L3, L4, and L5 Lagrangian points; quasi-satellites appear to orbit the planet (while not actually being bound to it). Quasi-satellites and horseshoe coorbitals are almost always unstable and thus temporary (e.g., Mikkola et al. 2006; Ćuk et al. 2012; Jedicke et al. 2018), with the exception of Saturn's moons Epimetheus and Janus, which are horseshoe coorbitals of each other (Fountain & Larson 1978). Greenstreet et al. (2020) and Li et al. (2018) discuss the existence of high-inclination (i > 90°) objects temporarily trapped in a 1:−1 retrograde "coorbital" resonance with Jupiter and Saturn, respectively, although these are not coorbitals in the traditional sense described above; since they orbit the Sun in the opposite direction than the planet, retrograde coorbitals are not protected from close approaches with the planet the way that prograde coorbitals are, nor do the resonant island librations (i.e., Trojan, horseshoe, quasi-satellite motion) behave in the traditional sense in the retrograde configuration.

For many planets, the coorbital phase space is unstable owing to perturbations from neighboring planets (e.g., Nesvorný & Dones 2002; Dvorak et al. 2010). Innanen & Mikkola (1989) first suggested, at a time when only the Jovian Trojans were known, that populations of objects in stable 1:1 resonance with each of the other giant planets may exist; their analysis showed that the exact Lagrangian points are unstable for Saturn, but that Trojans farther from the resonance center (featuring larger libration amplitudes) could be stable for at least 10 Myr. These results were confirmed by Holman & Wisdom (1993). Using longer timescales than previous studies, de la Barre et al. (1996) specifically studied the stability of Saturnian Trojans and found that Saturnian Trojans could only be long-term (>428 Myr) stable with very specific conditions: very small eccentricity (<0.028), ϕ₁₁ libration amplitude greater than 80°, ω libration about a point 45° ahead of Saturn's ω, and constraints on the timing of the maximum eccentricity relative to the timing of Jupiter's maximum eccentricity, so that Jupiter and the Trojans do not approach close enough to dislodge the Trojan from Saturn's 1:1 resonance. Nesvorný & Dones (2002) showed that while Neptunian Trojans may have only been depleted by a factor of 2 over the age of the solar system, the Saturnian Trojans would have been depleted by a factor of 100. Studying the cause of the instability of Saturnian Trojans, Marzari & Scholl (2000) and Hou et al. (2014) found that the instability is caused by interactions between mean motion and secular resonances. Huang et al. (2019) investigated the stability of retrograde Saturnian coorbitals and found that they are always unstable owing to an overlap with the ν₅ and ν₆ secular resonances. Given these destabilizing factors, causing any primordial population to have been mostly depleted and allowing only small niches to be long-term stable, it is not surprising that no long-term stable Saturnian Trojans have been discovered to date.

Only Mars, Jupiter, and Neptune have known populations of long-term (>Gyr) stable Trojans (which thus might be primordial; Wolf 1906; Bowell et al. 1990; Levison et al. 1997; Marzari et al. 2003; Scholl et al. 2005). These long-term stable Trojan populations are important for understanding planet formation processes. As a few examples, Polishook et al. (2017) suggested that the Martian Trojans are likely to be impact ejecta from Mars, and they used the mass of the current Trojan cloud to constrain how much the orbit of Mars could have evolved during the phase of collisions. Morbidelli et al. (2005) showed that in order to reproduce the wide inclination distribution of the Jovian Trojans, the Trojans must have been captured from an excited disk during a migration phase rather than having formed in place together with Jupiter. Nesvorný et al. (2013) demonstrated that a sudden displacement of Jupiter's semimajor axis can explain the asymmetry seen between the L4 and L5 clouds and use the mass of the Jovian Trojan clouds to estimate the mass of the primordial planetesimal disk. Gomes & Nesvorný (2016) used the observed mass of Neptunian Trojans to infer that Neptune migrated slightly past its current location and then back, destabilizing the cloud, as we would otherwise observe a more massive cloud. Parker (2015) demonstrated that if Neptune's migration and eccentricity damping was fast, the disk that it migrated into and captured Trojans from must already have been dynamically excited prior to Neptune's arrival in order to reproduce our observed orbital distribution.

While only three planets are known to have long-term stable Trojans, scattering objects (scattering trans-Neptunian objects (TNOs), Centaurs, and even some objects originating in the asteroid belt ¹⁴ ) can become temporary coorbitals, transiently captured into unstable resonance (Alexandersen et al. 2013; Greenstreet et al. 2020). All solar system planets except Mercury, Mars, and Jupiter now have known populations of temporary coorbitals on prograde (i < 90°) orbits (Wiegert et al. 1998; Karlsson 2004; Mikkola et al. 2004; Horner & Lykawka 2012; Alexandersen et al. 2013; Greenstreet et al. 2020). Temporary "sticking" like this also occurs in other resonances (e.g., Duncan & Levison 1997; Tsiganis et al. 2000; Alvarez-Candal & Roig 2005; Lykawka & Mukai 2007; Volk et al. 2018; Yu et al. 2018). While long-term stable Trojans inform us of conditions in the time of planet formation and migration, the temporarily captured coorbitals inform us about properties of the scattering population. For example, Alexandersen et al. (2013) confirmed the Shankman et al. (2013) finding that the size distribution of the scattering population must have a transition in order to explain the observed ratio of small, nearby scattering objects (including Uranian coorbitals) to larger, more distant ones (including Neptunian coorbitals).

Horner & Wyn Evans (2006) integrated the Centaurs known at the time, demonstrating that Centaurs do indeed get captured into temporary coorbital resonance with the giant planets, claiming that Jupiter should have by far the most temporary coorbitals, followed by Saturn and hardly any for Uranus and Neptune. Alexandersen et al. (2013) pointed out that using the known Centaurs as the starting sample is biased toward having more objects nearer the Sun, and thus more captures for the inner giant planets, resulting in a disagreement with the sample of temporary coorbitals known at the time; they instead used a model that started with scattering TNOs that scatter inward to become Centaurs and temporary coorbitals, to demonstrate that a TNO origin can explain the distribution of the temporary coorbitals of Neptune and Uranus.

In this paper we describe the discovery of the first known Saturnian horseshoe coorbital, 2013 VZ₇₀, and demonstrate its temporary nature (Section 2). Furthermore, we expand on the analysis of Alexandersen et al. (2013) to analyze the populations of temporary coorbitals of all four giant planets, in an attempt to demonstrate the likely origin of 2013 VZ₇₀ and similar objects. We use numerical integrations to construct a steady-state distribution model of the scattering TNOs and temporary coorbitals of the giant planets (Section 3). Lastly, we use survey simulations, exposing our model to the survey biases of a well-understood set of surveys, in order to compare our theoretical predictions to real detections of this population (Section 4).

2013 VZ₇₀ was discovered by the Outer Solar System Origins Survey (OSSOS; Bannister et al. 2016, 2018) in images taken on 2013 November 1 using the MegaCam wide-field imager (Boulade et al. 2003) on the Canada–France–Hawaii Telescope (CFHT). The object was subsequently measured in 37 tracking observations from 2013 August 9 to 2018 January 18 (for the full list of astrometric measurements, see MPEC 2021-Q55; Bannister et al. 2021). With 4.5 yr of high-accuracy astrometry, the orbit is very well known, being a = 9.1838 ± 0.0002 au, e = 0.097145 ± 0.000011, i = 1204110 ± 000006, Ω = 21522021 ± 000008, ω = 245754 ± 0006, M = 291425 ± 0006 for epoch = JDT 2,456,514.0. Here a, e, i are the barycentric semimajor axis, eccentricity, and inclination, respectively; the uncertainties are calculated from a covariance matrix using the orbit-fitting software Find_Orb (Gray 2011) and the JPL DE430 planetary ephemerides (Folkner et al. 2014). This orbit is very close to that of Saturn, although the two bodies are separated by ∼180° on the sky. From dynamical integrations we found that 2013 VZ₇₀ is in fact in the 1:1 mean-motion resonance with Saturn, in a horseshoe configuration (see Figure 1). However, the best-fit clone only remains resonant for about 11 kyr before leaving the resonance and rejoining the scattering population.

Zoom In Zoom Out Reset image size

We investigated whether the best-fit orbit could be near a stability boundary by generating orbit clones from appropriate resampling of our astrometry, allowing us to test whether any orbit consistent with the astrometry featured long-term stability. Each clone was produced by resampling all the astrometry (using a normal distribution with standard deviation equal to the mean residual of the best fit, 0146) and fitting a new orbit. This process was repeated 10,000 times using Find_Orb. The distribution of orbits generated by this process explicitly shows how the uncertainties of some of the orbital parameters are strongly coupled, as can be seen in Figure 2. From these 10,000 clones, we identified the most extreme orbits (largest and smallest value of each parameter) and integrated these eight clones (labeled in Figure 2), as well as the best-fit orbit. These dynamical integrations were done using Rebound (Rein & Liu 2012) with the WHFast (Kinoshita et al. 1991; Wisdom & Holman 1991; Rein & Tamayo 2015) symplectic integrator. The eight major planets and Pluto ¹⁵ were included as massive perturbers, and an integration step size of 5% of Mercury's initial orbital period (≃0.012 yr ≃ 4.39 days) was used, while output was saved approximately three times per year. As can be seen in Figures 3 and 4, the future evolution of all of the clones involves an initial period in coorbital resonance, but all clones leave the resonance between 6 and 26 kyr from now. The large range of resonance exit times is due to the highly chaotic nature of the orbit. We used a second set of numerical integrations, where clones were displaced infinitesimally (10⁻¹³ to 10⁻¹² au, or 1.5–15 cm) relative to each of the above clones, to estimate the object's current timescale for chaotic divergence (the Lyapunov timescale); we found this to be 410 ± 60 yr. ¹⁶ 2013 VZ₇₀ is thus definitely in coorbital resonance now, but the chaotic nature of the orbit means that the duration of this temporary resonance capture will likely not be constrained further, even with additional observations.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We proceed to investigate the potential origin of temporary coorbitals like 2013 VZ₇₀. We model the source of the giant planet temporary coorbitals and investigate their detectability in characterized surveys. We produced a steady-state distribution of scattering objects in the a < 34 au region from orbital integrations similar to those used in Alexandersen et al. (2013); the details of those integrations can be found in the supplementary material of that paper. We primarily outline the deviations from those used in the previous paper below.

To perform the dynamical integrations, we used the N-body code SWIFT-RMVS4 (provided by Hal Levison, based on the original SWIFT; Levison & Duncan 1994) with a base time step of 25 days and an output interval of 50 yr for the orbital elements of the planets and any particle that at the moment had a < 34 au. The gravitational influences of the four giant planets and the Sun were included. The system starts with 8500 particles, derived from the 34 au < a < 200 au scattering portion of the Kaib et al. (2011) model of the outer solar system. Particles were removed from the simulation when they hit a planet, they went outside 2000 au or inside 2 au from the Sun (since they would either interact with the terrestrial planets that are absent in our simulations or would rapidly be removed from the solar system by Jupiter), or the final integration time of 1 Gyr was reached. Since 1 Gyr is substantially longer than the dynamical lifetime of Centaurs and scattering TNOs, we thoroughly sample the a < 34 au phase space, despite the limited number of initial particles. The output is combined along the time axis to produce a distribution of approximately 300 million sets of orbital elements. As in Alexandersen et al. (2013), we confirm that the distribution in the first 100 Myr is similar enough to the distribution in the following 900 Myr (because the a < 34 au region is populated very quickly despite starting off empty) that we can treat the distribution as a whole as being in steady state. We also ran a similar simulation with particles drawn from the modified version of the Kaib et al. (2011) model also used in Alexandersen et al. (2013) and Shankman et al. (2013); this modified version was generated with an assumption of the primordial planetesimal disk being more dynamically excited. As in Alexandersen et al. (2013), we find that using the standard Kaib et al. (2011) model and the modified version as our starting condition makes little quantitative difference on the end results. For the rest of this work we will therefore only be referring to the results from the simulations using the standard Kaib et al. (2011) model.

The method for determining coorbital behavior in the particle histories is also very similar to that used in Alexandersen et al. (2013), with some small modifications. To diagnose whether particles are coorbital, the orbital histories (at 50 yr output intervals) were scanned using a running window 30 kyr long for Uranus/Neptune and 5 kyr long for Jupiter/Saturn; this window size was chosen to be several times longer than the typical Trojan libration period at the given planet (∼1 kyr for Saturn as seen in Figure 4). A particle was classified as a coorbital if, within the running window, both its average semimajor axis was less than 0.2 au from the average semimajor axis of a given planet and no individual semimajor axis value deviated more than R_H from that of the planet. Here R_H is the planet's Hill sphere radius (Murray & Dermott 1999), where R_H = 0.35 au for Jupiter, R_H = 0.44 au for Saturn, R_H = 0.47 au for Uranus, and R_H = 0.77 au for Neptune. Further determination of which resonant island a coorbital is librating in was made identically to the method used in Alexandersen et al. (2013).

Our results are in good agreement with those for Uranus and Neptune in Alexandersen et al. (2013). Table 1 contains the fraction of the steady-state population in coorbital motion with each of the giant planets at any given time, as well as the distribution of coorbitals between horseshoe, Trojan, and quasi-satellite orbits. The coorbital fractions for Uranus and Neptune are slightly higher than in Alexandersen et al. (2013), despite very similar methodology; however, these results agree within their expected accuracy. The capture fraction decreases from Neptune through to Jupiter (with almost 4000 times fewer coorbitals than Neptune), although this is unsurpring given that the source of the scattering objects is beyond Neptune and that the dynamical timescales (orbital period, libration timescale) are longer farther from the Sun. An interesting result is that while Neptune's and Uranus's coorbitals are roughly equally distributed between horseshoe and Trojan coorbitals, Saturn seems to very preferentially capture scattering objects into horseshoe orbits, and Jupiter has a much larger fraction of quasi-satellites than any of the other planets. Table 1 also shows the mean, median, and maximum duration of a capture in coorbital resonance with the planets, as well as the mean, median, and maximum number of captures experienced by particles with at least one episode of coorbital motion with a given planet. The mean and median coorbital lifetimes and number of captures for Uranus and Neptune are also within a factor of two of those in Alexandersen et al. (2013), which we thus adopt as the uncertainty. Note that the captures into coorbital motion with Saturn and Jupiter are typically significantly shorter than for Uranus and Neptune, although if the lifetimes are represented in units of orbital periods rather than years, the Jovian captures actually have the second-longest lifetimes, after Uranus. While the different coorbital resonances no doubt experience different interactions with secular resonances and experience different perturbations from neighboring planets, it is noteworthy that the median numbers of orbital periods for coorbital captures for all four planets are within a factor of 2.5 of each other. However, going from Neptune to Jupiter, particles are increasingly unlikely to have multiple captures, presumably due to the increasing ability of the planet to scatter the objects to large semimajor axes; this results in particles on average spending both more total time and more total orbital periods in coorbital motion with Uranus and Neptune than Saturn and Jupiter.

In order to compare our dynamical model to the real detections, we run the model through the OSSOS Survey Simulator (Bannister et al. 2018; Lawler et al. 2018a). The survey simulator generates one object at a time (with an orbit drawn from the dynamical model and an H magnitude drawn from a parametric model discussed later) and assesses whether the object would have been discovered by the input surveys. We used all of the characterized surveys with sufficient characterization available for use in the simulator: the Canada–France Ecliptic Plane Survey (CFEPS; Petit et al. 2011), the CFEPS High-Latitude extension (HiLat; J.-M. Petit et al. 2017), the Alexandersen et al. (2016) survey, and OSSOS (Bannister et al. 2018). These surveys combined will be referred to as OSSOS++.

4.1. Orbital Distribution

For our survey simulations, it is preferable to have orbital distribution functions rather than an orbital distribution composed of a fixed number of discrete particles. This allows for the simulator to be run for as long as necessary, without producing duplicate identical particles. We have set up independent distributions for the coorbitals and the scattering objects, as described below, both inspired by the distribution seen in the integrations discussed in Section 3. Our model files and scripts for use with the OSSOS Survey Simulator (Lawler et al. 2018a) are provided at DOI:10.11570/21.0008 for anybody curious to use this model distribution.

4.1.1. Scattering Objects

We use the output from the Section 3 integrations, taking every particle's orbit at every time step and binning them using bin sizes of 0.5 au, 0.02, and 20 in a, e, and i space, respectively. The survey simulator reads this binned table, randomly selects a bin weighted by the number of particles that went into the bin, and then randomly assigns a, e, and i from a uniform distribution within the bin. Ω, ω, and M are all assigned randomly from a uniform distribution from 0° to 360°, since the orientations of scattering objects' orbits are random. This process allows us to draw essentially infinite unique particles that follow a distribution consistent with the steady-state distribution from Section 3.

4.1.2. Coorbital Objects

We cannot simply bin the coorbital distributions as we did for the scattering distribution. The numbers of coorbitals in the Section 3 integrations are low (particularly for Jupiter), and the number of dimensions we would need to bin is higher since the resonant angle ϕ₁₁ is also important for the coorbital distribution. Instead, we opted to use parametric distributions, fitted to the distributions seen in Section 3.

In this simplified parametric model, the semimajor axis of the coorbital is always set equal to that of the planet, since the few tenths of au variability do not influence detectability by sky surveys as much as the details of the eccentricity and inclination distribution. The eccentricity is modeled with a normal distribution, centered at 0 with a width w_e , multiplied by ${\sin }^{2}(e)$, truncated to [e_min, e_max] ¹⁷ :

$$\begin{eqnarray}&&f(e| {e}_{\min }\leqslant e\leqslant {e}_{\max })=\displaystyle \frac{{\sin }^{2}(e)}{{w}_{e}\sqrt{2\pi }}\exp \left(\displaystyle \frac{-{e}^{2}}{2{w}_{e}^{2}}\right).\end{eqnarray} \tag{ 1 }$$

This functional form has little physical motivation and was merely chosen, as it provides in the end a good fit to the distribution seen in our integrations. The inclination is modeled as a normal distribution, with center at 0° and a width w_i , multiplied by $\sin (i)$, truncated at i_max:

$$\begin{eqnarray}&&f(i| 0^\circ \leqslant i\leqslant {i}_{\max })=\displaystyle \frac{\sin (i)}{{w}_{i}\sqrt{2\pi }}\exp \left(\displaystyle \frac{-{i}^{2}}{2{w}_{i}^{2}}\right).\end{eqnarray} \tag{ 2 }$$

This is simply a Normal distribution modified to account for the spherical coordinate system. Lastly, Ω and M are chosen randomly from a uniform distribution [0°, 360°), while ω is calculated from ϕ₁₁, the value of which depends on the type of coorbital. The different types of coorbitals are generated using the ratios in Table 1. The details of the selection of a ϕ₁₁ value are similar to those used in Alexandersen et al. (2013), accounting for a distribution of libration amplitudes and the fact that the center of libration is offset away from 60°/300° for Trojans with large libration amplitudes. The values of w_e , w_i , e_min, e_max, and i_max used in this work for each planet's coorbital population are shown in Table 2.

Table 2. Orbital Parameters Used for Each Planet in the Parametric Model Described in Section 4.1.2

Planet	we	emin	emax	wi (deg)	imax (deg)
Jupiter	0.188	0.0	0.523	16.1	40.6
Saturn	0.127	0.0	0.707	15.6	89.6
Uranus	0.134	0.0	0.998	19.4	51.3
Neptune	0.123	0.0	0.974	18.1	80.9

Download table as:  ASCIITypeset image

4.2. Absolute Magnitude Distributions

4.3. Population Estimate

4.4. Expected versus Detected Numbers

4.5. Caveat

2013 VZ₇₀ is the first known temporary Saturnian horseshoe coorbital, remaining resonant for 6–26 kyr; it likely originates in the trans-Neptunian region. Our simulations show that all the giant planets should have temporary coorbitals of TNO origin, although Jupiter has approximately a factor of 4000 fewer than Neptune; the durations of the coorbital captures are significantly more short-lived for Saturn and Jupiter than for Uranus and Neptune. Our simulations show that the Neptunian and Uranian coorbitals should be roughly equally distributed between horseshoe and Trojan coorbitals, Saturn very preferentially captures scattering objects into horseshoe orbits, and Jupiter should have a much larger fraction of its temporary coorbitals be quasi-satellites than any of the other planets. Accounting for observing biases in a set of well-characterized surveys (CFEPS, Petit et al. 2011; HiLat, J.-M. Petit et al. 2017; the Alexandersen et al. 2016 survey; and OSSOS, Bannister et al. 2018), we find that the fraction of a < 34 au scattering objects that are in temporary coorbital motion is higher in the real observations (13.7%) than in simulated observations (2.9%). However, for the distribution of the temporary coorbitals among the giant planets, we find that our predictions (∼0:1:2:2 for J:S:U:N) are consistent with the observations (0:1:1:2).

This work is based on TNO discoveries obtained with MegaPrime/MegaCam, a joint project of the Canada–France–Hawaii Telescope (CFHT) and CEA/DAPNIA, at CFHT, which is operated by the National Research Council (NRC) of Canada, the Institute National des Sciences de l'Universe of the Centre National de la Recherche Scientifique (CNRS) of France, and the University of Hawaii. A portion of the access to the CFHT was made possible by the Academia Sinica Institute of Astronomy and Astrophysics, Taiwan. This research used the facilities of the Canadian Astronomy Data Centre operated by the National Research Council of Canada with the support of the Canadian Space Agency.

The authors wish to recognize and acknowledge the very significant cultural role and reverence that the summit of Maunakea has always had within the indigenous Hawaiian community. We are most fortunate to have the opportunity to conduct observations from this mountain. We would also like to acknowledge the maintenance, cleaning, administrative, and support staff at academic and telescope facilities, whose labor maintains the spaces where astrophysical inquiry can flourish.

The authors thank Hanno Rein for useful discussions and help regarding how to estimate the Lyapunov timescale through numerical integrations. The authors also thank Jeremy Wood for the suggestion of looking at the capture durations in terms of orbital periods rather than only absolute time.

S.G. acknowledges support from the Asteroid Institute, a program of B612, 20 Sunnyside Ave., Suite 427, Mill Valley, CA 94941. Major funding for the Asteroid Institute was generously provided by the W. K. Bowes Jr. Foundation and Steve Jurvetson. Research support is also provided from Founding and Asteroid Circle members K. Algeri-Wong, B. Anders, R. Armstrong, G. Baehr, The Barringer Crater Company, B. Burton, D. Carlson, S. Cerf, V. Cerf, Y. Chapman, J. Chervenak, D. Corrigan, E. Corrigan, A. Denton, E. Dyson, A. Eustace, S. Galitsky, L. & A. Fritz, E. Gillum, L. Girand, Glaser Progress Foundation, D. Glasgow, A. Gleckler, J. Grimm, S. Grimm, G. Gruener, V. K. Hsu & Sons Foundation Ltd., J. Huang, J. D. Jameson, J. Jameson, M. Jonsson Family Foundation, D. Kaiser, K. Kelley, S. Krausz, V. Laas, J. Leszczenski, D. Liddle, S. Mak, G. McAdoo, S. McGregor, J. Mercer, M. Mullenweg, D. Murphy, P. Norvig, S. Pishevar, R. Quindlen, N. Ramsey, P. Rawls Family Fund, R. Rothrock, E. Sahakian, R. Schweickart, A. Slater, Tito's Handmade Vodka, T. Trueman, F. B. Vaughn, R. C. Vaughn, B. Wheeler, Y. Wong, M. Wyndowe, and nine anonymous donors. S.G. acknowledges the support from the University of Washington College of Arts and Sciences, Department of Astronomy, and the DIRAC Institute. The DIRAC Institute is supported through generous gifts from the Charles and Lisa Simonyi Fund for Arts and Sciences and the Washington Research Foundation. This work was supported in part by NASA NEOO grant NNX14AM98G to LCOGT/Las Cumbres Observatory.

K.V. acknowledges support from NASA (grants NNX15AH59G and 80NSSC19K0785) and NSF (grant AST-1824869).

This work was supported by the Programme National de Plantologie (PNP) of CNRS-INSU co-funded by CNES.

This research has made use of NASA's Astrophysics Data System Bibliographic Services. This research made use of SciPy (Virtanen et al. 2020), NumPy (van der Walt et al. 2011), and matplotlib (a Python library for publication quality graphics Hunter 2007).

Facility: CFHT(MegaCam). -

Software: Python (van Rossum & de Boer 1991), Matplotlib (Hunter 2007), NumPy (van der Walt et al. 2011), SciPy (Virtanen et al. 2020), Find_Orb (Gray 2011), Rebound (Rein & Liu 2012), fit_radec & abg_to_aei (Bernstein & Khushalani 2000) SWIFT-RMVS4 (Levison & Duncan 1994), Jupyter notebook (Kluyver et al. 2016).