# High-precision Orbit Fitting and Uncertainty Analysis of (486958) 2014 MU69

Simon Porter, Marc Buie, Alex Parker, John Spencer, Susan Benecchi, Paolo

Tanga, Anne Verbiscer, J. Kavelaars, Stephen Gwyn, Eliot Young, et al.

► To cite this version:

Simon Porter, Marc Buie, Alex Parker, John Spencer, Susan Benecchi, et al. High-precision Orbit Fitting and Uncertainty Analysis of (486958) 2014 MU69. The Astronomical Journal, 2018, 156 (1), pp.20. 10.3847/1538-3881/aac2e1. hal-03526041

HAL Id: hal-03526041

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

Submitted on 20 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished 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
# High-precision Orbit Fitting and Uncertainty Analysis of (486958) 2014 MU69

Simon B. Porter<fnref n="1" />, Marc W. Buie<fnref n="1" />, Alex H. Parker<fnref n="1" />, John R. Spencer<fnref n="1" />, Susan Benecchi<fnref n="2" />, Paolo Tanga<fnref n="3" />, Anne Verbiscer<fnref n="4" />, J. J. Kavelaars<fnref n="5" />, Stephen D. J. Gwyn<fnref n="5" />, Eliot F. Young<fnref n="1" />, H. A. Weaver<fnref n="6" />, Catherine B. Olkin<fnref n="1" />, Joel W. Parker<fnref n="1" />, and S. Alan Stern<fnref n="1" />

Received 2018 April 25; revised 2018 May 1; accepted 2018 May 3; published 2018 June 21

## Abstract

NASA's *New Horizons* spacecraft will conduct a close flyby of the cold-classical Kuiper Belt Object (KBO) designated (486958) 2014 MU69 on 2019 January 1. At a heliocentric distance of 44 au, “MU69” will be the most distant object ever visited by a spacecraft. To enable this flyby, we have developed an extremely high-precision orbit fitting and uncertainty processing pipeline, making maximal use of the *Hubble Space Telescope*'s Wide Field Camera 3 (WFC3) and pre-release versions of the ESA *Gaia* Data Release 2 (DR2) catalog. This pipeline also enabled successful predictions of a stellar occultation by MU69 in 2017 July. We describe how we process the WFC3 images to match the *Gaia* DR2 catalog, extract positional uncertainties for this extremely faint target (typically 140 photons per WFC3 exposure), and translate those uncertainties into probability distribution functions for MU69 at any given time. We also describe how we use these uncertainties to guide *New Horizons*, plan stellar occultations of MU69, and derive MU69's orbital evolution and long-term stability.

**Key words:** astrometry – celestial mechanics – Kuiper Belt: general – Kuiper Belt objects: individual (2014 MU69) – occultations

*Supporting material:* machine-readable table

## 1. Introduction

The cold-classical Kuiper Belt Object (KBO) (486958) 2014 MU69 is the primary target for NASA's *New Horizons* Kuiper Belt Extended Mission. The cold-classical Kuiper Belt consists of objects on low-eccentricity, low-inclination (<5° to the invariant plane) orbits (that is, dynamically “cold”) with heliocentric semimajor axes between about 40 and 50 au. The cold-classical objects were likely formed in-place and escaped perturbation from their initial orbits by giant planet migration (Batygin et al. 2011, and references therein), making them the most distant known remnants of the original protoplanetary disk.

NASA's *New Horizons* spacecraft was launched 2006 January 19, received a gravitational assist from Jupiter on 2007 February 28, and flew through the Pluto-Charon system on 2015 July 14 (Stern et al. 2015). Since the Pluto encounter, *New Horizons* has observed the 3:2 Neptune resonant (15810) Arawn (provisionally designated 1994 JR1) in 2016 as close at a distance of 0.7 au (Porter et al. 2016). *New Horizons* will encounter many other KBOs within 1 au, some as close as 0.1 au, and some (such as Quaoar and Haumea) that are much farther away, but all can be seen by *New Horizons* at much higher solar phase angles than is possible from Earth-based telescopes (S. B. Porter et al. 2018, in preparation, A. J. Verbiscer et al. 2018, in preparation). However, none of these KBOs will be seen as close as MU69, which will be within 3500 km of the spacecraft on the nominal trajectory. *New Horizons* will image the surface of MU69 at best resolutions of ≈35 m/pixel, and spectral maps at ≈1 km/pixel. In order to guide the spacecraft to such a close encounter with a KBO that only has a relatively short orbital

arc required a completely new approach to orbit determination and uncertainty analysis, which we describe in this paper.

2014 MU69 was discovered in 2014 July by the *Hubble Space Telescope (HST)* following eight years of dedicated searches for a second *New Horizons* encounter object (M. W. Buie et al. 2018, in preparation). After several ground-based searches down to $V \approx 26$, 194 *HST* orbits were allocated for a deeper, more systematic search for objects accessible to *New Horizons* (GO 13633, PI Spencer). MU69 was initially detected in 10 images acquired in two *HST* orbits, as were four other KBOs during the *HST* search. Three objects were potential targets for *New Horizons*: 2014 MU69, 2014 OS393, and 2014 PN70. In 2015 August, the *New Horizons* team selected 2014 MU69 as the potential *New Horizons* extended mission target. The spacecraft performed a series of four burns in 2015 October–November to target 2014 MU69. The *New Horizons* Kuiper Belt Extended Mission was approved by NASA after Senior Review in 2016 July, and its centerpiece is the flyby of 2014 MU69 on 2019 January 1.

In this paper, we will discuss our process of performing absolute astrometry on 2014 MU69 tied to a pre-release version of *Gaia* DR2, propagating that error forward to orbital uncertainty, and then using the orbital uncertainty to guide both occultations of the KBO and to guide the spacecraft to a close flyby. These techniques represent the highest-precision heliocentric orbit fitting of a KBO ever, and they can provide a basis for future applications of *Gaia*-driven astrometry to small bodies in the solar system.

## 2. Data Sources

2014 MU69 was discovered with the *HST* search program (GO 13633) described in (M. W. Buie et al. 2018, in

<footnote n="1">¹Southwest Research Institute, 1050 Walnut Street, Suite 300, Boulder, CO 80302, USA; porter@boulder.swri.edu</footnote>

<footnote n="2">²Planetary Science Institute, 1700 East Fort Lowell, Suite 106, Tucson, AZ 85719, USA</footnote>

<footnote n="3">³Université Côte d'Azur, Observatoire de la Côte d'Azur, CNRS, Laboratoire Lagrange, Bd de l'Observatoire, CS 34229, F-06304 Nice cedex 4, France</footnote>

<footnote n="4">⁴University of Virginia, P.O. Box 400325, Charlottesville, VA 22904, USA</footnote>

<footnote n="5">⁵National Research Council of Canada, Victoria, BC, Canada</footnote>

<footnote n="6">⁶Johns Hopkins University Applied Physics Laboratory, Laurel, MD, USA</footnote>
preparation). This program was designed to take five full-frame
370 s Wide Field Camera 3 (WFC3) UVIS images in one orbit
with the F350LP broadband filter, skip an *HST* orbit, and then
repeat the same observation. The images were tracked on a
nominal cold-classical orbit, producing streaked stars, but
nonstreaked KBOs. To perform the search, the images were
shift-stacked at 20 different representative cold-classical shift
rates. Objects that appeared in both search orbits for one of the
shift rates were identified and targeted for follow up *HST*
observations. The first object identified this way was
designated 1110113Y (*HST* orbit IDs 11/12, WFC3 CCD 1,
shift rate ID 011, random ID 3Y), later given the provisional
designation 2014 MU<sub>69</sub>. Four more KBOs were subsequently
detected, all of which were brighter, but all of which required
more fuel for *New Horizons* to reach.

The available data set for 2014 MU<sub>69</sub> has both a short temporal
arc (2014 July–2017 October) and an extremely high data quality,
making it ideal for the analysis described below. MU<sub>69</sub> is a very
faint object, with V ≈ 27.5 (S. Benecchi et al. 2018, in
preparation) and is in an extremely crowded star field (galactic
longitudes from −8° to −12°). These constraints have made it
effectively impossible to detect with ground-based telescopes, and
all observations of MU<sub>69</sub> have been with the *Hubble Space
Telescope*. A list of these *HST* observations is in Table 3.

Both the initial follow up observations and most observa-
tions conducted since then have adopted the search program's
basic format of five 367 to 370 s, F350LP filter, full-frame
UVIS WFC3 images. The key exceptions are the color
campaign in summer of 2016 (GO 14092, PI Benecchi), in
which four orbits each included two 348 s images using F606W
followed by three 373 s images with F814W. As shown in
Table 3, half of the follow up orbits were roughly evenly
spread over 2014 August–2017 October, while the other half
were spread over a roughly one week interval in 2017 June–
July. The latter was the light-curve campaign (GO 14627, PI
Benecchi), which was critical in successfully predicting the
2017 July occultation (see Section 6).

In addition to *HST* and the 2017 July occultation, the other data
source for this analysis was stellar astrometry from the ESA *Gaia*
project. Initially, the process in Section 3 was built using a custom
star catalog built from a deep composite of the MU<sub>69</sub> field
obtained using the Canada-France-Hawaii Telescope (Gwyn
2014). When the *Gaia* Data Release 1 (DR1) was made available
in 2016 September (Gaia Collaboration et al. 2016), we began
using that, applying a mean proper motion correction to images
significantly after the DR1 2015.0 epoch. The mean proper
motion was calculated from the *Gaia* TGAS catalog (Michalik
et al. 2015). The MU<sub>69</sub> observation fields were in areas of very
low coverage for TGAS, so the TGAS stars could not be used
directly. With the success of our application of DR1 and after a
special support request to the *Gaia* project, we were able to obtain
a sky patch from a pre-release version of Data Release 2 (DR2)
around the path of MU<sub>69</sub>. The major advances with DR2 are
proper motion for all catalog stars, obviating the need for a mean
proper motion correction, as well as a much more homogeneous,
bias-free distribution of errors on the sky, and much lower
uncertainty (one order of magnitude). This early version of DR2
was used to plan the three 2017 occultations, both for correcting
the *HST* absolute astrometry and for knowledge of the occultation
stars themselves. We obtained a second preview version of DR2
in our field of interest in 2017 October, and that version is used for
the astrometry in Table 3.

**3. Image Analysis**

A typical *HST* MU<sub>69</sub> image is shown in Figure 1. Over the course of 2014–2017, MU<sub>69</sub> has moved between galactic latitudes of −8° and −12°. Accordingly, the background star density has always been very high in any images of MU<sub>69</sub> (from both Earth and *New Horizons*), and this dense background field will persist through the *New Horizons* encounter. To mitigate this background, we developed a Python program called *warpy.py* to perform simple star subtraction. Because almost every observation sequence consists of five images of roughly the same field, *warpy.py* iterates through the images, warps the other four images of the visit to the frame of the fifth, median combines the four warped images, and subtracts them from the target image. The images are coregistered by matching sources detected with *Source Extractor* (Bertin & Arnouts 1996). Because the stars in the images are all smeared differently, they do not subtract cleanly. However, the star-subtracted images are useful for verification that the following steps are fitting the KBO and not a background source.

The smeared stars are also a problem for determining the
pointing of the images. As *HST* is tracking on the motion of the
object, a single Tiny Tim (Krist et al. 2011) point-spread
function (PSF) would not accurately describe the effective
PSFs of the stars. We thus built up “smear kernels” that
describe the motion of the stars relative to the KBO through the
370 s exposures. By shifting a Tiny Tim PSF to 400 discrete
times during the exposure and averaging them, we were able to
build up an exact model of each star’s PSF. We did this for
each star used for the World Coordinate System (WCS)
solution (effectively, all the *Gaia* stars in the field), as the Tiny
Tim PSF varies across the WFC3 field.

We next used the stellar PSFs to build up probability
distribution functions (PDFs) for the pixel location of each star
within each image. We did this task with a Markov-Chain
Monte Carlo (MCMC) algorithm implemented in IDL. This
minimal MCMC with a single “walker” was iterated for 2000
steps to build a PDF of both pixel position and total DN/
second flux. The result is 2000 equal-probability pixel positions
for the star, encompassing the true shape of the uncertainty
distribution. We found this number of steps to be sufficient, as
all the stars used in the astrometric solution had a high signal-
to-noise ratio and typically there were 70–100 stars used in a
given solution. The flux number was not used directly in the
fits, but provided diagnostics if the fits were successful. A
typical star had a 1σ pixel position uncertainty of <0.1 pixels,
equating to an angular uncertainty of <4 mas.

With these stellar PDFs, we could now build PDFs for the
pointing of *HST* in each image. WFC3 UVIS images use a
three-layer WCS to translate pixel to sky coordinates, as
described in Greisen et al. (2006). The first layer is the basic
pointing, roll, and trapezoidal warp, the second is a set of SIP
polynomials that describe the low-frequency distortion on the
chip, and the third layer is a look-up table that describes high-
frequency pixel distortion (e.g., by irregularities in the
lithography of the CCD). All of the distortion parameters are
highly calibrated, and we only needed to update the pointing
with deltas to the *CRVAL1* and *CRVAL2* keywords, and the roll
by multiplying the CD matrix by a rotation matrix. Thus, for
each image, we have three parameters for the WCS PDF: delta
R.A., delta decl., and delta roll. We built this WCS PDF by
selecting a pixel coordinate for each star from their PDFs,
fitting a best-fit WCS solution, and recording the deltas from
<figure><img src="image_2.png" /><figcaption>**Figure 1.** An image from the *HST* light-curve campaign (id8i16moq). The blue circles show the *Gaia* DR2 positions of stars, the red circles show the locations of star PDFs used for the WCS solution, and the outer yellow circles indicate a match (all PDF stars were matched in this case, which is typical). The inset green box shows the location of 2014 MU<sub>69</sub>; see Figure 2.</figcaption></figure>

the original WCS. This process was repeated 10,000 times to build a discretely sampled PDF of the WCS offsets. The resulting typical $1\sigma$ uncertainty in the *CRVAL1* and *CRVAL2* keywords (and thus in the pointing of *HST*) was <2 mas.

The pixel position PDFs for the KBO used the same basic MCMC algorithm as the star pixel PDFs, but with the single walker iterated 10000 times to build the PDF. Because *HST* tracked on the KBO, we did not need a smear kernel to fit the KBO and could use a Tiny Tim PSF directly. In addition, all but the discovery observations had MU<sub>69</sub> near the center of WFC3 chip 2 (FITS extension 1), making for a less distorted PSF than near the chip edges. Initially, the KBO fit was started with a manual click on the rough position of the KBO in the image. However, as the orbit improved, this manual position was replaced with a calculated initial position from prior orbit solutions and the WCS. We also made manual masks of stars and cosmic rays near the KBO that might adversely affect the PDF generation.

Finally, we needed to combine the KBO pixel PDFs with the WCS PDFs to make KBO sky PDFs. Similarly to the WCS PDFs, this was accomplished by selecting a randomly selected KBO pixel location with a randomly selected WCS PDF, translating to R.A. and decl., and then repeating 10000 times. An example of one of these PDFs is shown in Figure 2. The instrument magnitude in the pixel PDF was converted to AB apparent magnitude with the PHOTFLAM header keyword. While the magnitude was not used directly for astrometry, it was an important diagnostic for the quality of the pixel PDF. For the solutions presented, we filtered out any points with magnitude uncertainties larger than 0.5, which were generally failures to fit the object, or any points with uncertainties smaller than 0.1 mag, typically a cosmic ray close to the KBO causing spuriously high signal-to-noise ratios. We also only used the F350LP points for the orbit, as the narrower-band points had worse signal-to-noise ratios for both stars and MU<sub>69</sub>. This left 214 of the 264 images, which we used for an initial fit. An additional nine points were rejected because they had greater than 30 mas residuals relative to the initial fit. Our final *HST* astrometry thus used 205 of the 264 images (78%); these are shown in Table 3.

After the method described here was used to successfully predict the occultation of MU<sub>69</sub> on 2017 July 17 (See Section 6, M. W. Buie et al. 2018, in preparation), we were able to use the occultation itself as a high-quality occultation point. Because five solid-body chords were obtained on

<figure><img src="image_3.png" /><figcaption>**Figure 2.** The R.A./decl./Magnitude Probability Distribution Function (PDF) for MU<sub>69</sub> in the image id8i16moq (see Figure 1). R.A. and decl. are in milliarcseconds, so as to show the uncertainty in appropriate units. We generated similar PDFs for all of the astrometry shown in Table 3.</figcaption></figure>

July 17, we chose the mid-time of the longest chord and used it as the nominal center-of-figure. We could then combine this mid-time, the topocentric location of the portable telescope that obtained the longest chord, and the location and uncertainty of the occultation star from *Gaia* DR2 to produce an effective astrometric PDF. See M. W. Buie et al. (2018, in preparation) for more details about the circumstances and analysis of this occultation. This occultation PDF could then be combined with the *HST*-derived PDFs in the process described in Section 4.

## 4. Orbit Determination

Typically, small body orbits in the literature are described in either mean or osculating heliocentric elements, with error bars representing a normal error distribution. This is typically sufficient
for general dynamical studies and rough targeting from the ground, but not for spacecraft flybys or occultation planning. The actual uncertainty of an object's astrometry is rarely perfectly described by a normal distribution, and neither is that object's location and velocity in space. We thus sought to develop an orbit-fitting method that would accurately map the full astrometric uncertainty distribution into the ephemeris.

To perform these fits, we developed a high-precision few-body orbital integrator. As 2014 MU<sub>69</sub> is a cold-classical KBO, all of the planetary perturbations on it are interior, and tend to result in a slow precession. Nongravitational factors (i.e., YORP) and general relativity are not a factor for KBOs. We therefore developed a few-body conservative force integrator, capable of modeling the major planets and their perturbing forces on a massless test particle. This integrator (PyNBody<sup>7</sup>) is based on the 12/13 order Runge-Kutta-Nystrom integrator of Brankin et al. (1989) and was previously described in Porter et al. (2016). This integrator is not the fastest, but it is very accurate and can typically conserve system energy and momentum to within machine precision over the relevant timescales for orbit fitting.

The KBO's orbit is parameterized as a Cartesian state vector relative to the solar system barycenter at a fixed epoch. The inertial frame for the integrations the International Celestial Reference Frame (ICRF). Gaia DR1 is aligned to ICRF by matching optical detections of quasars with a subset of ICRF2, while Gaia DR2 uses several thousand quasars from ICRF3 and half a million AGNs to perform frame alignment (Gaia Collaboration et al. 2016). The integration epoch is set to 2014-06-01 00:00:00.000 UTC, a few weeks before the first observation (originally a safety factor in case of any precoveries in the HST search). To test the solution against the data, we propagate the state vector with the PyNBody integrator to the desired time and calculate its apparent ICRF R.A./decl. from HST, with appropriate light-time correction. We use the JPL NAIF HST and DE430 SPICE kernels to determine the location of HST relative to the solar system barycenter (Folkner et al. 2014).

We used the *emcee* Markov-Chain Monte Carlo package (Foreman-Mackey et al. 2013) to translate the astrometric uncertainty to orbital uncertainty. The *emcee* package provides a fast and natively multithreaded way to run MCMC from Python. As input to *emcee*, the fitting program calculates the likelihood for any solution by taking the predicted R.A./decl. for that solution and comparing them to the R.A./decl. PDFs. Because the PDFs are discretely sampled, we created a Kernel Density Estimator for each observations, using Silverman's Rule of Thumb (Silverman 1986) to choose the bandwidths, as most of the PDFs were roughly Gaussian. The log likelihoods for all the images could then be summed to provide a total log likelihood to *emcee*.

For any solution, the first step is to make an initial guess (typically an older solution) and minimize its $\chi^2$ with a downhill simplex method (Nelder & Mead 1965). This polished solution is then used to create 200 slightly perturbed state vectors as the initial “walkers” for *emcee* to use. We then run the 200 walkers for 100 iterations to “burn-in” and allow them to move away from the artificial initial distribution. We then reset *emcee* and run it for 500 iterations to produce the full PDF cloud of 10,000 state vectors at the fitting epoch. These numbers of iterations were arrived at after much testing and are typically more burn-in than is actually necessary, so as to ensure that the solutions are well-distributed. We save the

<table><caption>**Table 1**
The “rd2b” Orbit Solution for 2014 MU<sub>69</sub></caption><thead><tr><th></th><th>Value</th><th></th><th>1σ</th><th></th></tr></thead><tbody><tr><td>x</td><td>+1.163133074444e+09 ± 2.80233e+02</td><td>km</td><td></td><td></td></tr><tr><td>y</td><td>-6.385039581373e+09 ± 1.52754e+03</td><td>km</td><td></td><td></td></tr><tr><td>z</td><td>+2.373261916929e+08 ± 5.87015e+01</td><td>km</td><td></td><td></td></tr><tr><td>v<sub>x</sub></td><td>+4.461378977476e+00 ± 5.92714e-06</td><td>km s<sup>-1</sup></td><td></td><td></td></tr><tr><td>v<sub>y</sub></td><td>+9.619622770583e-01 ± 2.45488e-05</td><td>km s<sup>-1</sup></td><td></td><td></td></tr><tr><td>v<sub>z</sub></td><td>-1.066958207821e-01 ± 9.45150e-07</td><td>km s<sup>-1</sup></td><td></td><td></td></tr><tr><td>a</td><td>44.23555350 ± 0.00003999</td><td>au</td><td></td><td></td></tr><tr><td>e</td><td>0.03787388 ± 0.00000476</td><td></td><td></td><td></td></tr><tr><td>I</td><td>2.44993086 ± 0.00000203</td><td>deg</td><td></td><td></td></tr><tr><td>Ω</td><td>159.04712465 ± 0.00006746</td><td>deg</td><td></td><td></td></tr><tr><td>ω</td><td>183.74800591 ± 0.00469779</td><td>deg</td><td></td><td></td></tr><tr><td>M</td><td>301.30454775 ± 0.00406885</td><td>deg</td><td></td><td></td></tr></tbody></table>

*Note.* State vector and orbit are relative to the solar system barycenter and in the ICRF ecliptic frame at the epoch 2014-06-01 00:00:00 UTC.

resulting state vector PDF in a format that can then be propagated to any time of interest. The state vector and orbit for our “rd2b” orbit solution are presented in Table 1.

While the full PDF cloud of 10,000 states encapsulates the uncertainty in the location of MU<sub>69</sub> any any time, it is rather unwieldy to use in most circumstances. We therefore calculated the state clouds of 2014 MU<sub>69</sub> at 2000 discrete times between 2004 January 1 and 2024 January 1, and averaged the states at each time. We then generated order-27 Chebyshev polynomials (Tchebychev 1853) for the positions of the KBO, and saved them in a JPL SPICE Type 02 SPK kernel. JPL SPICE could then be used to rapidly interpolate the location of MU<sub>69</sub> at any time along the interval. Testing the kernel at 769 random points over the interval returned a root-mean-squared residual of 20 m, well below the uncertainty in the orbit.

## 5. New Horizons Trajectory Planning

The primary reason to determine the orbit of 2014 MU<sub>69</sub> to very high precision is to ensure the success of the *New Horizons* flyby. *New Horizons* performed the major Trajectory Correction Maneuver to guide it to MU<sub>69</sub> over a series of four burn segments in 2015 October and November, after all Pluto observations had finished. The initial orbit used to target the spacecraft was based on the first year of data, from 2014 June to 2015 July (GO 13633 and GO/DD 14053, PI Spencer). After that burn, early versions of the analysis described here showed that significantly more *HST* observations would be required to enable a close flyby of 2014 MU<sub>69</sub>. We thus proposed and were awarded six *HST* orbits in 2016, and five in 2017 (GO/DD 14485, GO 14629, and GO 15158, PI Buie). In addition, 24 *HST* orbits were used in 2017 June/July to measure the light curve of MU<sub>69</sub> (GO 14627, PI Benecchi). The orbit presented here uses data from all of these *HST* programs, in addition to the July 17 occultation.

The *New Horizons* spacecraft will nominally fly closest to 2014 MU<sub>69</sub> at 05:33 2019 January 1 UTC. This time was chosen to enable both the Goldstone and Canberra Deep Space Network (DSN) 70 m dishes to uplink to the spacecraft simultaneously for an attempted bistatic radar experiment (as was performed at Pluto; Linscott et al. 2016). *New Horizons* will not be able to acquire MU<sub>69</sub> any earlier than 2018 August. Because of *New Horizons*'s almost radial trajectory out of the solar system, the KBO will move very slowly against the background stars until a just few weeks before encounter.

<footnote n="7">7 https://github.com/ascendingnode/PyNBody</footnote>
While the spacecraft can use the LOnG Range Reconnaissance Imager (LORRI; Conard et al. 2005) to well constrain the location of 2014 MU<sub>69</sub> in the “B-Plane” (the plane perpendicular to the spacecraft’s motion and containing the flyby target), the time-of-flight (ToF) uncertainty along the direction of the spacecraft’s motion is constrained only by the Earth-based orbital solution. The distance to the spacecraft from Earth will be well-constrained by Doppler radio measurements on approach to MU<sub>69</sub>, and so the uncertainty in the absolute location of MU<sub>69</sub> relative to the solar system barycenter will determine the flyby ToF uncertainty.

**6. Application to Occultation Planning**

In addition to guiding *New Horizons*, we also used our orbit solution to predict three stellar occultations by 2014 MU<sub>69</sub> in 2017 and one in 2018. The 2017 occultation campaign is comprehensively described in M. W. Buie et al. (2018, in preparation) and here we detail only the procedures used to predict the occultations. MU<sub>69</sub> is a small object, with an absolute magnitude of H<sub>V</sub> ≈ 11 S. Benecchi et al. (2018, in preparation), corresponding to a size likely smaller than 50 km diameter. We therefore knew that the occultations would only be successful if we had very high-quality orbital estimates and uncertainty models. Thankfully, that is exactly what we had developed for guiding *New Horizons*.

Stellar occultations occur when a solar system object passes in front of a star from the perspective of an observer. They have been used to discover the atmosphere of Pluto (Elliot et al. 1989) and the rings around Uranus, Chariklo, and Haumea (Elliot et al. 1977; Braga-Ribas et al. 2014; Ortiz et al. 2017, respectively). The latter is most important for planning the *New Horizons* flyby of MU<sub>69</sub>, as occultations provide the only way of detecting rings or other opacity structures around the KBO before the spacecraft is close enough to see them directly. In addition, occultations can (and in this case did) provide estimates of the size and shape of a body. Knowledge of the approximate size of MU<sub>69</sub> enabled estimates of its bulk albedo, and therefore allowed mission planners to better estimate the correct exposure times for the flyby images.

Because of the motion and rotation of the Earth, stellar occultations sweep across Earth from west to east. The typical approach to observe an occultation of an object with an uncertain orbit is therefore to set up a north-south “picket fence” of portable telescope teams perpendicular or “crosstrack” to the occultation path. It is therefore important to know both the crosstrack uncertainty of the prediction, and what the crosstrack distance of any station is. The latter often requires some iteration, as finding a logistically viable site with the proper crosstrack can be challenging, especially in an unfamiliar country. We thus developed tools to export KML files to Google Maps with lines showing the target crosstracks for each observing team. These could be used for planning site reconnaissance, and with GPS-enabled smartphones, used to see in real time where a potential site was located compared to the desired crosstrack line. To estimate the crosstrack uncertainty, we propagated the full 10000 states to the occultation time and calculated the geometry for all of those states. This produced a full PDF of the occultation uncertainty, which we could use to plan the crosstracks to maximize the chance of success.

The first MU<sub>69</sub> occultation of 2017 was on June 3. The ground track for this event passed over both South America and South Africa, and 25 portable telescope teams were deployed to both Mendoza, Argentina and the Northern Cape and Western Cape Provinces of South Africa. Two *HST* astrometric observations of

MU<sub>69</sub> had been planned for the spring of 2017, on March 16 and in mid May. However, the March observation failed due to a safing event on *HST*, and the observation could not be rescheduled in April, because MU<sub>69</sub> was passing through quadrature and did not have enough sky motion from *HST* to be detected. Thus, the first MU<sub>69</sub> observations of 2017 were acquired by *HST* on May 1 and 25. An initial solution though May 1, “may1c,” was used to plan the deployments to Argentina and South Africa. This was superseded by the “may25a” solution with data through May 25, which was produced after the orbit fitters (Porter and Buie) had deployed, and was used to plan the actual ground tracks. The “may25a” solution purported to have a crosstrack uncertainty of 44 km, though the subsequent 2017 June–July observations showed that the “may25a” ground track had been too far north by almost 2-sigma. This offset precluded a solid-body occultation on June 3, though the high signal-to-noise ratio observations of the event at the South African Astronomical Observatory 74 inch telescope and at Gemini South did exclude optically thick rings around MU<sub>69</sub>. See M. W. Buie et al. (2018, in preparation) for more details.

The second MU<sub>69</sub> occultation of 2017 was on July 10. This event occurred mainly over the Pacific Ocean with a much dimmer star (V ≈ 15.5) and nearly full Moon, thus preventing a ground-based observation campaign. However, the NASA-DLR Stratospheric Observatory for Infrared Astronomy (SOFIA) airborne observatory was able to reach the occultation track from its southern deployment base in Christchurch, New Zealand, and NASA awarded a flight to observe the July 10 occultation (PI E. Young). Between the June 3 and July 10 events, *HST* observed 2014 MU<sub>69</sub> over 24 orbits between June 25 and July 4 (GO 14627, PI Benecchi). This program provided a wealth of new images to integrate into the MU<sub>69</sub> orbit solution in a very short amount of time, and time was especially critical, as the last six orbits worth of data was downlinked from *HST* after the orbit fitting team (Porter and Buie) had arrived in Christchurch. Thus, the orbit solution used to guide the SOFIA flight was necessarily determined at the United States Antarctic Program Christchurch facility, and delivered to the SOFIA mission planners 36 hours before the flight. This orbit solution, “lc1,” used all of the light-curve campaign points, plus the highest-quality preceding *HST* MU<sub>69</sub> observations.

The final MU<sub>69</sub> occultation of 2017 on July 17 was observed with portable ground stations in the Chubut and Santa Cruz provinces in Argentina. No additional observations of *HST* MU<sub>69</sub> were made between the July 10 and July 17 occultations, but we did perform a more thorough filtering of low-quality points. This resulted in the “lc1gr” solution that was used to guide the placement of stations for the July 17 occultation. The “lc1gr” solution had a 1σ uncertainty in crosstrack of 13 km, much tighter than for June 3. This solution allowed a tighter picket fence of stations up and down the Patagonian coast, centered a few kilometers north of the city of Comodoro Rivadavia. Despite high winds on the occultation night, 22 of the 24 deployed stations successfully observed the occultation. Five of those stations observed the solid-body occultation, with the southernmost being the predicted centerline from the “lc1gr” solution.

This work predicts an additional stellar occultation opportunity on 2018 August 4. The ground track for this event passes over western Africa (Mail, Mauritania, and Senegal) and northern South America (Guyana, Venezuela, and Colombia). With the “rd2b” solution presented in Table 1, the 1σ crosstrack uncertainty is 12 km. This uncertainty will decrease somewhat with additional *HST* observations in 2018.
## 7. Long-term Orbital Evolution

2014 MU<sub>69</sub> is a cold-classical KBO. Elliot et al. (2005) identified the “Classical” KBOs as nonresonant objects with eccentricities smaller than 0.2. This classification was further refined as a “Cold Classical” or “Kernel” population by Petit et al. (2011) with *a* ≈ 44 au, *e* ≈ 0.05, *i* < 5° to the invariant plane. MU<sub>69</sub> has *a* = 44.2 au, *e* = 0.03, *i* = 2°4, making it an phenotype of the cold-classical population. Batygin et al. (2011) showed the orbits of cold-classical objects were likely formed in-place and survived being disturbed from their initial orbits by giant planet migration. The unusually high binary fraction of cold-classical KBOs (Noll et al. 2008) is an additional line of evidence that they are mostly undisturbed from their original orbits. Indeed, the observed cold-classical KBO binary fraction is high enough that nearly all must have originally formed as binaries or higher-order multiple systems (Fraser et al. 2017).

With a few small modifications to the PyNBody code, we were able to integrate the orbit of MU<sub>69</sub> over sufficiently long timescales to test this stability and determine mean orbital elements. Specifically, we changed the time unit in the integration from seconds to years to allow for longer integrations without worry of overflows and removed the terrestrial planets as perturbers (instead dropping their masses into the Sun). The results of integrations forward and back 10<sup>8</sup> years can be seen in Figure 3 and Table 2, projected in both the mean solar system plane defined by the *de430*.bsp planets in the ICRF2000 Ecliptic frame, *i*<sub>m</sub> = 1°6 and Ω<sub>m</sub> = 72°4, and in the mean Kuiper Belt plane at 44 au as determined from known KBOs Volk & Malhotra (2017), *i*<sub>m</sub> = 1°8 and Ω<sub>m</sub> = 77°0. The mean, free, and forced elements of MU<sub>69</sub>'s orbit are shown in Table 2. The forced inclination of MU<sub>69</sub> to the mean solar system plane is 0°26, but only 0°0012 to the mean Kuiper Belt at 44 au. Likewise, the forced eccentricity of MU<sub>69</sub> is less than 0.0001. The apparent lack of any forced inclination or eccentricity to the mean Kuiper Belt is strong evidence that MU<sub>69</sub> has not suffered any significant orbital evolution beyond secular perturbations. MU<sub>69</sub> should therefore represent a truly pristine fossil of the Sun's protoplanetary disk, an object unlike any other previously visited by a spacecraft.

## 8. Summary

We have described the process we have used to fit the orbit of 2014 MU<sub>69</sub>, as of the start of 2018. This process combines *Gaia* DR2 and HST/WFC3 to produce extremely high-precision absolute astrometry of MU<sub>69</sub>, and translates that uncertainty into a Cartesian state vector probability distribution function that can be evolved to any time of interest. The results of this analysis were used to successfully predict and observe a solid-body stellar occultation of MU<sub>69</sub> on 2017 July 17, predict a stellar occultation on 2018 August 4, and to guide the *New Horizons* spacecraft to a close (3500 km) flyby of MU<sub>69</sub> on 2019 January 1.

The process described here should enable high-precision orbit determination for future occultations and spacecraft missions. 2014 MU<sub>69</sub> presents the extreme case of a very interesting object that is both faint and in a very crowded star field. Now that the *Gaia* DR2 catalog has been released, solar system objects with higher signal-to-noise ratios should benefit even more from this technique, enabling a substantial improvement in orbital uncertainty and increasing the number of objects that might be observed with stellar occultations.

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Free and forced elements for best-fit 2014 MU<sub>69</sub>; the centers of the circles show forced inclination/eccentricity, while the radii show free inclination/eccentricity. The forced inclination is 0°26 from the mean solar system angular momentum vector (left), but almost perfectly fits the mean Kuiper Belt at 44 au pole from (Volk & Malhotra 2017, right). The lack of forced eccentricity or inclination implies that MU<sub>69</sub> has not experienced any significant orbital evolution since formation.</figcaption></figure>

<table><caption>**Table 2**
Mean, Free, and Forced Elements for Best-fit MU<sub>69</sub></caption><thead><tr><th>Best Fit</th><th>10<sup>8</sup> years</th></tr></thead><tbody><tr><td>Mean <i>a</i></td><td>44.23 au</td></tr><tr><td>Forced <i>e</i></td><td>6 × 10<sup>-5</sup></td></tr><tr><td>Free <i>e</i></td><td>0.037</td></tr><tr><td>Forced <i>i</i><sub>mss</sub></td><td>0°26</td></tr><tr><td>Forced <i>i</i><sub>mkb</sub></td><td>0°0012</td></tr><tr><td>Free <i>i</i></td><td>2°54</td></tr></tbody></table>

This work was supported by NASA's *New Horizons* mission and HST programs GO 13633, GO/DD 14053, GO 14092, GO/DD 14485, GO 14627, GO 14629, and GO 15158. Support for this program was provided by NASA through a grant from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555. Special thanks to ESA for providing pre-release versions of *Gaia* DR2 over the relevant regions. Special thanks to Bill Folkner for providing his expertise and verification of the orbit fitting results.

**Facilities:** HST(WFC3), *Gaia*, SOFIA.

**Software:** astropy (The Astropy Collaboration et al. 2018), scipy (Jones et al. 2001), emcee (Foreman-Mackey et al. 2013), PyNBody (Porter et al. 2016), Matplotlib (Hunter 2007), photutils (Bradley et al. 2017), SpiceyPy (Annex et al. 2018).

## Appendix

The Appendix comprises Table 3.
<table><caption>**Table 3**
HST/WFC3 Astrometry for MU<sub>69</sub></caption><thead><tr><th>WFC3<br>Data Set</th><th>UTC Time</th><th>R.A.</th><th>Decl.</th><th>σ R.A.<br>(mas)</th><th>σ decl.<br>(mas)</th></tr></thead><tbody><tr><td>icii1r7q</td><td>2014 Jun 26T08:51:42.4042</td><td>18h45m10s67179</td><td>-20°53'03"0565</td><td>6.496</td><td>5.874</td></tr><tr><td>icii1r8q</td><td>2014 Jun 26T09:00:33.4117</td><td>18h45m10s63684</td><td>-20°53'03"1070</td><td>8.280</td><td>8.519</td></tr><tr><td>icii1raq</td><td>2014 Jun 26T09:09:24.4036</td><td>18h45m10s59521</td><td>-20°53'03"1753</td><td>9.150</td><td>9.504</td></tr><tr><td>icii1rcq</td><td>2014 Jun 26T09:18:15.3964</td><td>18h45m10s55800</td><td>-20°53'03"2431</td><td>7.130</td><td>6.348</td></tr><tr><td>icii1req</td><td>2014 Jun 26T09:27:06.4048</td><td>18h45m10s52223</td><td>-20°53'03"3121</td><td>6.935</td><td>5.378</td></tr></tbody></table>

(This table is available in its entirety in machine-readable form.)

## ORCID iDs

Simon B. Porter https://orcid.org/0000-0003-0333-6055

Marc W. Buie https://orcid.org/0000-0003-0854-745X

Alex H. Parker https://orcid.org/0000-0002-6722-0994

Susan Benecchi https://orcid.org/0000-0001-8821-5927

Anne Verbiscer https://orcid.org/0000-0002-3323-9304

J. J. Kavelaars https://orcid.org/0000-0001-7032-5255

H. A. Weaver https://orcid.org/0000-0003-0951-7762

Catherine B. Olkin https://orcid.org/0000-0002-5846-716X

## References

Annex, A., Carcich, B., Badger, T. G., et al. 2018, AndrewAnnex/SpiceyPy: SpiceyPy, v2.1.1, Zenodo, doi:10.5281/zenodo.1228513, https://github.com/AndrewAnnex/SpiceyPy

Batygin, K., Brown, M. E., & Fraser, W. C. 2011, ApJ, 738, 13

Bertin, E., & Arnouts, S. 1996, A&AS, 117, 393

Bradley, L., Sipocz, B., Robitaille, T., et al. 2017, astropy/photutils, v0.4, Zenodo, doi:10.5281/zenodo.1039309, https://photutils.readthedocs.io/

Braga-Ribas, F., Sicardy, B., Ortiz, J. L., et al. 2014, Natur, 508, 72

Brankin, R. W., Gladwell, I., Dormand, J. R., Prince, P. J., & Seward, W. L. 1989, ACM Transactions on Mathematical Software, 15, 31

Conard, S. J., Azad, F., Boldt, J. D., et al. 2005, Proc. SPIE, 5906, 407

Elliot, J. L., Dunham, E., & Mink, D. 1977, Natur, 267, 328

Elliot, J. L., Dunham, E. W., Bosh, A. S., et al. 1989, Icar, 77, 148

Elliot, J. L., Kern, S. D., Clancy, K. B., et al. 2005, AJ, 129, 1117

Folkner, W. M., Williams, J. G., Boggs, D. H., Park, R. S., & Kuchynka, P. 2014, IPNPR, 196, 1

Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306

Fraser, W. C., Bannister, M. T., Pike, R. E., et al. 2017, NatAs, 1, 0088

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

Greisen, E. W., Calabretta, M. R., Valdes, F. G., & Allen, S. L. 2006, A&A, 446, 747

Gwyn, S. D. J. 2014, JInst, 9, C04003

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

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

Krist, J. E., Hook, R. N., & Stoehr, F. 2011, Proc. SPIE, 8127, 81270J

Linscott, I., Protopapa, S., Hinson, D. P., et al. 2016, in AAS/DPS Meeting 48 Abstracts, 213.04

Michalik, D., Lindegren, L., & Hobbs, D. 2015, A&A, 574, A115

Nelder, J. A., & Mead, R. 1965, CompJ, 7, 308

Noll, K. S., Grundy, W. M., Chiang, E. I., Margot, J.-L., & Kern, S. D. 2008, in The Solar System Beyond Neptune, ed. M. A. Barucci et al. (Tucson, AZ: Univ. Arizona Press), 345

Ortiz, J. L., Santos-Sanz, P., Sicardy, B., et al. 2017, Natur, 550, 219

Petit, J.-M., Kavelaars, J. J., Gladman, B. J., et al. 2011, AJ, 142, 131

Porter, S. B., Spencer, J. R., Benecchi, S., et al. 2016, ApJL, 828, L15

Silverman, B. W. 1986, Density Estimation for Statistics and Data Analysis (London: Chapman and Hall)

Stern, S. A., Bagenal, F., Enrico, K., et al. 2015, Sci, 350, aad1815

Tchebychev, P. L. 1853, Théorie des mécanismes Connus sous le nom de parallélogrammes (St.-Pétersbourg: Imprimerie de l'Académie impériale des sciences)

The Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, arXiv:1801.02634

Volk, K., & Malhotra, R. 2017, AJ, 154, 62