https://orcid.org/0000-0002-0246-3303
Abstract
Asteroid (514107) 2015 BZ509 was discovered recently in Jupiter’s co-orbital region with a retrograde motion around the Sun. The known chaotic dynamics of the outer Solar system have so far precluded the identification of its origin. Here, we perform a high-resolution statistical search for stable orbits and show that asteroid (514107) 2015 BZ509 has been in its current orbital state since the formation of the Solar system. This result indicates that (514107) 2015 BZ509 was captured from the interstellar medium 4.5 billion years in the past as planet formation models cannot produce such a primordial large-inclination orbit with the planets on nearly coplanar orbits interacting with a coplanar debris disc that must produce the low-inclination small-body reservoirs of the Solar system such as the asteroid and Kuiper belts. This result also implies that more extrasolar asteroids are currently present in the Solar system on nearly polar orbits.
1 INTRODUCTION
Centaurs, the asteroids that roam the space between the giant planets of the Solar system, have a chaotic dynamical evolution governed for some by the close encounters with the giant planets, and for others by their successive hopping in and out of the outer planets web of mean motion resonances. The mean lifetime of the first group ranges from 1 to 10 Myr whereas that of the second group extends to 100 Myr as their resonant status provides them with some dynamical protection (Tiscareno & Malhotra 2003; Bailey & Malhotra 2009; Volk & Malhotra 2013). Capture in resonance may occur for prograde or retrograde motion with a greater likelihood for the latter making the retrograde resonant Centaurs possibly the oldest asteroid residents of the outer Solar system (Namouni & Morais 2015, 2017a). A number of Centaurs are currently known to be in retrograde resonance between the outer planets’ orbits (Morais & Namouni 2013b) but it is the discovery of asteroid (514107) 2015 BZ509, inside Jupiter’s co-orbital region, with an orbit of moderate eccentricity of 0.38 and a retrograde inclination of 163° (Wiegert, Connors & Veillet 2017) that has produced so far the most puzzling example of retrograde resonance in the Solar system.
Jupiter’s co-orbital region is known to host the Trojan asteroids that were captured mostly permanently by the planet during the late stage of Solar system formation (Tsiganis, Varvoglis & Dvorak 2005; Robutel & Gabern 2006; Nesvorný, Vokrouhlický & Morbidelli 2013; Jewitt 2018). Asteroid 2015 BZ509 shares the co-orbital region with the Trojans but moves in the opposite orbital direction. It sits on Jupiter’s peak of capture probability (Namouni & Morais 2017c) and would have an indefinitely stable orbit if the Solar system contained only Jupiter (Morais & Namouni 2013a, 2016). It is thought to be trapped temporarily as standard planet formation models do not produce Centaurs in situ with long-term stable retrograde orbits. In this framework, 2015 BZ509 could originate from distant reservoirs such as the scattered disc or the Oort cloud like typical Centaurs (Brasser et al. 2012). To reach Jupiter’s orbit, 2015 BZ509 would have had to cross the giant planets’ space whose chaotic dynamics and its 100 Myr lifetime time-scale suggest a recent capture at Jupiter’s orbit and preclude the identification of the asteroid’s origin during the final stage of planet formation some 4 billion years ago.
However, when the motion of 100 asteroid clones with orbits that differ slightly from that of 2015 BZ509 by amounts compatible with the orbital elements’ error bars, was simulated over 1 Myr, it was found to have a stable evolution (Wiegert et al. 2017). This time-scale is about two orders of magnitude longer than those of temporarily captured retrograde resonant asteroids (Morais & Namouni 2013b) hinting to a different origin from that of most Centaurs.
In this Letter, we report on a high-resolution statistical search for stable orbits for 2015 BZ509 that allows us to trace the asteroid’s origin back to the epoch of Solar system formation. In Section 2, we describe our numerical approach of simulating the evolution of one million clones of 2015 BZ509. In Section 3, we present our results that show 2015 BZ509 has been a co-orbital of Jupiter since the end of planet formation and has a strongly stable orbit that can live theoretically at least 43 billion years. In Section 4, we explain that this finding implies that 2015 BZ509 was captured from the interstellar medium and that there should be more extrasolar asteroids bound to the Solar system on nearly polar orbits.
2 ONE MILLION CLONE SIMULATION
The motion of the asteroid and the planets constitutes an N-body dynamical system whose phase space structure is complex and chaotic but does not preclude the presence of long-term stability islands even at or between the outer giant planets. Probing phase space to identify such islands requires a significant number of initial conditions especially since 2015 BZ509’s precise orbit is unknown and only a representation thereof in parameter space by a nominal orbit and a covariance matrix exists (Knežević & Milani 2012). We therefore simulated the evolution of one million asteroid clones that interact with the giant planets and the Galactic tide back to 43 billion years in the past to search for the most stable orbits. The existence of stable orbits over the age of the Solar system is by no means guaranteed especially as we probe one of the most chaotic regions in the Solar system. However if such orbits exist then they are the ones that correspond to the actual motion of the asteroid and not the short-lived unstable orbits. Choosing the former orbits over the latter is motivated by the copernican principle that 2015 BZ509 is not being observed at a preferred epoch in Solar system history.
The nominal orbit of 2015 BZ509 (Table 1) and its equinoctial covariance matrix were obtained from the AstDys data base1 for the Julian date 2457 800.5. The orbital elements of the planets were obtained from NASA JPL’s Horizons ephemeris system2 for the same epoch. Clone orbits were generated using the Cholesky method for multivariate normal distributions (Thomopoulos 2013). This method consists in writing the equinoctial covariance matrix as C = LLt where L is a lower triangular matrix. A sample of one million clones is generated from the equinoctial nominal elements e0i where 1 ≤ i ≤ 6 as ei = e0i + rjLij; summation over the repeated j index is implied with j ≤ i, and rj is a six-dimensional vector with components generated independently from a normal distribution with mean 0 and variance 1. The one million sample achieves 10−7 to 10−5 relative error in reproducing the observational covariance matrix whereas smaller samples achieve larger relative error (e.g. 10−2 to 10−1 for 1000 clones) because of the large dimension of the initial conditions’ space.
Nominal orbital elements of asteroid (514107) 2015 BZ509 and their standard deviations.
| Julian date | 2457 800.5 |
|---|---|
| Semimajor axis (au) | 5.140 344 420 ± 5.7871 × 10−5 |
| Eccentricity | 0.3806 619 168 ± 6.2258 × 10−6 |
| Inclination (°) | 163.004 558 786 ± 2.0009 × 10−5 |
| Longitude of ascending node (°) | 307.376 837 106 ± 3.5503 × 10−5 |
| Argument of perihelion (°) | 257.44 919 796 ± 2.7616 × 10−4 |
| Mean anomaly (°) | 32.46 169 240 ± 4.4177 × 10−4 |
| Julian date | 2457 800.5 |
|---|---|
| Semimajor axis (au) | 5.140 344 420 ± 5.7871 × 10−5 |
| Eccentricity | 0.3806 619 168 ± 6.2258 × 10−6 |
| Inclination (°) | 163.004 558 786 ± 2.0009 × 10−5 |
| Longitude of ascending node (°) | 307.376 837 106 ± 3.5503 × 10−5 |
| Argument of perihelion (°) | 257.44 919 796 ± 2.7616 × 10−4 |
| Mean anomaly (°) | 32.46 169 240 ± 4.4177 × 10−4 |
Nominal orbital elements of asteroid (514107) 2015 BZ509 and their standard deviations.
| Julian date | 2457 800.5 |
|---|---|
| Semimajor axis (au) | 5.140 344 420 ± 5.7871 × 10−5 |
| Eccentricity | 0.3806 619 168 ± 6.2258 × 10−6 |
| Inclination (°) | 163.004 558 786 ± 2.0009 × 10−5 |
| Longitude of ascending node (°) | 307.376 837 106 ± 3.5503 × 10−5 |
| Argument of perihelion (°) | 257.44 919 796 ± 2.7616 × 10−4 |
| Mean anomaly (°) | 32.46 169 240 ± 4.4177 × 10−4 |
| Julian date | 2457 800.5 |
|---|---|
| Semimajor axis (au) | 5.140 344 420 ± 5.7871 × 10−5 |
| Eccentricity | 0.3806 619 168 ± 6.2258 × 10−6 |
| Inclination (°) | 163.004 558 786 ± 2.0009 × 10−5 |
| Longitude of ascending node (°) | 307.376 837 106 ± 3.5503 × 10−5 |
| Argument of perihelion (°) | 257.44 919 796 ± 2.7616 × 10−4 |
| Mean anomaly (°) | 32.46 169 240 ± 4.4177 × 10−4 |
The evolution of an asteroid clone back in time was followed in the system composed of the four giant planets and the Sun whose mass was augmented by those of the inner Solar system’s planets. The full three-dimensional Galactic tide (Heisler & Tremaine 1986) and relative inclination of the ecliptic and Galactic planes were taken into account. The Oort constants (A = 15.3 km s−1 kpc−1, B = −11.9 km s−1 kpc−1) and star density in the solar neighbourhood (ρ0 = 0.119 M⊙ pc−3) were taken from the recent Gaia DR1 determinations (Bovy 2017; Widmark & Monari 2017). The five-body problem with the Galactic tide is adequate to analyse the stability of 2015 BZ509 as stable clones surviving in the co-orbital region near the nominal orbit have a perihelion at 3.6 au far outside Mars’s orbit thus precluding close encounters with the inner Solar system’s planets. Numerical integration was carried out using the Bulirsch and Stoer algorithm with an error tolerance of 10−11. More standard symplectic-based alternatives such as the hybrid Mixed-Variable Symplectic integrator (MVS) (Chambers 1999) were not used because of the peculiar geometry of the retrograde co-orbital resonance that implies the asteroid encounters the planet twice per period. Tests with the hybrid MVS show that code switches to the Bulirsch and Stoer algorithm on a large portion of the orbit because of that geometry. To avoid such systematic and frequent algorithm switching, we opted for the Bulirsch and Stoer algorithm as it is also adequate to model large eccentricity orbit evolution (Wiegert & Tremaine 1999). Orbital evolution was monitored for the following events: collision with the Sun, collision with the planets, ejection from the Solar system, and reaching the inner 1 au semimajor axis boundary. No event at the inner boundary was registered.
3 LONG-TERM STABLE ORBITS
The simulation shows that the clone minimum and median lifetimes are respectively 0.29 and 6.48 Myr. Unstable clones that exit the co-orbital region undergo close encounters with the planets and are temporarily captured in mean motion resonances much like the general behaviour of Centaurs. They also follow a distinct path towards polar inclinations into a dynamical structure that extends to the Oort cloud and that we term ‘the polar corridor’, before being removed from the system. Fig. 1 illustrates this structure with a system snapshot after an evolution back in time of 40 Myr when 43 628 clones are present. The dynamical structure starting from the current location of 2015 BZ509 is centred around the curve of the asteroid’s Tisserand relation with Jupiter. The clones follow this curve in the first few million years of evolution only to be dispersed around it by their encounters with the other giant planets. The location breakdown at 40 Myr is as follows: 35 372 in retrograde co-orbital resonance with Jupiter, 52 in the inner Solar system, 3361 between Jupiter’s and Neptune’s orbits, 4296 trans-Neptunians with semimajor axes smaller than 1000 au, and 547 with semimajor axes larger than 1000 au that extend to the Oort cloud.
Distribution of 2015 BZ509’s clones at 40 Myr in the past in the semimajor axis-inclination plane (a) and the semimajor axis-eccentricity plane (b). The curve in (a) is given by the asteroid’s Tisserand relation with Jupiter. The four V-shaped curves in (b) denote the intersection at aphelion and perihelion of the clone’s orbit with those of Jupiter, Saturn, Uranus, and Neptune respectively from left to right.
At the 100 Myr Centaur maximum instability time-scale, 6577 clones remain stable, 75 per cent of which are sheltered by Jupiter’s co-orbital resonance. At the end of planet formation, 4.5 billion years in the past, most clones are lost: 553 811 increased their eccentricities to unity thereby reaching the Sun’s surface, 445 678 were ejected from the Solar system and 465 collided with a planet. The surviving clones number 46 of which 27 are in co-orbital resonance with Jupiter whereas the remaining 19 are dispersed between the current locations of the scattered disc and the inner Oort cloud as shown in Fig. 2. The 10 clones in the scattered disc region reside in the polar corridor and have high-inclination prograde orbits with perihelia at Uranus’s or Neptune’s orbits. Most of the nine clones in the Oort cloud region had their orbits extracted from the polar corridor by the Galactic tide and do not suffer close encounters with the planets at the corresponding epoch.
Distribution of 2015 BZ509’s clones at 4.5 billion years in the past in the semimajor axis-inclination plane (a) and the semimajor axis–eccentricity plane (b). The intersections at perihelion with the orbits of Uranus and Neptune are shown by the left and right curves, respectively.
The clustering of 60 per cent of long-term stable clones at 4.5 billion years in Jupiter’s co-orbital region shows that 2015 BZ509 has been in its current state since planet formation. For the clones in the current regions of the Oort cloud and the scattered disc, all semimajor axes are widely spaced and have each a 2 per cent chance of being the original semimajor axis.
The 27 long-term stable co-orbital clones share a number of orbital properties regarding their final states at 4.5 billion years illustrated by the example shown in Fig. 3. First, all but one clone have increased their semimajor axis above Jupiter’s. The average orbital elements and their standard deviations are semimajor axis 5.3319 ± 0.0531 au, eccentricity 0.2888 ± 0.0294, and inclination 161.9939 ± 3.0144°. Secondly, none of the clones librate with the 1:1 resonance arguments ϕ = λ − λJupiter or ϕ⋆ = λ − λJupiter − 2ω, where λ, λJupiter, and ω are respectively the mean longitudes of the clone and Jupiter and the clone’s argument of perihelion (Morais & Namouni 2013a). Thirdly, they are all solidly trapped in the Kozai–Lidov secular resonance with ω = 0° or 180°3 (Kozai 1962; Lidov 1962; Morais & Namouni 2016).
Evolution of a stable co-orbital clone’s orbital elements near 4.5 billion years in the past. The location of the semimajor axis is shown with respect to Jupiter’s co-orbital region (dashed lines). Inclination is normalized to 180° and corresponds to the top curve in the corresponding panel whereas eccentricity is the bottom curve. The possible resonant arguments ω, ϕ, and ϕ⋆ are defined in Section 3.
The phase space structure near 2015 BZ509’s orbit may be visualized through the distribution of the 61 clones’ initial conditions that were present in the first billion years after planet formation (i.e. in the interval [−4.5:−3.5] billion years).4 At 3.5 billion years in the past, there were 36 clones in the co-orbital region, 15 in the scattered disc region, and 10 in the Oort cloud region. The distribution is shown in Fig. 4 with respect to the nominal orbit and the confidence levels of the covariance matrix according to the clone’s state at 4.5 billion years in the past. The clones do not cluster in terms of co-orbital orbits, non-co-orbital ones, ejected or Sun-colliding orbits indicating, in mathematical terms, that the set of long-term stable orbits, and the set of unstable orbits are dense in parameter space around the nominal orbit. In physical terms, this means that a stability island does exist but different nearby orbits within it have different lifetimes. Dynamical instability and removal from the co-orbital region is likely caused by a slow chaotic diffusion from the secular and near-mean motion resonances similar to that of the Trojan asteroids (Tsiganis et al. 2005; Robutel & Gabern 2006).
Distribution of the 61 initial conditions of 2015 BZ509s clones that were present from 3.5 to 4.5 billion years in the past. Panel (a) shows the initial conditions in the (esin ϖ, ecos ϖ) plane where e is the eccentricity and ϖ is the longitude of perihelion. Panel (b) shows the initial conditions in the (tan (I/2)sin Ω, tan (I/2)cos ϖ) plane where I is the inclination and Ω is the longitude of the ascending node. Panel (c) shows the initial conditions in the (a, λ) plane, where a is the semimajor axis and λ is the mean longitude. In each panel the 1σ, 2σ, and 3σ confidence levels of the covariance matrix are shown around the nominal orbit (large central cross). A clone’s state is shown at 4.5 billion years in the past: co-orbital clones with blue filled circles, non-co-orbital clones with empty green circles, ejected clones with small red crosses, and the only Sun-colliding clone with a medium orange cross. The location of the longest lived clone (42.91 billion years) is denoted by a large black filled circle.
Running the simulation further back in time to test orbital stability shows that the clones in the co-orbital region, scattered disc and Oort cloud number (14, 3, 8) at 9 billion years, (8, 3, 5) at 14 billion years, and (2, 2, 1) at 30 billion years. At 42.91 billion years, the last co-orbital clone leaves Jupiter’s orbit with a stable motion into the scattered disc region. The slow number decay of stable co-orbitals is further indication of possible chaotic diffusion. However, on such very long integration timespans, accumulation of numerical error could contribute to the clone’s earlier exit from the co-orbital region. That is why the 42.91 billion year estimate is likely a lower bound on the clone’s lifetime. A finer analysis of the dynamics in the co-orbital region is required to ascertain the various diffusion times associated with the secular and near-mean motion resonances as well as the role of numerical error in the asteroid’s evolution on very long time-scales. The decreasing number of Oort cloud clones is caused by the Galactic tide as it conserves the vertical component of angular momentum forcing slow cyclic oscillations of the clone’s eccentricity and inclination that lead to collisions with the planets (Heisler & Tremaine 1986).
4 INTERSTELLAR ORIGIN
The presence of 2015 BZ509 in retrograde co-orbital resonance early in the Solar system’s timeline is unexpected from the standpoint of Solar system formation theory as retrograde Centaurs are believed to originate in the scattered disc or the Oort cloud well after the planets settled down dynamically (Brasser et al. 2012). Furthermore, planet formation models cannot produce such a primordial large inclination orbit as that of 2015 BZ509 with the planets on nearly coplanar orbits interacting with a coplanar debris disc that must produce the low inclination small body reservoirs of the Solar system such as the asteroid and Kuiper belts (Pfalzner et al. 2015). This implies that 2015 BZ509 was captured from the interstellar medium. In this respect, even the low-probability long-lived orbits found in the scattered disc and Oort cloud regions (Fig. 2) should have an interstellar origin because of their location and high inclinations at the end of planet formation. Interstellar capture events can occur during planet formation in a tightly packed star cluster whose relaxation was more violent than the one that was thought to have formed the Oort cloud in the Sun’s birth cluster early in the Solar system’s history (Levison et al. 2010) as the smallest captured semimajor axis was only 1100 au.
The one million clone simulation provides further evidence that there are currently more extrasolar asteroids in the Solar system. In effect, if more objects were captured along 2015 BZ509 by Jupiter early in the Solar system’s history, the less stable orbits must have left the co-orbital region by way of chaotic diffusion into the polar corridor. This occurs because the N-body problem is time-reversible and unstable clones of 2015 BZ509 that are followed into the future exit the co-orbital region and end up in the polar corridor. The prominent presence of the polar corridor in the simulation over the age of the Solar System mainly in the trans-Neptunian region implies that it is currently populated by extrasolar asteroids. Interestingly, a structure similar to the polar corridor was observed in the known nearly polar trans-Neptunian objects (TNOs) (Gladman et al. 2009; Chen et al. 2016; Morais & Namouni 2017). Integrations of a 1000 clones of TNOs (471325) and 2008 KV42 over 1 billion years have shown their orbits to evolve towards larger semimajor axes while they cluster around 90° inclination (see fig. 1 in Chen et al. 2016).
The presence of extrasolar asteroids bound to the Solar system early in its history implies the need for a revision of planet formation theory as such interstellar contamination of small body reservoirs will affect not only the dynamics of small bodies but also their physical properties. Observed discrepancies such as that of Trojan colours may originate in a different Trojan origin at different planets (Jewitt 2018).
The stability search method presented here is new and is aimed at beating dynamical unpredictability in one of the most ch aotic regions of the Solar system using large statistics and intensive computing. It is similar in principle to the orbit determination of newly discovered multiplanet systems that are systematically vetted for stability upon discovery and only the stable orbits are chosen from the available range of observational error bars. In both problems, a stable configuration is preferable to an unstable one as the opposite would imply that the system is being observed at a preferred epoch. Applying systematically our new method to Centaurs and TNOs will help constrain their origin and improve our understanding of Solar system formation.
Footnotes
http://Hamilton.dm.unipi.it/astdys (Knežević & Milani 2012).
http://ssd.jpl.nasa.gov (Giorgini et al. 1996).
The values ω = 0° and 180° represent the same equilibrium as the Kozai–Lidov resonance is 180° periodic. Note that the eccentricity and inclination oscillate with the same phase because motion is retrograde. In the Kozai–Lidov cycle, eccentricity is maximal when motion is nearly co-planar. For retrograde motion this occurs when inclination is near 180° (Namouni & Morais 2017b).
We do not show similar data for a more recent interval to avoid an overcrowded plot.
ACKNOWLEDGEMENTS
We thank the anonymous reviewer for helpful comments. The numerical simulation was done at the Mésocentre SIGAMM hosted at the Observatoire de la Côte d’Azur. MHMM was supported by grant 2015/17962-5 of São Paulo Research Foundation (FAPESP).
REFERENCES
