# Dynamical Analysis of Three Distant Trans-Neptunian Objects with Similar Orbits

T. Khain, J.C. Becker, F.C. Adams, D.W. Gerdes, S.J. Hamilton, K. Franson, L. Zullo, M. Sako, K. Napier, Hsing Wen Lin, et al.

► To cite this version:

T. Khain, J.C. Becker, F.C. Adams, D.W. Gerdes, S.J. Hamilton, et al. Dynamical Analysis of Three Distant Trans-Neptunian Objects with Similar Orbits. Astron.J., 2018, 156 (6), pp.273. 10.3847/1538-3881/aaeb2a. hal-02066940

HAL Id: hal-02066940

https://hal.science/hal-02066940v1

Submitted on 29 Sep 2025

**HAL** is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire **HAL**, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Dynamical Analysis of Three Distant Trans-Neptunian Objects with Similar Orbits

T. Khain¹, J. C. Becker²,<fnref n="37" />, F. C. Adams¹,², D. W. Gerdes¹,², S. J. Hamilton¹,<fnref n="37" />, K. Franson¹, L. Zullo¹, M. Sako³,

K. Napier¹, Hsing Wen Lin (林省文)¹, L. Markwardt²,<fnref n="37" />, P. Bernardinelli³, T. M. C. Abbott⁴, F. B. Abdalla⁵,⁶, J. Annis⁷,

S. Avila⁸, E. Bertin⁹,¹⁰, D. Brooks⁶, A. Carnero Rosell¹¹,¹², M. Carrasco Kind¹³,¹⁴, J. Carretero¹⁵, C. E. Cunha¹⁶,

L. N. da Costa¹¹,¹², C. Davis¹⁶, J. De Vicente¹⁷, S. Desai¹⁸, H. T. Diehl¹, P. Doel⁶, T. F. Eifler¹⁹,²⁰, B. Flaugher⁷, J. Frieman⁷,²¹,

J. García-Bellido²², D. Gruen¹⁶,²³, R. A. Gruendl¹³,¹⁴, G. Gutierrez⁷, W. G. Hartley⁶,²⁴, D. L. Hollowood²⁵,

K. Honscheid²⁶,²⁷, D. J. James²⁸, E. Krause¹⁹, K. Kuehn²⁹, N. Kuropatkin⁷, O. Lahav⁶, M. A. G. Maia¹¹,¹², F. Menanteau¹³,¹⁴,

R. Miquel¹⁵,³⁰, B. Nord⁷, R. L. C. Ogando¹¹,¹², A. A. Plazas²⁰, A. K. Romer³¹, E. Sanchez¹⁷, V. Scarpine⁷,

R. Schindler²³, M. Schubnell¹, I. Sevilla-Noarbe¹⁷, M. Smith³², M. Soares-Santos³³, F. Sobreira¹¹,³⁴, E. Suchyta³⁵,

M. E. C. Swanson¹⁴, G. Tarle¹, V. Vikram³⁶, A. R. Walker⁴, W. Wester⁷, and Y. Zhang⁷

(DES Collaboration)

¹Department of Physics, University of Michigan, Ann Arbor, Mi 48109, USA; talikh@umich.edu

²Department of Astronomy, University of Michigan, Ann Arbor, Mi 48109, USA

³Department of Physics and Astronomy, University of Pennsylvania, Philadelphia, PA 19104, USA

⁴Cerro Tololo Inter-American Observatory, National Optical Astronomy Observatory, Casilla 603, La Serena, Chile

⁵Department of Physics and Electronics, Rhodes University, PO Box 94, Grahamstown, 6140, South Africa

⁶Department of Physics & Astronomy, University College London, Gower Street, London, WC1E 6BT, UK

⁷Fermi National Accelerator Laboratory, P.O. Box 500, Batavia, IL 60510, USA

⁸Institute of Cosmology & Gravitation, University of Portsmouth, Portsmouth, PO1 3FX, UK

⁹Sorbonne Universités, UPMC Univ Paris 06, UMR 7095, Institut d'Astrophysique de Paris, F-75014, Paris, France

¹⁰CNRS, UMR 7095, Institut d'Astrophysique de Paris, F-75014, Paris, France

¹¹Laboratório Interinstitucional de e-Astronomia—LIneA, Rua Gal. José Cristino 77, Rio de Janeiro, RJ—20921-400, Brazil

¹²Observatório Nacional, Rua Gal. José Cristino 77, Rio de Janeiro, RJ—20921-400, Brazil

¹³Department of Astronomy, University of Illinois at Urbana-Champaign, 1002 W. Green Street, Urbana, IL 61801, USA

¹⁴National Center for Supercomputing Applications, 1205 West Clark St., Urbana, IL 61801, USA

¹⁵Institut de Física d'Altes Energies (IFAE), The Barcelona Institute of Science and Technology, Campus UAB, E-08193 Bellaterra (Barcelona) Spain

¹⁶Kavli Institute for Particle Astrophysics & Cosmology, P.O. Box 2450, Stanford University, Stanford, CA 94305, USA

¹⁷Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas (CIEMAT), Madrid, Spain

¹⁸Department of Physics, IIT Hyderabad, Kandi, Telangana 502285, India

¹⁹Department of Astronomy/Steward Observatory, 933 North Cherry Avenue, Tucson, AZ 85721-0065, USA

²⁰Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Dr., Pasadena, CA 91109, USA

²¹Kavli Institute for Cosmological Physics, University of Chicago, Chicago, IL 60637, USA

²²Instituto de Física Teórica UAM/CSIC, Universidad Autónoma de Madrid, E-28049 Madrid, Spain

²³SLAC National Accelerator Laboratory, Menlo Park, CA 94025, USA

²⁴Department of Physics, ETH Zurich, Wolfgang-Pauli-Strasse 16, CH-8093 Zurich, Switzerland

²⁵Santa Cruz Institute for Particle Physics, Santa Cruz, CA 95064, USA

²⁶Center for Cosmology and Astro-Particle Physics, The Ohio State University, Columbus, OH 43210, USA

²⁷Department of Physics, The Ohio State University, Columbus, OH 43210, USA

²⁸Harvard-Smithsonian Center for Astrophysics, Cambridge, MA 02138, USA

²⁹Australian Astronomical Observatory, North Ryde, NSW 2113, Australia

³⁰Institució Catalana de Recerca i Estudis Avançats, E-08010 Barcelona, Spain

³¹Department of Physics and Astronomy, Pevensey Building, University of Sussex, Brighton, BN1 9QH, UK

³²School of Physics and Astronomy, University of Southampton, Southampton, SO17 1BJ, UK

³³Brandeis University, Physics Department, 415 South Street, Waltham MA 02453, USA

³⁴Instituto de Física Gleb Wataghin, Universidade Estadual de Campinas, 13083-859, Campinas, SP, Brazil

³⁵Computer Science and Mathematics Division, Oak Ridge National Laboratory, Oak Ridge, TN 37831, USA

³⁶Argonne National Laboratory, 9700 South Cass Avenue, Lemont, IL 60439, USA

Received 2018 October 5; revised 2018 October 22; accepted 2018 October 22; published 2018 November 21

Abstract

This paper reports the discovery and orbital characterization of two extreme trans-Neptunian objects (ETNOs), 2016 QV89 and 2016 QU89, which have orbits that appear similar to that of a previously known object, 2013 UH15. All three ETNOs have semimajor axes *a* ≈ 172 au and eccentricities *e* ≈ 0.77. The angular elements (i, ω, Ω) vary by 6°, 15°, and 49°, respectively, between the three objects. The two new objects add to the small number of TNOs currently known to have semimajor axes between 150 and 250 au, and they serve as an interesting dynamical laboratory to study the outer realm of our solar system. Using a large ensemble of numerical integrations, we find that the orbits are expected to reside in close proximity in the (a, e) phase plane for roughly 100 Myr before diffusing to more separated values. We find that an explanation for the orbital configuration of the bodies as a collision product is disfavored. We then explore other scenarios that could influence their orbits. With aphelion distances over 300 au, the orbits of these ETNOs extend far beyond the classical Kuiper Belt and an order of magnitude beyond Neptune. As a result, their orbital dynamics can be affected by the proposed new solar system member, referred to as Planet Nine in this work. With perihelion distances of 35–40 au, these orbits are also influenced by resonant interactions with Neptune.

<footnote n="37">³⁷ NSF Graduate Research Fellow.</footnote>
A full assessment of any possible new solar system planets must thus take into account this emerging class of TNOs.

*Key words:* Kuiper belt objects: individual (2016 QV89, 2016 QU89, 2013 UH15)

# 1. Introduction

In addition to its major planets, the solar system contains a vast number of small, rocky bodies with a variety of orbital elements. The orbits of these minor bodies provide an important record of the history of our solar system. In particular, the history of impacts and binary dissociations, which leave minimal observable traces, can be discerned through dynamical techniques in the asteroid belt (e.g., Nesvorný et al. 2002), among Jovian satellites (e.g., Nesvorný et al. 2003), and in the Kuiper Belt (e.g., the case of Haumea; Brown et al. 2007).

This paper reports the discovery and dynamical characterization of two new extreme trans-Neptunian objects and a third previously known body, all of which exhibit similar orbital elements. Unfortunately, the discovery and characterization of objects with common histories becomes more difficult in the outer solar system (beyond Neptune), where the surface density of known objects is much lower than in the asteroid belt, and objects are often observable only near perihelion. Although an initial velocity dispersion of a few hundred meters per second will result in orbits that disperse quickly, the semimajor axes, eccentricities, and inclinations of such objects are expected to remain differentiable from the background (see Figures 2 and 3 of Marcus et al. 2011). As a result, groups of objects with similar orbital elements require additional study to discern whether or not they actually have a common origin. Haumea is the best example of a Kuiper Belt collisional family (Ragozzine & Brown 2009; Schlichting & Sari 2009; Volk & Malhotra 2012). Since its initial discovery, more family members (which have orbits similar to each other) have been found and confirmed.

Among the population of trans-Neptunian objects (TNOs), unidentified families certainly exist. Recent progress toward identifying such associated objects in the outer solar system comes from the work of de la Fuente Marcos & de la Fuente Marcos (2018), who performed a correlation analysis on outer solar system objects and found a number of potentially associated objects. One leading candidate for an associated pair of extreme TNOs (ETNOs) beyond 150 au is the case of 2004 VN112 and 2013 RF98 (de León et al. 2017), which have been proposed to have come from a binary dissociation event (de la Fuente Marcos et al. 2017). The identification of additional associated objects in the solar system beyond Neptune will enable a better understanding of binary dissociation mechanisms at all orbital locations.

In recent years, dedicated (Bernstein et al. 2004; Petit et al. 2017) and serendipitous surveys (Gaudi 2004; Liu et al. 2015; Gerdes et al. 2017) have allowed for the discovery of many new objects with more distant orbits than was previously thought possible (Allen et al. 2001), allowing for the identification of new trends. Trujillo & Sheppard (2014) observed an alignment in argument of perihelion for the most distant TNOs, and Batygin & Brown (2016) subsequently pointed out an additional alignment in longitude of perihelion. This clustering was used as evidence for the Planet Nine hypothesis. Since then, an additional group of eight ETNOs in this class (with $a > 250$ au) has been found. Less distant

TNOs, with semimajor axes between 150 and 250 au, can also provide insight toward the Planet Nine hypothesis, as it remains unclear exactly where the demarcation between the TNOs influenced mainly by Neptune and those influenced primarily by Planet Nine should lie. As a result, the identification of new ETNOs with semimajor axes in the range $a = 150–250$ au not only expands our census of the outer solar system, but also provides further constraints on the Planet Nine hypothesis.

In this paper, we use data from the Dark Energy Survey (DES) to discover two new objects in the outer Kuiper Belt with large semimajor axes. Motivated by the importance of identifying potential family candidates in the outer solar system, we examine the apparent similarities between their orbits with each other and with another previously discovered (Sheppard & Trujillo 2016b) Kuiper Belt Object (KBO). With semimajor axes $a \sim 170$ au and aphelion distances $\sim 300$ au, these new objects add to the growing inventory of distant objects in the solar system. Motivated by the uniqueness of the large orbital distances (as these objects are the 10th and 11th known objects to have semimajor axis between 150 and 250 au), we also evaluate their orbital dynamics.

Section 2 describes the observational results, including the methods used in DES and the analysis that specifies the orbital elements of the newly discovered bodies. The dynamics of these objects are studied in Section 3, which considers whether the apparent similarity in their orbits is potentially due to random chance. In Section 4, we investigate the role of mean-motion resonances in the dynamics of these objects, considering both possible resonances with Neptune and with the proposed Planet Nine. The paper ends in Section 5 with a summary and a look to the future.

# 2. Discovery

The DES (Dark Energy Survey Collaboration et al. 2016) is an optical survey that observes nearly 5000 square degrees of the southern sky using the 4 m Blanco telescope at the Cerro Tololo Inter-American Observatory in Chile. DES had a nominal survey allocation of 525 nights over five years from 2013 September through 2018 February; a portion of a sixth season is planned for late 2018. DES uses the Dark Energy Camera (DECam; Flaugher et al. 2015), a prime-focus camera with a 3 square degree field of view and a focal plane consisting of 62 separate $2K \times 4K$ red-sensitive CCDs. There are two distinct DES survey modes: the wide survey and the supernova survey (Bernstein et al. 2012). Wide-field survey fields are observed roughly two to four times per six-month observing campaign in each of the five *grizY* photometric bands, to nominal single-exposure $10\sigma$ depths of $g = 23.6$, $r = 23.3$, $i = 22.8$, $z = 22.1$, and $Y = 20.7$ (Morganson et al. 2018), with the result being eventual periodic complete coverage of the 5000 square degree survey area in each wavelength band. In contrast, the supernova survey comprises 10 distinct DECam pointings where repeated observations in the *griz* bands are made roughly every 6 days. The 10 supernova fields are small (3 square degrees each) compared to the wide field (which makes up the remainder of the 5000 square degrees of the survey area), but their dense observing
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Paths of 2016 QV<sub>89</sub> and 2016 QU<sub>89</sub> over the course of the entire (to-date) DES survey. Nights on which the objects were detected by DES are shown as dots along each trajectory. The objects were initially identified in the S1, S2, and X1 supernova fields, which are indicated by outlines of the DECam focal plane. Additional observations were subsequently added from wide-survey exposures outside these fields.</figcaption></figure>

cadence and somewhat greater depth make them well suited for
moving-object searches. So far, TNO discoveries have come
from both supernova survey (Gerdes et al. 2016) and wide
survey (Gerdes et al. 2017; Becker et al. 2018; Lin et al. 2018)
DES data.

## 2.1. Detection of 2016 QV<sub>89</sub> and 2016 QU<sub>89</sub>

The original detection of these objects came in supernova
survey data from the DES 2016–17 observing campaign (DES
Year 4; Diehl et al. 2018). Data from these fields were
processed with the DES difference imaging pipeline (Kessler
et al. 2015). A template image was subtracted from each new
search image, and statistically significant sources in the
subtracted image were identified. Artifacts and non-point
spread function-like sources were rejected using the techniques
in Goldstein et al. (2015), resulting in a relatively clean catalog
of single-epoch transients. We then identified pairs of
detections within 20 nights of each other whose angular
separation was consistent with what would be expected given
the predicted earth-reflex motion of a distant object. Once a
database of pairs had been constructed, we linked the pairs into
chains of observations by testing the goodness of fit of the best-
fit orbit to each chain and requiring the reduced chi-squared
χ²/N < 2 (Bernstein & Khushalani 2000) to qualify as a
detection. The reduced chi-squared of the fit was 1.75 for 2016
QV<sub>89</sub>, 1.3 for 2016 QU<sub>89</sub>, and 1.1 for 2013 UH<sub>15</sub>.

After the initial barycentric orbit was determined for both
objects, we searched for additional detections in wide-survey
and supernova exposures from epochs both before and after the
discovery opposition. Indeed, 2016 QU89 was found to have
appeared in a supernova field on two nights in 2013 October,
less than a month after the start of the survey. The trajectories
of these objects over the full five years of the survey are shown
in Figure 1. The full five-year DES data set allows both orbital
arcs over multiple oppositions, resulting in total arc lengths<fnref n="38" />
of 1481 and 1537 days for 2016 QV89 and 2016 QU89,
respectively.

## 2.2. Physical Parameters

After data processing, we found the best-fit orbital elements
for the two newly discovered bodies using the fitting algorithm
from Bernstein & Khushalani (2000). We computed refined
astrometric positions for each wide-field observation using the
WCSFIT software (Bernstein et al. 2017), which provides
astrometric solutions referenced to the GAIA DR1 catalog
(Gaia Collaboration et al. 2016). This process included
corrections for the effects of distortions on the DECam CCDs,
as well as for chromatic terms from lateral color and differential
atmospheric refraction. The fit for 2016 QV89 used a series of
99 total exposures taken between 2013 December 21 and 2018
January 11. Object 2016 QU89 was fit using a series of 39
observations taken between 2013 October 12 and 2017
December 26. For consistency, we also refit the orbital
elements for 2013 UH15 using the 10 available observations.
The best-fit values, along with magnitudes and colors derived
from DES photometric measurements, are presented in Table 1.
The colors for 2016 QV89 and 2016 QU89 were computed by
subtracting simultaneous (defined as taken on the same night,
to account for systematic variability between nights that would
lead to subtracted colors drifting when measurements from
multiple nights were combined) measurements for each color
and averaging across all measurements for an object. We
ignore object rotation because each supernova field exposure
sequence analyzed here consists of a set of griz exposures and
lasts just 14 minutes in total. This is much less than the typical
rotation period of a TNO, so the color measurements at each
epoch can be regarded as essentially simultaneous. The colors
of 2016 QV89 and 2016 QU89 are consistent with each other to
within 1σ for all four wavebands.

The orbital elements of 2016 QV<sub>89</sub> and 2016 QU<sub>89</sub> are
remarkably similar to each other and also to known object 2013
UH<sub>15</sub> (Sheppard & Trujillo 2016a). These three objects have
quantities (*a*, *e*) varying by less than about 2% among the three
bodies, and the remaining angles (i, ω, Ω) are also similar: in
particular, the orbits of 2016 QV<sub>89</sub> and 2013 UH<sub>15</sub> appear to
be nearly perfectly aligned. Figure 2 shows the orbits of the
three bodies. The top panel shows all of the known ETNOs
with semimajor axes *a* > 150 au and perihelion distances

<footnote n="38">³⁸ The arc length of an object is the number of days between the earliest and latest observations of the object.</footnote>
<table><caption>**Table 1**

Orbital Elements of the Three ETNOs Considered in This Work</caption><thead><tr><th>Parameter</th><th>2013 UH<sub>15</sub></th><th>2016 QV<sub>89</sub></th><th>2016 QU<sub>89</sub></th></tr></thead><tbody><tr><td><i>a</i></td><td>173.6 ± 1.7 au</td><td>171.70 ± 0.05 au</td><td>171.40 ± 0.02 au</td></tr><tr><td><i>e</i></td><td>0.798 ± 0.002</td><td>0.76731 ± 0.00007</td><td>0.79439 ± 0.00002</td></tr><tr><td><i>i</i></td><td>26°081 ± 0°001</td><td>21°38750 ± 0°00007</td><td>16°97552 ± 0°00002</td></tr><tr><td><i>ω</i></td><td>282°87 ± 0°06</td><td>281°093 ± 0°004</td><td>303°337 ± 0°004</td></tr><tr><td>Ω</td><td>176°543 ± 0°001</td><td>173°2150 ± 0°0002</td><td>102°8996 ± 0°0002</td></tr><tr><td>Epoch</td><td>2456594.5804 JD</td><td>2456647.6445 JD</td><td>2456575.6372 JD</td></tr><tr><td>Time of perihelion</td><td>2472269.15 ± 14.12 JD</td><td>2469915.40 ± 0.21 JD</td><td>2459260.79 ± 0.60 JD</td></tr><tr><td>Perihelion</td><td>35.0 ± 0.7 au</td><td>39.95 ± 0.02 au</td><td>35.249 ± 0.007 au</td></tr><tr><td>Aphelion</td><td>312 ± 3 au</td><td>303.45 ± 0.09 au</td><td>307.63 ± 0.03 au</td></tr><tr><td>Absolute magnitude</td><td>7.7</td><td>5.9</td><td>7.95</td></tr><tr><td><i>g</i>–<i>r</i> (mag)</td><td>...</td><td>0.66 ± 0.10</td><td>0.64 ± 0.09</td></tr><tr><td><i>r</i>–<i>z</i> (mag)</td><td>...</td><td>0.43 ± 0.12</td><td>0.44 ± 0.17</td></tr><tr><td><i>i</i>–<i>z</i> (mag)</td><td>...</td><td>0.15 ± 0.18</td><td>0.17 ± 0.18</td></tr><tr><td><i>r</i>–<i>i</i> (mag)</td><td>...</td><td>0.31 ± 0.11</td><td>0.23 ± 0.17</td></tr></tbody></table>

**Note.** The solution for 2013 UH<sub>15</sub> was computed using data from JPL's SSDG SBDB calculated at epoch JD 2458000.5 as written in Table 1 of de la Fuente Marcos & de la Fuente Marcos (2017). The arc lengths of DES observations were 1481 and 1537 days for 2016 QV<sub>89</sub> and 2016 QU<sub>89</sub>, respectively. Colors of 2016 QV<sub>89</sub> and 2016 QU<sub>89</sub> come from DES's multiwaveband observations. As 2013 UH<sub>15</sub> has not been observed by DES and the observations used to compute its orbit were all taken in the *r* band, colors do not yet exist for this object.

<figure><img src="image_4.png" /></figure>

**Figure 2.** The top panel shows the orbits of all currently known ETNOs with $q > 35$ au and $a > 150$ au. The objects with $150 < a < 250$ au are shown in light blue, and objects with $a > 250$ au are in pink. The three objects we discuss in this work (which we call the triplet objects) are boxed and in blue, and the dashed object is the high-inclination and recently discovered extreme TNO 2015 BP<sub>519</sub> (Becker et al. 2018). The bottom panel shows a closer view of the orbital geometries of the triplet objects. The filled-in circles indicate their position on the orbit at discovery (near perihelion).

$q > 35$ au. A close-up of the orbits of the two new objects and their previously discovered counterpart is shown in the bottom panel.

The orbits of these objects are not identical, and similar orbits do not prove a common origin. However, considered in the context that these objects reside in the outer solar system, in a sparsely populated region of orbital parameter space ($a > 150$ au), such similar orbits are unusual and merit further study. The size ratios between the objects is also intriguing: assuming albedos of 10%, 2016 QV<sub>89</sub>, 2016 QU<sub>89</sub>, and 2013 UH<sub>15</sub> have sizes of 280 km, 110 km, and 120 km, respectively. If all three objects are assumed to have the same albedo and density, the two smaller objects (2013 UH<sub>15</sub> and 2016 QU<sub>89</sub>) have about 15% of the mass of the largest object, 2016 QV<sub>89</sub>.

## 2.3. Observational Bias

Due to the small survey area (the relatively small sky-area three-square-degree supernova fields) where these objects were initially detected, it is possible that the apparent similarity in their orbital elements comes from bias in their detection locations. If this were the case, we would expect TNOs with the (*a*, *e*, *i*) of these objects to be more easily recovered by our detection pipeline. This would mean that the association between the objects is only due to a bias in the sensitivity of our pipeline.

To test this, we simulate roughly 440,000 clones of objects dynamically similar to those considered in this work. The orbital elements of these clones are drawn from uniform distributions with semimajor axis $150 < a < 1000$ au, perihe- lion distance $q > 30$ au, and inclination *i* uniformly distributed between 0° and 180°. Choosing this population of clones allows us to identify the selection biases of DES; this initial distribution is not meant to represent a realistic Kuiper Belt population. Using the computed orbits of these objects, we then determine which orbits could be detected in DES observations in the first three observing campaigns. First, we eliminate objects whose sky positions or apparent magnitudes made them unobservable by DES. Then, of the remaining clones, we test which would be observable and linkable by our pipeline. Clones are considered linkable in the supernova data if they are detectable in three or more DES exposures (on different nights), with neighboring exposures separated by 20 nights or fewer.
<figure><img src="image_5.png" /><figcaption>**Figure 3.** DES selection function for objects with perihelion distances greater than 30 au and semimajor axis greater than 150 au and smaller than 1000 au. The orbital elements of 2016 QV<sub>89</sub> and 2016 QU<sub>89</sub> are denoted by triangles. Note that the two triangles are overlapping in the top panel, due to the objects' similarity in semimajor axis.</figcaption></figure>

A total of 6446 unique clones were determined to be observable using these criteria, and those surviving clones are plotted in Figure 3. These histograms show the distribution of instantaneous orbital elements for linked clones, as a fraction of the input uniform distributions. Since the resulting distributions are fairly flat and do not show any strong preference for orbital elements around the (a, e, i) values of 2016 QV<sub>89</sub> and 2016 QU<sub>89</sub>, we can be confident that our detection of two objects with such similar orbits is not due to observational bias.

## 3. Dynamics and Orbital Similarities

In order to study the dynamics of these three objects, we employ a numerical model, as these ETNOs fall in the Neptune scattering region (with perihelion distances in the $q \sim 35-45$ au range). As a result, their behavior can be driven by resonant effects (Saillenfest et al. 2017; Volk et al. 2018), and their energy diffusion can be driven by close encounters with Neptune (Duncan et al. 1987). Because the evolution of these objects depends on the aforementioned short-period effects, numerical simulations are an effective tool for understanding their dynamics (for comparisons of numerical and analytic treatments in related contexts, see, e.g., Becker & Adams 2017; He & Petrovich 2018). In this work, we use N-body code *Mercury6* to evaluate the evolution of these three objects in the presence of the four giant planets. We exclude the terrestrial planets and use a time step of 20 days, with a hybrid symplectic and Bulirsch–Stoer integrator. In each integration, we include test particle clones of each of the three TNOs, with orbital elements for each drawn from the covariance matrices resulting from our fits and corresponding to the values in Table 1. We integrate the orbital elements for each TNO to a common epoch before beginning the simulations. Simulations are run in two batches, one forward and one backward in time for 4.5 Gyr each. These simulations allow us to study the orbital evolution of these three objects (which we call the “triplet” objects) and assess their similarity. Before considering the potential association between the three TNOs, however, we explore their dynamics individually.

### 3.1. Individual Dynamics

Although all three of the triplet objects currently have long-period orbits, these ETNOs are not decoupled from the rest of the solar system as their perihelia are bound to Neptune. This effect is evident in Figure 4, which shows the evolution of the semimajor axis, eccentricity, and inclination of a representative subset of triplet object clones. Due to scattering interactions with Neptune, the semimajor axis of the triplets changes rapidly and their eccentricities grow. As a result, not all clones survive the entire solar system lifetime. In fact, only 36% of the clones of 2016 QU<sub>89</sub> survive the full 4.5 Gyr simulation; similarly, only 45% of the clones of 2013 UH<sub>15</sub> survive. Interestingly, the clones of 2016 QV<sub>89</sub> are significantly more stable, with 89% of clones surviving, possibly because of its larger perihelion distance. The stability fractions for the backward integrations are similar, as expected, with 29% of clones surviving for 2016 QU<sub>89</sub>, 47% for 2013 UH<sub>15</sub>, and 82% for 2016 QV<sub>89</sub>.

To characterize this behavior in more detail, we consider how the orbital elements of clones of a single object change over time. For example, in Figure 5, we present several snapshots of the eccentricity distribution of the surviving clones of 2016 QU<sub>89</sub>. Initially, this distribution is tightly centered around the best-fit value. As the system evolves, the eccentricity distribution spreads apart, as shown by the three panels in Figure 5. In order to further characterize this behavior, we consider the time evolution of the distribution width of several parameters.

In particular, we focus on the distribution width of the parameter $1/a$ (which scales with orbital energy), orbital eccentricity *e*, and inclination *i*. We compute the standard deviation $\sigma$ of the parameter of interest among all surviving clones of each object at every time step. The resulting time evolution of $\sigma$ is shown in Figures 6 (for $1/a$ and *e*) and 7 (for *i*) on a log-log plot.

The initial distribution of clones of each orbital element resembles a Gaussian, centered at the best-fit value of the parameter. In Figure 6, we show the behavior of $\sigma$ for the parameters $1/a$ and *e*. For the first few million years, the distribution hardly changes, as the orbits of the clones do not evolve on such short timescales. Over the span of four billion years, however, the spread in the clones increases. Linear fits to the blue (2016 QU<sub>89</sub>) and purple (2013 UH<sub>15</sub>) $\sigma$ curves find a slope of $m \approx 0.4$, which indicates that these clones spread apart in a subdiffusive manner. Note that diffusion can be described as a random walk, so quantities should increase with time $\propto t^{1/2}$. One might expect the distribution of orbital elements to behave similarly, as a random walk in phase space, due to scattering interactions with Neptune (subject to an absorbing boundary at $e=1$, where objects are ejected from the solar system). If this were, in fact, a diffusion process, we would expect to see a straight line with slope $m=0.5$ on the log-log plot. However, it is evident that the observed behavior is not a pure random walk, as the time evolution of the $\sigma$ values is somewhat slower.

The behavior of the green $\sigma$ curve (2016 QV<sub>89</sub>) differs from the two others. It appears that the dynamics of 2016 QV<sub>89</sub> are not as dominated by Neptune scattering as are the other two objects. This trend is evident in the much higher number of surviving clones of 2016 QV<sub>89</sub> cited above. The reason for this difference is likely due to the higher perihelion distance of 2016 QV<sub>89</sub> as compared to the other objects. With a perihelion
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Time evolution of the semimajor axis, eccentricity, and inclination for representative clones of each object. The clones are initialized with sets of orbital elements that are statistically equivalent to the best-fit parameters. The clones of 2016 QV<sub>89</sub> are shown in green, 2016 QU<sub>89</sub> in blue, and 2013 UH<sub>15</sub> in purple.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Snapshots of the eccentricity distributions for the surviving clones of 2016 QU<sub>89</sub> at three epochs. At early times (1 Myr), the eccentricities of the clones are tightly centered around the best-fit value (left panel). At later times, this distribution spreads apart, with distributions shown for 100 Myr (middle panel) and 1 Gyr (right panel). The full time evolution of the distribution widths is shown in Figure 6.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Time evolution of the width $\sigma$ of the inclination distribution, where $\sigma$ is defined as in Figure 6. The guiding dashed line has a slope of 1 on this log-log plot, showing that initially the width of the distribution grows, but soon plateaus to a roughly constant value. This indicates that the inclination distribution has reached a point after which it does not significantly change; this distribution is shown in the inset plot for 2016 QU<sub>89</sub>.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 6.** Time evolution of the widths $\sigma$ of the distributions of the orbital parameters $1/a$ and $e$. Here $\sigma$ is defined to be the standard deviation of the distribution of a given orbital element for a population of clones, which have been integrated numerically forward in time. The dashed lines show linear fits to the curves (in the log-log plot) for the time span 0.1 Myr–4.5 Gyr for all three objects. The slopes for 2016 QU<sub>89</sub> and 2013 UH<sub>15</sub> are $m \approx 0.4$, which is somewhat shallower than the benchmark value $m = 0.5$ expected for diffusive behavior. Object 2016 QV<sub>89</sub> does not appear to fit the linear trend well.</figcaption></figure>

of $q \approx 40$ au, this object is more detached from Neptune and experiences fewer (and less disruptive) close encounters. In addition, the evolution of 2016 QV<sub>89</sub> may be affected more strongly by mean-motion resonances with Neptune, a possibility that we discuss in Section 4.1.

In Figure 7, we show the time evolution of the width of the inclination distribution for the three objects. Initially, the width grows, roughly following a line with a slope of unity on the log-log plot, but appears to saturate on a timescale of order 0.1 Gyr, indicating that the inclination distribution remains approximately the same for the rest of the evolution. Notice that this last decade and a half in time corresponds to most of the age of the solar system. This final inclination distribution is shown for 2016 QU<sub>89</sub> in the inset plot; most of the distribution falls in the range 15°–25°.

These results (in particular, Figures 4, 6, and 7) show that the future orbits of these objects diverge steadily on relatively rapid timescales. As a result, any similarity in their current-day orbits must be due to either a recent event (such as fragmentation or a binary dissociation event), the past attainment of an orbital resonance that would force the objects to maintain similar orbits, or random chance. Taking this dynamical evolution into account, we now consider possible explanations for the similarity of the orbits of these three objects. In the next section, we evaluate the probability that the orbits of these objects are similar by random chance.
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Plot showing the current metric distance between the objects 2016 QV89, 2016 QU89, and 2013 UH15 (red vertical line) as compared to the distances between random sets of three objects drawn from several distributions. The green population (MPC Set 1) contains all of the TNOs reported to the Minor Planet Center with $a > 50$ au, $e > 0.45$, and $q > 30$ au; the purple population (MPC Set 2) contains all of the TNOs reported to the Minor Planet Center; and the blue population is the test particle distribution from de la Fuente Marcos et al. (2017), where $150 < a < 800$ au, $0.7 < e < 0.95$, $i < 55^{\circ}$, and $30 < q < 90$ au. As we can see, the triplet objects discussed in this work are currently closer together in phase space than most objects, for all three distributions considered. For all MPC objects, we used barycentric orbital elements with solution epoch 2458200.5 JD, which were converted from heliocentric orbital elements coming from an MPC file downloaded on 2018 June 9.</figcaption></figure>

## 3.2. Comparison of Orbits

To evaluate the likelihood that the orbits of these objects are similar merely by coincidence, we compare the orbital similarity of the triplets to that of randomly chosen groups of three objects drawn from several control distributions.

First, we define a dimensionless distance metric to characterize a set of three objects. This function $d(t)$ is taken to be the sum of the pairwise distances in the space of orbital elements. Specifically, we use the sum of Euclidean metrics acting on the scaled elements $a$, $e$, and $i$ for each pair of TNOs:

$$d(t) = \sum_{j \neq k} \left[ \left( \frac{a_j(t) - a_k(t)}{\frac{1}{2}(a_{j0} + a_{k0})} \right)^2 + \left( \frac{e_j(t) - e_k(t)}{\frac{1}{2}(e_{j0} + e_{k0})} \right)^2 + \left( \frac{i_j(t) - i_k(t)}{\frac{1}{2}(i_{j0} + i_{k0})} \right)^2 \right], \quad (1)$$

where $j$ and $k$ denote the objects, and the subscript 0 denotes the initial orbital elements of the objects.

Next, we choose a clone of each of 2016 QV89, 2016 QU89, and 2013 UH15 from the measured orbit posteriors and compute the above metric distance at the current epoch. Repeating this process for 125,000 distinct combinations of clones, we compute the current average distance in our triplet.

From this analysis, we have a measure of the current-day orbital similarity of these three objects. To explore whether this distance is remarkable among the larger set of objects in the outer solar system, we compare the current average distance of the triplet objects to several different representative distributions of TNOs, shown in Figure 8. Each of the distributions was created by drawing 125,000 sets of three objects from the specified population of TNOs and computing their current-day metric distance. The green distribution (MPC Set 1) contains all of the TNOs reported to the Minor Planet Center with $a > 50$ au, $e > 0.45$, and $q > 30$ au; the purple population (MPC Set 2) contains all of the distant objects reported to the Minor Planet Center (we note that this includes Neptune trojans); and the blue population is the uniform test particle distribution from de la Fuente Marcos et al. (2017), with $150 < a < 800$ au, $0.7 < e < 0.95$, $i < 55^{\circ}$, and $30 < q < 90$ au.

The red vertical line in Figure 8 shows the current-day average distance between our triplet objects. It is clear that this red line falls on the far left of all three TNO distributions considered. For MPC Set 1, 4023 members of the 125,000 sets in the control population, or 3.2%, were more similar than the observed triplet. For MPC Set 2, 96 of 125,000 sets, or 0.077%, were more correlated; finally, for the distribution drawn from de la Fuente Marcos et al. (2017), 1137 out of 125,000 sets, or 0.9096%, were more correlated than the triplet. The robust nature of this result leads us to conclude that the current-day similarity of 2016 QV89, 2016 QU89, and 2013 UH15 is highly unusual, even when compared to a variety of background populations.

Considering the distance between objects, as outlined above, may provide an effective method of identifying other triplet sets from the whole population of KBOs. Toward this end, we apply our analysis to the dwarf planet Haumea and its associated family members. Although a full characterization of the dynamics of the Haumea system is beyond the scope of this work, it provides a useful test of our method. The Haumea system contains two moons, which are known from previous literature to be associated (Brown et al. 2007; Levison et al. 2008; Schaller & Brown 2008; Ragozzine & Brown 2009; Schlichting & Sari 2009; Leinhardt et al. 2010; Volk & Malhotra 2012).

For this test, we compute the current-day metric distance among all possible sets of three objects drawn from the Haumea family. We find that the Haumea triplets have distances directly comparable to the values computed for our triple system. The metric distance for triplets of objects drawn from the Haumea family ranges from 0.0024 to 1.28, showing that the correlated Haumea family is much more similar by our metric than the general control populations considered in Figure 8. One should keep in mind, however, that the dynamical environment of the Haumea family is significantly different from that of the new triplet. Jupiter has a much larger effect on the evolution of Haumea due to its closer proximity. Moreover, it is possible that the breakup of the original dwarf planet into its present family members was caused by rotational fission (Snodgrass et al. 2010; Ortiz et al. 2012). As a result, comparisons between these two sets of objects can only be made at the order-of-magnitude level.
<figure><img src="image_11.png" /><figcaption>Figure 9. Plot showing the current metric distance for semimajor axis only (Equation (2)) for the three objects considered in this work, as compared to the distances between random sets of three objects drawn from our control distributions (MPC Set 1 contains all of the TNOs reported to the Minor Planet Center with $a > 50$ au, $e > 0.45$, $q > 30$ au; the de la Fuente Marcos (dlFM) distribution is the test particle distribution from de la Fuente Marcos et al. 2017, where $150 < a < 800$ au, $0.7 < e < 0.95$, $i < 55^{\circ}$, and $30 < q < 90$ au). The inset shows a zoom-in on the left side of the distribution, with the red vertical line indicating the current distance among our three triplet objects.</figcaption></figure>

The main conclusion from the distance analysis of this section is that the triplet objects (2016 QV89, 2016 QU89, and 2013 UH15) are statistically more correlated than the general population of TNOs. Using the distance metric from Equation (1), we compare their orbital elements to those derived for three different control populations. For each population considered, the triplet objects are among the most correlated; however, we cannot exclude that the similarity in overall orbits is due to random chance.

Since the similarity between the semimajor axes of the triplets is the most unusual of all the parameters, we consider one more test. Recomputing the metric we used previously for only semimajor axis as

$$d_a(t) = \sum_{j \neq k} \left( \frac{a_j(t) - a_k(t)}{\frac{1}{2}(a_{j0} + a_{k0})} \right)^2 \quad (2)$$

yields the histograms in Figure 9. In this figure, we compare MPC Set 1 and the de la Fuente Marcos distribution, which have eccentricities similar to our objects.

We find that, compared to the present-day distance between our triplet objects, 0.08% of trials in Set 1 are more similar than our triplet and 0.01% of trials in the dlFM set are more similar. We exclude Set 2 because it includes Neptune trojans, which we expect to have fixed, identical values of semimajor axis.

The analysis based on semimajor axis alone describes these objects as more remarkable than does the full metric using all orbital elements. The existing population of TNOs with semimajor axis values greater than 150 au remains small, with these three objects representing roughly 10% of the known total at the time this paper was written. As such, it seems striking that these three would reside in the same ~2 au of such a large parameter space. One explanation that can force objects to attain particular values of semimajor axis (but affect the objects' eccentricities and inclinations less directly) is orbital resonance, which we discuss in the next section.

## 4. 2016 QV89, 2016 QU89, and 2013 UH15 as Resonant Objects

To explore the possibility that the semimajor axis similarity of the triplets is explained by orbital resonances, we evaluate

the likelihood of these objects falling into resonances with other solar system bodies, specifically, Neptune and the proposed new solar system member Planet Nine. We find that the resonant dynamics in each of these two cases are distinct.

### 4.1. Resonances with Neptune

Thus far, we have discussed only the scattering interactions that these three objects experience due to Neptune. In theory, however, these ETNOs, although distant, may be affected by mean-motion resonances with Neptune. In fact, recent work by Volk et al. (2018) has identified two objects with semimajor axes of 129.8 au and 129.9 au as living in the 9:1 resonance with Neptune. These two objects are currently the most distant TNOs known to reside in Neptune resonances.

Inspired by this finding, we perform a similar analysis on the clones of our three objects. We identify intervals of time in which the period ratio between a clone and Neptune is approximately constant, and we then test resonances up to 29th order, considering resonances that fall into the period ratio interval ($P_{\text{Nep}}/P_{\text{TNO}}$) ± 0.1. For each resonance argument, we generate plots of the time evolution and identify intervals of libration by finding low-point-density regions in the plots.

The output of the above analysis allows us to identify the time intervals during which a clone is securely librating in a resonance. Interestingly, we find that these three objects often experience resonant interactions with Neptune, which usually last tens of millions of years. Summing over the evolution of all of the clones, we find that 2016 QV89 is resonant most often, with 21% of the total integration time spent in resonance, while 2016 QU89 and 2013 UH15 spend 8% and 14% of the time in resonance, respectively (see Figure 10).

The most populated resonances for these three objects are shown in Figure 11. Here we take the 10 most populated resonances for each of the three objects and then consider their union. Since the most populated resonances are not the same for all three bodies, we plot 17 resonances in total. As expected, most of these resonances are quite high order. For instance, the most populated resonance for 2016 QV89 is the 2:27 commensurability, which corresponds to a semimajor axis of about $a = 170$ au, close to its current-day value. An example of such a resonant instance is shown in Figure 12.
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Fraction of time each object considered in this work spends in a Neptune resonance in a 4.5 Gyr backward integration. To compute the fraction of time, we sum the total time spent in resonance by all of the clones of an object and divide by the total survival time of those clones. The most resonant object of the three is 2016 QV<sub>89</sub>, which spends about 20% of its time in resonance.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 11.** Histograms showing the time spent in resonant configurations with Neptune. The bins on the *x* axis represent the union of the 10 most populated resonances for each object. The fraction of time spent in a resonance is computed by summing over all clones of an object, using the same procedure as in Figure 10 (see text).</figcaption></figure>

It is important to note, however, that these intervals of resonance are transient. None of the clones of these objects are currently in resonance, so they cannot be classified as resonant TNOs according to the standards of Gladman et al. (2008). Nonetheless, the fact that these objects are likely to have been in Neptune resonances in the past provides an additional explanation for their orbital similarity. As shown in Figure 12, even after true resonance is lost, a TNO may reside in a semimajor axis value close to the (previously) resonant value for some time. It is possible that these objects were recently in the same Neptune resonance and are presently evolving away from that state, explaining their current similarity in semimajor axis.

<figure><img src="image_14.png" /><figcaption>**Figure 12.** Example of the 2:27 resonance in the backward integrations for 2016 QU<sub>89</sub>. The top panel shows the semimajor axis evolution, and the bottom panel plots the appropriate resonance argument over the same time interval. The time spent in resonance is brief, less than 20 Myr. Note that time interval for the libration of $\phi_{res}$ (bottom) clearly corresponds to the regime of constant semimajor axis evolution (top).</figcaption></figure>

## 4.2. Resonances with Planet Nine

Recently, observed correlations in argument of perihelion (Trujillo & Sheppard 2014) and then in longitude of perihelion (Batygin & Brown 2016) among the most distant TNOs have led to the proposal of a ninth solar system planet, often called Planet Nine or Planet X. If this proposed planet exists, it would also provide a network of mean-motion resonances in which TNOs could reside. In this section, we consider how the three ETNOs discussed in this work fit into the context of the Planet Nine hypothesis.

As shown in Figure 2, the orbital orientation (as defined by the longitude of perihelion, $\varpi$) of these three objects appears to be antialigned with the proposed orientation of the Planet Nine orbit. This clustering is thus consistent with the orbits of the extreme TNOs first used to infer the existence of Planet Nine and predict its orbit. However, these objects have semimajor axes ($a \sim 171–173$ au) significantly shorter than those of TNOs thought to be dominated by Planet Nine-induced dynamics. (Note that different models for Planet Nine imply different inner boundaries (in *a*) for the region affected by the putative planet. The most recent analysis, Batygin & Morbidelli 2017, shows that objects with $a > 250$ au are likely to reside in the regime influenced by Planet Nine.) As a result, despite the consistent orbital properties of these three ETNOs, it is unclear whether these objects could be part of the securely antialigned population.

In addition to maintaining the apsidally antialigned population, Planet Nine has been shown to generate TNO orbits that are apsidally aligned with the orbit of Planet Nine and others that experience apsidal circulation (Batygin & Brown 2016; Batygin & Morbidelli 2017; Khain et al. 2018). In the presence of Planet Nine, these three objects of interest could fall into either the antialigned or circulating categories. In this section,
we discuss the dynamics of these TNOs in the presence of representative Planet Nine candidates and consider possible mean-motion resonances.

To differentiate between the behaviors outlined above, we define the difference between the longitude of perihelion of Planet Nine and a TNO as

$$ \Delta\bar{\omega} \equiv \bar{\omega} - \bar{\omega}_9, \qquad (3) $$

where $\bar{\omega}$ is the longitude of perihelion of the TNO, and the subscript 9 denotes the orbital elements of Planet Nine. Previous studies have shown that objects that are heavily influenced by Planet Nine do not retain a constant $\Delta\bar{\omega}$ (Beust 2016; Batygin & Morbidelli 2017). Instead, TNOs aligned with Planet Nine experience libration around $\Delta\bar{\omega} \sim 0^{\circ}$, antialigned objects librate about $\Delta\bar{\omega} \sim 180^{\circ}$, and circulating objects have $\Delta\bar{\omega}$ that circulates (by definition) through all values in $[0^{\circ}, 360^{\circ}]$.

To determine which class of behavior a given TNO experiences, it is not enough to know only its current $\Delta\bar{\omega}$; we need to integrate its orbit forward in the presence of Planet Nine and the current solar system and analyze the time evolution of its $\Delta\bar{\omega}$. To evaluate the effect of Planet Nine on the orbital similarity of these three ETNOs, we use numerical simulations similar to those run in the previous section, but now including Planet Nine.

Different incarnations of Planet Nine lead to distinct behaviors for these three ETNOs. Since the orbital elements of Planet Nine are not yet well constrained, we vary its orbit over a range of parameter space (as in Becker et al. 2017) to examine each class of possible interactions between Planet Nine and the three ETNOs under consideration. As shown below, we examine each of the possible dynamical classes that these objects could belong to by adjusting Planet Nine’s orbit, and we focus on the question of the orbital similarity of the triplet objects in each of these cases.

It is important to note that an exhaustive survey of parameter space, as well as determining the exact limits of the Planet Nine parameters that cause each type of behavior considered, is an interesting task that is beyond the scope of this work. Instead, we consider two test cases to illustrate—but not comprehensively study—the regimes of possible behavior.

In Section 4.2.1, we study a Planet Nine candidate with semimajor axis $a_9 = 315$ au, eccentricity $e_9 = 0.5$, and inclination $i_9 = 20^{\circ}$. This orbit lies near the inner edge of the parameter space proposed for Planet Nine and is found to induce apsidal antialignment in the triplet orbits. In Section 4.2.2, we then consider a Planet Nine candidate with orbital elements $a_9 = 505$ au, eccentricity $e_9 = 0.5$, and inclination $i_9 = 20^{\circ}$. This orbit lies near the center of the proposed parameter space and allows for the triplets to be apsidally circulating.

For the numerical work carried out in this section, all simulation parameters (time step and so on) are identical to those of earlier integrations, with the exception of the introduction of a massive body in the form of the one proposed within the framework of the so-called Planet Nine hypothesis.

### 4.2.1. ETNOs with Antialigned $\Delta\bar{\omega}$

Independent of considerations of their orbital similarity, these three objects are currently roughly antialigned with the orbit of the proposed Planet Nine, with a longitude of

perihelion $\bar{\omega} = \omega + \Omega$ that orients their orbits between Sedna and 2012 VP113 (Figure 2). As mentioned above, however, an instantaneous antialigned $\Delta\bar{\omega}$ is not sufficient to sort these objects into the antialigned class, as they could just be opportunistly observed ETNOs whose $\Delta\bar{\omega}$ are truly circulating. In the presence of an appropriately chosen Planet Nine, however, these objects do experience librations in the offset of longitude of perihelion around $\Delta\bar{\omega} \sim 180^{\circ}$.

Given the parameters of Planet Nine, it is possible to approximate the semimajor axis threshold at which the antialignment of the distant Kuiper Belt begins (Batygin & Morbidelli 2017; Becker et al. 2017). That is, depending on the orbit of Planet Nine, we can estimate which outer solar system objects are expected to experience librations about $\Delta\bar{\omega} \sim 180^{\circ}$. Assuming a fixed mass for Planet Nine, the location of this threshold is determined (in part) by the perihelion distance of Planet Nine and is thus a function of $a_9$ and $e_9$. Using these estimates as a guideline, we run a 1 Gyr integration of the clones of the three ETNOs in the presence of the four gas giants and a Planet Nine with orbital elements $a_9 = 315$ au, $e_9 = 0.5$, $i_9 = 20^{\circ}$, $\omega_9 = 150^{\circ}$, $\Omega_9 = 120^{\circ}$, and mean anomaly $M_9 = 180^{\circ}$, which are parameters that fall in the antialigned region for these objects. Indeed, these simulations show that the ETNOs experience libration in the offset of longitude of perihelion around $\Delta\bar{\omega} \sim 180^{\circ}$.

Due to their antialignment, these objects reside in an orbit-crossing region with Planet Nine. In order to avoid close encounters and consequent instability, such ETNOs are often in mean-motion resonances with Planet Nine (Malhotra et al. 2016; Millholland & Laughlin 2017). In fact, they exhibit a more complicated behavior now known as “resonance hopping” (Batygin & Morbidelli 2017; Becker et al. 2017; Hadden et al. 2017), where the objects are often locked into mean-motion resonance but transition into different resonances over the age of the system. It is important to note that mean-motion resonances, which correspond to libration of resonant arguments $\phi$ (e.g., $\phi = (p + q)\lambda_9 - p\lambda - q\bar{\omega}_9$), are distinct from apsidal resonances, which correspond to libration of $\Delta\bar{\omega}$.

As one working example, Figure 13 shows the behavior of a clone of 2016 QV89 over a time interval for which it resides in a 5:2 resonance with Planet Nine. The top panel shows the nearly constant semimajor axis evolution that is characteristic of resonance, and the bottom panel shows the librating resonance argument $\phi$. In this plot, we only show the time interval during which this object is truly in resonance. However, 2016 QV89 often resides in near resonance for billions of years, most likely experiencing nodding behavior (Ketchum et al. 2013) in which it transitions back and forth from a librating to a circulating resonance argument (note that similar behavior has been reported in Millholland & Laughlin 2017 and Becker et al. 2017 for some of the longer-period ETNOs, such as Sedna).

An example of resonance hopping behavior is shown in the bottom panel of Figure 14, which displays the evolution of the triplets in the presence of Planet Nine. It is clear that in this dynamical regime, these objects spend more time with constant semimajor axis than in the transient, highly variable regimes in between. In other words, the triplet objects effectively “hop” from one value of semimajor axis to another.

In contrast, the top panel shows the evolution of the triplet objects in the presence of the current solar system only (without Planet Nine). Here we see that the objects tend to spend more
<figure><img src="image_15.png" /><figcaption>**Figure 13.** Possible 5:2 resonance between 2016 QV<sub>89</sub> and Planet Nine. The top panel shows the time evolution of the semimajor axis for a clone of 2016 QV<sub>89</sub> that lands in a 5:2 resonance with Planet Nine. The bottom panel shows the corresponding resonance angle given by the argument φ = 5λ<sub>9</sub> − 2λ − 3ω<sub>9</sub>. The resonance angle is clearly librating and thus indicates a mean-motion resonance with Planet Nine.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Figure 14.** The top panel shows the rapid semimajor axis evolution of a few representative clones of the triplet objects in the current solar system. The bottom panel shows the constant-a evolution induced by the addition of Planet Nine with a<sub>9</sub> = 315 au, e<sub>9</sub> = 0.5, and i<sub>9</sub> = 20° to the system.</figcaption></figure>

of their time scattering, with more highly variable semimajor
axes, and only spend short time intervals with nearly constant
orbital elements. In other words, rather than continually
hopping between semimajor axis values, the objects only
occasionally get “stuck” into a resonance. As such, this plot

highlights the difference between the effects of “resonance
sticking” (top) and “resonance hopping” (bottom).<fnref n="39" />

In light of this dynamical behavior, perhaps it is not as
surprising to find a collection of ETNOs with similar
semimajor axes. Rather than invoking collisions to explain
the origin of this associated ETNO triplet, an alternate
explanation is that these three objects are currently in the same
mean-motion resonance with Planet Nine. Since resonances
result in semimajor axis oscillations with a finite width, it is not
necessary for the three ETNOs to have identical values of *a* as
long as they are similar. In fact, the observed differences of a
few percent correspond to a relatively narrow libration
amplitude of the resonance argument.

Of course, we cannot know whether these objects are in fact
in resonance with Planet Nine until the planet is discovered (or
ruled out). For a given orbit of Planet Nine, however, we can
search for librating resonance angles among the clones in our
simulations. It is important to note that to achieve antialigned
behavior for these triplet objects, we have chosen Planet Nine
parameters that fall inside the typically accepted range of (a<sub>9</sub>,
e<sub>9</sub>). The original authors—as well as subsequent analyses—
have shown that the perihelion distance of Planet Nine is likely
to be q<sub>9</sub> ~ 250 au, which is significantly larger than our choice
of q<sub>9</sub> ~ 150 au. Taking this difference into account, the next
section considers a more typical set of parameters for Planet
Nine. As a consequence, the triplet objects lose their
antialignment in Δω space.

4.2.2. ETNOs with Circulating Δω

Given the current arguments ϖ of the three objects under
consideration, it is possible that their orbits are in fact
circulating in Δμ and do not remain confined in physical
space. Given their observed semimajor axes, which are smaller
than the typical values for which Planet Nine dominates TNO
dynamics, the Δμ for these ETNO orbits are likely to circulate
even in the presence of Planet Nine. In this case, the currently
observed antialignment with the proposed Planet Nine orbit is
due to chance. Even with circulating Δμ, the dynamical
behavior of the ETNOs could be driven by (1) resonant or
(2) secular interactions with Planet Nine, as discussed below.

**Resonant Case.** In theory, some ETNOs with circulating Δω could be trapped in mean-motion resonances with Planet Nine. Despite their potentially orbit-crossing behavior, these objects would be able to remain stable in the presence of Planet Nine due to the associated phase-space protection mechanism. In this case, the form of the librating resonance angle is different from the antialigned case (Batygin & Morbidelli 2017). An example of an object with circulating Δω but which also has a librating resonance argument is shown in Figure 15.

If the triplet objects are in this dynamical class, then the
similarity of their orbits could be explained by their mean-
motion resonance with Planet Nine, as in Section 4.2.1 for the
antialigned class.

**Nonresonant Case.** Although objects circulating in Δω
may be in mean-motion resonances with Planet Nine, not all
must be in this synchronized state in order to avoid close
encounters with Planet Nine. Depending on the exact orbit of
the planet, it is possible for the objects with circulating Δω
values to follow trajectories that avoid orbit-crossing

<footnote n="39">³⁹ For a recent numerical characterization of resonance sticking in the Kuiper Belt, see the analysis by Yu et al. (2018).</footnote>
<figure><img src="image_17.png" /><figcaption>**Figure 15.** Possible 5:1 resonance between 2016 QV<sub>89</sub> and Planet Nine. The top panel shows the time evolution of the semimajor axis for one clone of 2016 QV<sub>89</sub> that lands in a 5:1 resonance with Planet Nine. The bottom panel shows the corresponding resonance angle given by the argument $\varphi = 5\lambda_0 - \lambda - 3\bar{\varpi}_9 - \bar{\bar{\varpi}}$. The resonance angle is clearly librating and thus indicates a mean-motion resonance with Planet Nine.</figcaption></figure>

configurations. In this case, the behavior of these objects is driven by secular interactions with Planet Nine, and although these objects are not in true resonance, their dynamics are marked by regions of constant semimajor axis evolution.

**5. Conclusion**

This paper has reported the discovery of two new extreme trans-Neptunian objects, 2016 QV<sub>89</sub> and 2016 QU<sub>89</sub>. Our main findings can be summarized as follows:

[1] These two objects may be associated with each other and with the ETNO 2013 UH<sub>15</sub>. The three objects have orbital elements ($a, e$) that differ by less than about 2% among the three bodies, and the remaining angular orbital elements ($i, \omega, \Omega$) fall within ranges of 9°, 20°, and ~70° (see Table 1). For 2016 QV<sub>89</sub> and 2016 QU<sub>89</sub>, the observed colors are statistically indistinguishable (see Table 1). However, we cannot state with certainty that the apparent similarities in orbits are not due to random chance.

[2] The existence of this triplet ETNO system has important ramifications for the dynamics of the outer solar system. Numerical integrations indicate that the three bodies tend to stay together in parameter space for timescales of order of 100 Myr, but diverge over longer spans of time (see Figure 4). In addition, the probability that three random ETNOs happen to lie so close in orbital element space is low (see Figure 8), but this possibility cannot be excluded with current data. These objects should be considered a candidate system for a common origin in future analyses.

[3] These objects are not likely to be in true resonances with Neptune at the current epoch. However, it is possible that they were in resonances in the recent past; this possibility could explain the present-day similarity in their semimajor axes. These objects appear to exhibit “resonance sticking” in our

simulations, generally residing outside of Neptune resonances but “sticking” to resonances for short periods of time.

[4] The existence of this triplet ETNO system has important implications for the Planet Nine hypothesis. The orbits of the three ETNOs under study should be sufficiently distant from the Sun that they would be influenced by the proposed Planet Nine, even though their orbital distances fall below the $a \sim 250$ au cutoff used previously in some works (e.g., Batygin & Morbidelli 2017; Becker et al. 2017; Millholland & Laughlin 2017) to describe the population most susceptible to Planet Nine's influence. In particular, in the presence of Planet Nine, with canonical orbital elements and mass, the long-term evolution of semimajor axes for the three bodies is markedly different, as they exhibit “resonance hopping” whereby long stints in resonance are disrupted by short bursts of migration in semimajor axis. In addition, if Planet Nine exists with a semimajor axis near the low end of its proposed range, then it can cause the ETNO orbits to be apsidally antialigned (as observed). The apparent alignment of the physical orbits of these objects with that predicted by Planet Nine requires additional work to determine the true extent of Planet Nine's dynamical influence, as the 250 au cutoff may not be sufficient to describe the population of objects shepherded by Planet Nine.

Future discoveries and their associated dynamical studies will expand the census of known ETNOs and allow for a better determination between the various scenarios presented in this work.

This material is based upon work supported by the National Aeronautics and Space Administration under grant NNX17AF21G issued through the SSO Planetary Astronomy Program, and by NSF grant AST-1515015. We thank Konstantin Batygin, Mike Brown, Sarah Millholland, Andrew Vanderburg, John Monnier, and Kat Volk for useful conversations. We thank the anonymous referee for valuable comments that increased the clarity of the manuscript. J.C.B. and S.J.H. are also supported by the NSF Graduate Research Fellowship Grant DGE 1256260. The computations for this work used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation grant ACI-1053575. This research was done using resources provided by the Open Science Grid, which is supported by the National Science Foundation and the U.S. Department of Energy Office of Science.

Funding for the DES projects has been provided by the U.S. Department of Energy, the U.S. National Science Foundation, the Ministry of Science and Education of Spain, the Science and Technology Facilities Council of the United Kingdom, the Higher Education Funding Council for England, the National Center for Supercomputing Applications at the University of Illinois at Urbana-Champaign, the Kavli Institute of Cosmological Physics at the University of Chicago, the Center for Cosmology and Astro-Particle Physics at the Ohio State University, the Mitchell Institute for Fundamental Physics and Astronomy at Texas A&M University, Financiadora de Estudos e Projetos, Fundação Carlos Chagas Filho de Amparo à Pesquisa do Estado do Rio de Janeiro, Conselho Nacional de Desenvolvimento Científico e Tecnológico and the Ministério da Ciência, Tecnologia e Inovação, the Deutsche Forschungsgemeinschaft, and the Collaborating Institutions in the Dark Energy Survey.
The Collaborating Institutions are Argonne National Labora-
tory, the University of California at Santa Cruz, the University
of Cambridge, Centro de Investigaciones Energéticas, Med-
ioambientales y Tecnológicas-Madrid, the University of
Chicago, University College London, the DES-Brazil Con-
sortium, the University of Edinburgh, the Eidgenössische
Technische Hochschule (ETH) Zürich, Fermi National Accel-
erator Laboratory, the University of Illinois at Urbana-
Champaign, the Institut de Ciències de l'Espai (IEEC/CSIC),
the Institut de Física d'Altes Energies, Lawrence Berkeley
National Laboratory, the Ludwig-Maximilians Universität
München and the associated Excellence Cluster Universe, the
University of Michigan, the National Optical Astronomy
Observatory, the University of Nottingham, The Ohio State
University, the University of Pennsylvania, the University of
Portsmouth, SLAC National Accelerator Laboratory, Stanford
University, the University of Sussex, Texas A&M University,
and the OzDES Membership Consortium.

This research is based in part on observations at Cerro
Tololo Inter-American Observatory, National Optical Astron-
omy Observatory, which is operated by the Association of
Universities for Research in Astronomy (AURA) under a
cooperative agreement with the National Science Foundation.

The DES data management system is supported by the National Science Foundation under grants AST-1138766 and AST-1536171. The DES participants from Spanish institutions are partially supported by MINECO under grants AYA2015-71825, ESP2015-66861, FPA2015-68048, SEV-2016-0588, SEV-2016-0597, and MDM-2015-0509, some of which include ERDF funds from the European Union. IFAE is partially funded by the CERCA program of the Generalitat de Catalunya. Research leading to these results has received funding from the European Research Council under the European Union's Seventh Framework Program (FP7/2007-2013) including ERC grant agreements 240672, 291329, and 306478. We acknowledge support from the Australian Research Council Centre of Excellence for All-sky Astrophysics (CAASTRO), through project number CE110001020, and the Brazilian Instituto Nacional de Ciência e Tecnologia (INCT) e-Universe (CNPq grant 465376/2014-2).

This manuscript has been authored by the Fermi Research
Alliance, LLC, under contract DE-AC02-07CH11359 with the
U.S. Department of Energy, Office of Science, Office of High
Energy Physics. The United States Government retains and the
publisher, by accepting the article for publication, acknowl-
edges that the United States Government retains a nonexclu-
sive, paid-up, irrevocable, worldwide license to publish or
reproduce the published form of this manuscript, or allow
others to do so, for United States Government purposes.

Software: WCSfit (Bernstein et al. 2017), mp_ephem
(https://github.com/OSSOS/liborbfit), pandas (McKinney
2010), IPython (Pérez & Granger 2007), matplotlib (Hunter
2007), scipy (Jones et al. 2001), numpy (Oliphant 2006),
Jupyter (Kluyver et al. 2016).

ORCID iDs

T. Khain https://orcid.org/0000-0001-7721-6457

J. C. Becker https://orcid.org/0000-0002-7733-4522

F. C. Adams https://orcid.org/0000-0002-8167-1767

D. W. Gerdes https://orcid.org/0000-0001-6942-2736

S. J. Hamilton https://orcid.org/0000-0002-6126-8487

Hsing Wen Lin (林省文) https://orcid.org/0000-0001-7737-6784

J. Annis https://orcid.org/0000-0002-0609-3987

D. Brooks https://orcid.org/0000-0002-8458-5047

A. Carnero Rosell https://orcid.org/0000-0003-3044-5150

M. Carrasco Kind https://orcid.org/0000-0002-4802-3194

J. García-Bellido https://orcid.org/0000-0002-9370-8360

D. Gruen https://orcid.org/0000-0003-3270-7644

R. A. Gruendl https://orcid.org/0000-0002-4588-6517

D. L. Hollowood https://orcid.org/0000-0002-9369-4157

K. Kuehn https://orcid.org/0000-0003-0120-0808

R. Miquel https://orcid.org/0000-0002-6610-4836

B. Nord https://orcid.org/0000-0001-6706-8972

A. A. Plazas https://orcid.org/0000-0002-2598-0514

A. K. Romer https://orcid.org/0000-0002-9328-879X

E. Sanchez https://orcid.org/0000-0002-9646-8198

M. Smith https://orcid.org/0000-0002-3321-1432

M. Soares-Santos https://orcid.org/0000-0001-6082-8529

F. Sobreira https://orcid.org/0000-0002-7822-0658

M. E. C. Swanson https://orcid.org/0000-0002-1488-8552

G. Tarle https://orcid.org/0000-0003-1704-0781

A. R. Walker https://orcid.org/0000-0002-7123-8943

References

Allen, R. L., Bernstein, G. M., & Malhotra, R. 2001, ApJL, 549, L241

Batygin, K., & Brown, M. E. 2016, AJ, 151, 22

Batygin, K., & Morbidelli, A. 2017, AJ, 154, 229

Becker, J. C., & Adams, F. C. 2017, MNRAS, 468, 549

Becker, J. C., Adams, F. C., Khain, T., Hamilton, S. J., & Gerdes, D. 2017, AJ, 154, 61

Becker, J. C., Khain, T., Hamilton, S. J., et al. 2018, arXiv:1805.05355

Bernstein, G., & Khushalani, B. 2000, AJ, 120, 3323

Bernstein, G. M., Armstrong, R., Plazas, A. A., et al. 2017, PASP, 129, 074503

Bernstein, G. M., Trilling, D. E., Allen, R. L., et al. 2004, AJ, 128, 1364

Bernstein, J. P., Kessler, R., Kuhlmann, S., et al. 2012, ApJ, 753, 152

Beust, H. 2016, A&A, 590, L2

Brown, M. E., Barkume, K. M., Ragozzine, D., & Schaller, E. L. 2007, Natur, 446, 294

Dark Energy Survey Collaboration, Abbott, T., Abdalla, F. B., et al. 2016, MNRAS, 460, 1270

de la Fuente Marcos, C., & de la Fuente Marcos, R. 2017, MNRAS, 471, L61

de la Fuente Marcos, C., & de la Fuente Marcos, R. 2018, MNRAS, 474, 838

de la Fuente Marcos, C., de la Fuente Marcos, R., & Aarseth, S. J. 2017, Ap&SS, 362, 198

de León, J., de la Fuente Marcos, C., & de la Fuente Marcos, R. 2017, MNRAS, 467, L66

Diehl, H. T., Neilsen, E., Gruendl, R. A., et al. 2018, Proc. SPIE, 10704, 107040D

Duncan, M., Quinn, T., & Tremaine, S. 1987, AJ, 94, 1330

Flaugher, B., Diehl, H. T., Honscheid, K., et al. 2015, AJ, 150, 150

Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2016, A&A, 595, A2

Gaudi, B. S. 2004, ApJ, 610, 1199

Gerdes, D. W., Jennings, R. J., Bernstein, G. M., et al. 2016, AJ, 151, 39

Gerdes, D. W., Sako, M., Hamilton, S., et al. 2017, ApJ, 839, L15

Gladman, B., Marsden, B. G., & Vanlaerhoven, C. 2008, in The Solar System Beyond Neptune, ed. M. A. Barucci et al. (Tucson, AZ: Univ. Arizona Press), 43

Goldstein, D. A., D'Andrea, C. B., Fischer, J. A., et al. 2015, AJ, 150, 82

Hadden, S., Li, G., Payne, M. J., & Holman, M. J. 2017, arXiv:1712.06547

He, M. Y., & Petrovich, C. 2018, MNRAS, 474, 20

Hunter, J. D. 2007, CSE, 9, 90

Jones, E., Oliphant, T., Peterson, P., et al. 2001, SciPy: Open source Scientific Tools for Python, https://www.scipy.org/

Kessler, R., Marriner, J., Childress, M., et al. 2015, AJ, 150, 172

Ketchum, J. A., Adams, F. C., & Bloch, A. M. 2013, ApJ, 762, 71

Khain, T., Batygin, K., & Brown, M. E. 2018, AJ, 155, 250

Kluyver, T., Ragan-Kelley, B., Pérez, F., et al. 2016, in Positioning and Power in Academic Publishing: Players, Agents and Agendas, ed. F. Loizides & B. Schmidt (Clifton, VA: IOS Press), 87

Leinhardt, Z. M., Marcus, R. A., & Stewart, S. T. 2010, ApJ, 714, 1789
Levison, H. F., Morbidelli, A., Vokrouhlický, D., & Bottke, W. F. 2008, AJ, 136, 1079

Lin, H.-W., Gerdes, D. W., Hamilton, S. J., et al. 2018, arXiv:1806.09696

Liu, C.-Y., Doressoundiram, A., Roques, F., et al. 2015, MNRAS, 446, 932

Malhotra, R., Volk, K., & Wang, X. 2016, ApJL, 824, L22

Marcus, R. A., Ragozzine, D., Murray-Clay, R. A., & Holman, M. J. 2011, ApJ, 733, 40

McKinney, W. 2010, in Proc. of the 9th Python in Science Conf., ed. S. van der Walt & J. Millman, 51, http://conference.scipy.org/proceedings/scipy2010/pdfs/mckinney.pdf

Millholland, S., & Laughlin, G. 2017, AJ, 153, 91

Morganson, E., Gruendl, R. A., Menanteau, F., et al. 2018, PASP, 130, 074501

Nesvorný, D., Alvarellos, J. L. A., Dones, L., & Levison, H. F. 2003, AJ, 126, 398

Nesvorný, D., Bottke, W. F., Jr., Dones, L., & Levison, H. F. 2002, Natur, 417, 720

Oliphant, T. E. 2006, Guide to NumPy (Provo, UT)

Ortiz, J. L., Thirouin, A., Campo Bagatin, A., et al. 2012, MNRAS, 419, 2315

Pérez, F., & Granger, B. E. 2007, CSE, 9, 21

Petit, J.-M., Kavelaars, J. J., Gladman, B. J., et al. 2017, AJ, 153, 236

Ragozzine, D., & Brown, M. E. 2009, AJ, 137, 4766

Saillenfest, M., Fouchard, M., Tommei, G., & Valsecchi, G. B. 2017, CeMDA, 127, 477

Schaller, E. L., & Brown, M. E. 2008, ApJL, 684, L107

Schlichting, H. E., & Sari, R. 2009, ApJ, 700, 1242

Sheppard, S. S., & Trujillo, C. 2016a, Minor Planet Electronic Circulars, 2016-Q40

Sheppard, S. S., & Trujillo, C. 2016b, AJ, 152, 221

Snodgrass, C., Carry, B., Dumas, C., & Hainaut, O. 2010, A&A, 511, A72

Trujillo, C. A., & Sheppard, S. S. 2014, Natur, 507, 471

Volk, K., & Malhotra, R. 2012, Icar, 221, 106

Volk, K., Murray-Clay, R. A., Gladman, B. J., et al. 2018, AJ, 155, 260

Yu, T. Y. M., Murray-Clay, R., & Volk, K. 2018, arXiv:1805.08228