Following the Pluto flyby of the New Horizons spacecraft the mission provided a unique opportunity to explore the Kuiper Belt in situ. The possibility existed to fly by a Kuiper Belt object KBO as well as to observe additional objects at distances closer than are feasible from Earthorbit facilities. However at the time of launch no KBOs were known about that were accessible by the spacecraft. In this paper we present the results of 10 yr of observations and three uniquely dedicated effortstwo groundbased using the Subaru Suprime Camera the Magellan MegaCam and IMACS Cameras and one with the Hubble Space Telescopeto find such KBOs for study. In this paper we overview the search criteria and strategies employed in our work and detail the analysis efforts to locate and track faint objects in the Galactic plane. We also present a summary of all of the KBOs that were discovered as part of our efforts and how spacecraft targetability was assessed including a detailed description of our astrometric analysis which included development of an extensive secondary calibration network. Overall these efforts resulted in the discovery of 85 KBOs including 11 that became objects for distant observation by New Horizons and 486958 Arrokoth which became the first postPluto flyby destination.
['Trans-Neptunian objects', 'Astrometry', '1705', '80', 'Astrophysics - Earth and Planetary Astrophysics', 'Astrophysics - Instrumentation and Methods for Astrophysics']
The New Horizons Extended Mission Target Arrokoth Search and Discovery
{'The New Horizons Extended Mission Target: Arrokoth Search and Discovery': '<!-- image -->  \nMarc W. Buie  \n,  \n1  \nJohn R. Spencer  \n,  \n1  \nSimon B. Porter  \n,  \n1  \nSusan D. Benecchi  \n,  \n2  \net al.  \n1 Southwest Research Institute, 1301 Walnut St., Suite 400, Boulder, CO 80302 2 Planetary Science Institute, 1700 East Fort Lowell, Suite 106, Tucson, AZ 85719', 'ABSTRACT': 'Following the Pluto fly-by of the New Horizons spacecraft, the mission provided a unique opportunity to explore the Kuiper Belt in-situ. The possibility existed to fly-by a Kuiper Belt object (KBO) as well as to observe additional objects at distances closer than are feasible from earth-orbit facilities. However, at the time of launch no KBOs were known about that were accessible by the spacecraft. In this paper we present the results of 10 years of observations and three uniquely dedicated efforts - two ground-based using the Subaru Suprime Camera, the Magellan MegaCam and IMACS Cameras, and one with the Hubble Space Telescope - to find such KBOs for study. In this paper we overview the search criteria and strategies employed in our work and detail the analysis efforts to locate and track faint objects in the galactic plane. We also present a summary of all of the KBOs that were discovered as part of our efforts and how spacecraft targetability was assessed, including a detailed description of our astrometric analysis which included development of an extensive secondary calibration network. Overall, these efforts resulted in the discovery of 89 KBOs including 11 which became objects for distant observation by New Horizons and (486958) Arrokoth which became the first post-Pluto fly-by destination.', '1. INTRODUCTION': "NASA's New Horizons mission (Stern 2008) flew by Pluto in July 2015 and then deeper into the Kuiper Belt, providing the opportunity of close encounters with other objects. This is the only opportunity for close-up observations of Kuiper Belt objects (KBOs) for the foreseeable future. Prior to the launch of New Horizons , no KBOs were known that were accessible within the propulsion limits of the spacecraft. Finding suitable targets became a critical need for the extended mission. We estimated that objects numerous enough to be targetable with available ∆V were likely to be in the 50 km size range, with apparent R magnitudes near 26 (Spencer et al. 2003). The search area lay in the Milky Way, in Sagittarius, because of the location of Pluto and the trajectory of the spacecraft, requiring KBOs to be found against a very high density of background stars.  \nbuie@boulder.swri.edu  \nFigure 1. Map showing the search areas on the sky during the timespan of this project. The sky is shown with an overlay of the equatorial sky-plane coordinates. Superimposed are ellipses for the area that includes 67% of the encounterable cold-classical KBO population at opposition for the indicated year. The different colors correspond to the three distinct search epochs discussed in this paper (color coded by common effort). No searches were conducted in years with white ellipses. Also shown are curves for the ecliptic and the galactic plane.  \n<!-- image -->  \nA key element affecting our search efforts was the size and orbit distribution of the KBO population. Based on the orbit distribution known even in the early 2000's, it was clear that the spacecraft trajectory would be passing through the cold-classical KBOs (CCKBOs) that are characterized by low inclination and low eccentricity orbits and represent a dynamically cold and presumably undisturbed primordial population. The nature of this component of the KBO population was thus very important to the outcome of the search. The hot classical component provided an additional small chance for an encounterable option. No other dynamical types of object were considered when designing our search program due to the extremely large search area required for even a miniscule probability of detecting an encounterable object from the non-classical populations.  \nThe size distribution also strongly affected the outcome of our search efforts. Our search needed to reach to fainter brightness and thus smaller size than was fully characterized at the time. There were already hints of a break in the slope of the size-frequency distribution when we started the project, but a firm consensus had not yet been reached (eg., Bernstein et al. 2004; Fraser et al. 2008). Ongoing work has provided much better characterization of this important attribute of the CCKBOs (Fraser et al. 2014; Kavelaars et al. 2021). Our earliest estimate of the sky-density of flyby targets used a single-slope power law based on objects in the 20-23 mag brightness range. Now we know that this provided a very optimistic prediction for the outcome of the search.  \nWe will not review all the work and improvements to the characteristics of the general KBO population that was happening during this project. But, we were well aware of such work going on and folded such new results into our search efforts as they came available. Equally important to our work and the orbit estimation needed to support New Horizons mission operations was the astrometric catalog produced by the Gaia mission (Gaia Collaboration et al. 2016a,b; Lindegren et al. 2016). By the end of this project, the astrometric support catalog was a solved problem but along the way this was a major challenge to overcome prior to the Gaia DR2 catalog release (Brown et al. 2018; Lindegren et al. 2018).  \nIn this paper we describe the full extent of our efforts to locate objects for both a fly-by visit by the spacecraft as well as observations of objects at distances that provide higher resolution than the Hubble Space Telescope ( HST ). Fig. 1 shows the search area by year as well as the background field. We first describe our initial wide and shallow search for  \nobjects which began in 2004 and set the stage for determining what actually needed to be done to ensure finding a fly-by target. We then describe our dedicated ground-based searches from 2011-2014 and our HST large program which ran during the summer of 2014 and resulted in the discovery of (486958) Arrokoth which became the New Horizons fly-by target. We describe our various reduction techniques for each of the different observing configurations, groundand space-based, and our efforts to assess and determine the targetability of potential objects both in our search area and then of the discovered objects themselves. Discussion details include the astrometric catalogs used in our work and the development of this process over time. An appendix provides extensive documentation of the mathematics, library routines, and other support tools developed for this project as well.", '2. WIDE AND SHALLOW SEARCH - 2004-2005': 'Our initial search, intended as a pilot study to assess the effort required for a future deeper search, was carried out in 2004 and 2005 with the then newly commissioned SuprimeCam prime focus imager at the Subaru telescope on Mauna Kea. The search area to cover was then about 7 deg 2 . The field-of-view (FOV) of SuprimeCam was 0.3 deg 2 , then the largest FOV on any of the sufficiently large telescopes. The readout time was about 20 seconds, crucial for part of our observing strategy. Additional details about this camera system can be found in Miyazaki et al. (2002). As expected, the biggest problem for this search was the large number of stars in the background. Our intent from the beginning was to use optimal image subtraction as described by Alard & Lupton (1998) to remove the stellar background, leaving the solar system objects to be identified and measured. Our observing strategy for this project was to collect a series of relatively short exposures (either 90 or 120 seconds depending on sky conditions) at multiple distinct epochs, or visits. Using short exposures minimized saturated stars and by stacking many images we built up longer integration times to reach a fainter limit. At this stage in our efforts we did not attempt to shift the images to track target motions. As a result, the length of time for any set of frames to be stacked was kept below the threshold where trailing losses could begin to affect our sensitivity. The limiting magnitude targeted for this search was R = 25 and was largely dictated by the relatively large search area and the amount of telescope time we could obtain. In total, we were allocated 16 nights of time in 2004 and 6 in 2005 to this project and realized 67 hours of integration time on 58 distinct search fields. Table 1 gives an overview of the these runs and Table 2 provides the pointing details and the estimated relative probability of finding a targetable KBO on that field. In a given run a field was imaged on at least 2 nights with two visits on one night and one visit on an adjacent night.  \nFigure 2. Map showing the field pointings for the 2004 search plotted on the plane of the sky in equatorial coordinates. The solid blue line is the galactic equator and the center of the displayed segment is just 10 · from the galactic center. The solid green line at the bottom is the ecliptic. The green boxes are the size of the SuprimeCam FOV and show the positions from Table 2. The solid ellipse shows the 1σ search area in 2004 and the dashed ellipse shows the search area in 2005.  \n<!-- image -->  \nThe search pattern for this phase of the project is shown in Fig. 2. Each tile in this pattern represents a single FOV on the sky for a single pointing of the telescope (Table 2). The gridded area covers 87% of the distribution of targetable objects on the sky. The highest probability tile contained 4.2% of the distribution. The goal for each surveyed field  \nBuie et al.  \nTable 1. Summary of 2004 Search Runs  \nNote - Total time is the aggregate on-target exposure time in hours.  \nwas to get five visits over a small number of days to serve as the discovery dataset. This choice was dictated by the software requirements of the Japanese search team members who already had software for the search based on this required number of visits. Additional observations in 2004 and 2005, consisting of pairs or triplets of visits were used for followup observations. The coverage of fields in the search area was such that covering the search area also served for followup without needing to specifically target the newly discovered objects. This part of our strategy was driven by the realized expectation that data analysis would take a very long time to complete compared to the observational phase of the project.', '2.1. Image calibration and characterization': 'The raw images were processed using standard techniques prior to analysis. These steps include overscan corrections, bias-frame corrections, and flat-fielding corrections. Our flat fields were generated from the imaging data directly with a robust stack of all data. Building flats and other image stacking operations required a good characterization of the sky signal in each image. We developed a new routine ( skyfit.pro 1 ), that includes robust estimation and simultaneous fitting of a low-order 2-D polynomial to the sky background that uses incomplete information on the distribution of sky signal due to the crowded fields to reconstruct the full noise distribution. See Appendix B for more details about this process. The spatially variable sky signal was subtracted from the images, and the mean value near the center of the array was added back to roughly preserve the original detected signal level.  \nThe required data processing was facilitated by collecting key derived quantities for each image and saving all of this in a relational database to prevent having to recompute these values in each step of the processing. Information collected in this way was the image quality, photometric depth, sky signal properties, detector gain. Some derived information was too voluminous to be practical for storage in the database, such as lists of detected sources and the full astrometric solution for each image at every step in that process. This derived information was very valuable for examining data from a given night as a whole to help detect those images where the processing failed for some reason.', '2.2. Astrometric support catalogs': "Accurate astrometric solutions for the images is essential to the stacking and subtraction process. SuprimeCam exhibits a large amount of spatial distortion in the focal plane between the center and corners of the mosaic. The positions near the corner are offset compared to what a linear focal plane would provide. It was not practical to control the pointing for the search images to the sub-pixel accuracy needed for stacking so the distortion had to be directly addressed.  \nWe were well aware of the need to avoid saturation, hence the plan to take multiple short exposures and stack. Under good sky conditions, the limiting magnitude was 23.9 mag. More importantly, stars brighter than 18.2 were saturated and of no use for astrometric calibration. This brightness limit meant that the best data had an insufficient number of astrometric reference stars to permit determining an astrometric solution of the distorted field.  \nTable 2. Summary of 2004 Search Fields  \nNote -(1) The coordinates are all equinox J2000. (2) 'Prob' is the estimated relative probability of finding a targetable KBO.  \nWe used the USNO-B1.0 catalog (Monet et al. (2003)) for this phase of the project. This catalog includes proper motions, but we used the form of the catalog extracted for use with the Deep Ecliptic Survey (Elliot et al. 2005) where the format of the earlier USNO-A1.0 catalog files was used with all stars saved at an epoch of 2005. Our first attempt was to use this catalog on a chip-by-chip reduction, but the combination of strong distortion and minimal overlap in brightness between the catalog and images made this simple approach insufficient. Obtaining good astrometric  \nsolutions required two adjustments in procedure. First, we had to create a secondary astrometric catalog that went to a fainter limiting magnitude. Second, we developed tools to treat the entire mosaic as a single astrometric device.  \nTo create the secondary catalog we identified the data taken with the worst transparency from each field. Every field had at least a few images taken with several magnitudes of extinction. These low-throughput images could then be processed using the USNO-B catalog since most of the catalog stars were now not saturated. More importantly, even though these images were photometrically shallower, they still pushed significantly deeper than the USNO-B catalog and we could extract the fainter stars for a deeper secondary catalog. Note that the secondary catalog was constructed from the average of all data. No attempt was made to correct or otherwise adjust these positions for proper motion.  \nObtaining an astrometric solution was also complicated. The large number of frames pushed us to automatic pipeline processing rather than a supervised manual operation. The automatic tools available then were effective on only about half of the data. To overcome this problem, we used the focal plane distortion map provided by NAOJ (Miyazaki et al. 2002) for an initial mapping of pixel coordinates onto the sky plane. This map was provided in the form of software rather than as formulas and coefficients.  \nFor each image, we used the NAOJ tool to convert the pixel positions of all measured sources in the image to RA and Dec given the knowledge of the telescope pointing. The result of this step served to linearize the positions. From here, the RA and Dec estimates were converted to tangent plane coordinates using the header position. The overlapping secondary star catalog positions were also converted to tangent plane coordinates using the same tangent point. These two lists were close to each other in position, but always had offsets to be determined as well as occasionally needing to be rotated up to a couple of degrees.  \nOnce the lists were correlated, a new astrometric solution was relatively easy to derive. To support a full-mosaic coordinate system, we developed our own method for describing the focal plane distortion. This mathematical description and its mapping onto a set of FITS keywords is described more fully in Appendix A. Deriving a full solution required an iterative approach for our fourth order fit. Given a solution, we calculated the new location of the center of the mosaic, thus changing the tangent point. This required a new solution, leading to a new tangent point. This process converged within 3 iterations to provide a self-consistent tangent point and astrometric fit.  \nAnother component of the astrometric fit is the relative offset and orientation of all the detectors. One chip (#2) was selected as the reference device and its coordinate system is extended to cover the full focal plane. In principle, the choice of reference chip is arbitrary, but we found it convenient to use a device in the center of the mosaic where the distortion is lower. The original distortion map for SuprimeCam was used to initially determine these chip offsets and orientations. For all images, the fit is held constant and the offset and orientation is adjusted to find the best fitting solution.  \nIn principle, this method assumes that the astrometric solution can be different on each image. The constant and linear terms are always seen to vary with time. The quadratic and higher-order terms were never seen to exhibit statistically significant variations within a run. We found the results to be more internally consistent and stable when the mean values of the high-order terms were imposed on the fit, leaving only the constant and linear terms to be fitted for each image.  \nAdditionally, the relative locations of the chips was not expected to vary. This expectation seems like a reasonable assumption, but we chose a more conservative assumption that the devices don't perceptibly move during an observing run. Finding the offsets is not entirely straightforward. In general, there is a coupling between the astrometric solution for the entire mosaic and the values of the offsets and rotations used for each chip. We did not have the computing resources and tools at that time to attempt a simultaneous determination of the solution and the chip positions. Instead, we iteratively determined one while fixing the other and repeating until the answer converged. This proved to be a very effective method if somewhat tedious. Doing this required a starting guess for the offsets and rotations which we had from the NAOJ documentation. Then, we fixed these chip locations and solved for the chip-by-chip astrometric solution. These solutions where then averaged to get mean high-order coefficients and the linear and constant terms were refit. From this point, the solutions were fixed and the chip offset and rotation relative to the reference was scanned to find the values most consistent with the fixed solution. This is done for all images in a run. Given these values, a mean offset and rotation is computed using a robust mean estimator ( robomean.pro ). These mean values were then adopted, and the solutions were iterated again. This process was repeated two or three times for each run.  \nThe measure of success in all of this intense astrometric calibration was in looking at multiple views of the same patch of sky. A high-speed animation loop looking at the same region registered by their astrometric solutions should only reveal the usual variations in transparency and image quality induced by the atmosphere, but there should be no  \napparent motion in the positions of the sources. In the early stages when chip-by-chip independent solutions were used, there was a lot of residual motion - sometimes many pixels in the corner chips. After this more involved calibration process, the images were seen to be internally consistent. Given that this process is referenced to an astrometric catalog, the final absolute astrometric positions should be as good as the underlying catalog.", '2.3. Photometric calibration': 'The usual process of flat-fielding an image preserves surface brightness. This is useful for removing pixel-to-pixel variations, over small and large spatial scales. In the presence of optical distortions, flat fielding can introduce systematic shifts to point-source photometry. The reason for this is that the area of a pixel, on the plane of the sky, varies with such distortion. To correct for the effect of this distortion on fixed-aperture photometry, we compute the area of each pixel in arcsec 2 . The map is normalized by the largest area. Then, the flux of any source is divided by the value from the normalized map at the location it was measured before converting to an instrumental magnitude. For SuprimeCam, this correction is as large as 30% towards the more heavily distorted corners of the mosaic. The computation of this map is time-consuming but straightfoward given a good astrometric description of an image. We did not apply this correction to the images themselves, but used it to correct the extracted photometry.', '2.4. Image registration': 'Given an accurate astrometric solution, image registration is straightforward. In the early phase of the project we tried to determine a polynomial expansion approximation for the mapping from one image to another image. The goal in this case was to pick the best image of a set and then resample the other images to match. Operationally, this method should work but the high-order differences between two distorted images required too many terms and that fit was numerically unstable for a significant number of images. A better methodology turned out to be to resample all images to a linear grid defined by the mean position and orientation of all images of the same field within a run. Image scale for the output images was set to the scale of the image at the center of the mosaic. The astrometric fit already determined was sufficient to map the observed images onto the idealized grid. The output grid was sized to bound all of the input images, meaning that there could be pixels that map onto only a subset of the input images.  \nThe images were all re-sampled with a flux conserving interpolation that is functionally equivalent to sinc interpolation. This process used the Buie library routine dewarp.pro that itself uses the built-in IDL image interpolation tool, congrid .', '2.5. Image subtraction': 'The power of image subtraction for picking out changes in images cannot be denied. What sets this project apart is its time-dependent nature. We routinely needed to work on data with a much larger range of image quality than is typically tolerated by other investigations of a largely static sky. The quality of subtraction was seen to degrade as the difference in seeing between two images increased.  \nThe first step in the subtraction process is to take all of the images for a given field that are consecutive in time and average them together. This averaging process rejects outliers such as cosmic ray strikes and fast moving objects in the field but preserves stars and KBOs. Prior to the stacking process, each image had its spatially variable sky level subtracted and then flux scaled to match the flux of the best image in the set. The constant term of the sky fit from the best image was then added to the average to return the image to pixel values closely approximating the original to facilitate tracking the noise properties of the images. All images in a set were used to build the stack.  \nAny pair of stacked images can be used for image subtraction at this point, but not all pairs were equally effective. For each image, we computed the difference against the best quality stack. The final difference pairs were computed and saved for the next step of searching.', '2.6. Searching for objects': "After the visit difference images were completed, we visually examined the data to look for slow-moving objects. This process was essentially the same as used by the Deep Ecliptic Survey in using color composites for comparing two epochs of data at a time. Traditional image blinking could be used if desired. When an object was noted, its position was measured in the difference image. Since the astrometric grid is known, it then becomes a simple matter to convert from pixel coordinates to sky coordinates. This tool used a low-order polynomial short-arc predictor to project the position of a measured object onto each image not yet measured. This system allowed finding an object across the chip gaps.  \nMuch later in the project, we implemented a parallel system for finding objects as a citizen-science project. A broad discussion of what citizen science entails can be found at https://www.citizenscience.gov/. In our case, we engaged our participants by providing image data to search for transient sources that included variable stars and moving objects with distinctions between main-belt asteroids and KBOs. This system was called 'Ice Hunters' and hosted as part of the 'Zoo' collection of citizen science arenas https://www.zooniverse.org/. After analyzing data from a number of observing runs we found no difference in the detection list between our scanning and Ice Hunters. In the end we found that technology constraints at the time served to make the citizen-science approach much slower and more effort to manage than just using a local tool with a couple of people working to scan the data. More details about this experiment can be found in Appendix C.", '2.7. Wide and shallow search results': "Our search effort was able to find 24 objects. Table 3 summarizes our results. The faintest object seen was at R=24.9, exactly as expected from our initial estimates prior to the search. Of those objects detected, 14 (58%) were likely cold-classical and the lowest ∆V for New Horizons was 1.1 km/s, or about 10 times bigger than the available amount. The closest candidate is one of two objects that are currently designated (though clearly there are others with good enough orbits for designations). NI09 (or 2004 LW 31 ) now has a much better orbit than was available from just our search data. It now has a 15-year arc with 39 unique observations (including ours). The ∆V estimate shown in Table 3 for this object is good to 10% showing that even short-arc data are sufficient to strongly constrain the list of potential targets. The orbit estimates for these objects as well as those in following efforts evolved a great deal during this project. The initial tools used were based on the approach of Bernstein & Khushalani (2000). Our tools evolved to include a more direct probabilistic estimate meant to improve on the uncertainties and be more ameable to addressing encounterability. By the end of this project, we moved to a full Markov-Chain Monte-Carlo (MCMC) approach that is described elsewhere (Porter et al. 2018, 2022). The orbital elements tabulated in this work should not be taken as the final word on these objects. They were vetted only enough to assess the encounter probability.  \nEven though this early search work did not result in the detection of a mission target, important lessons were learned that were carried forward to the later search efforts. Here we summarize the more important of these lessons. Also, one of these objects (2004 LW 31 ) made it to the list of bodies to be observed with long-range encounters by New Horizons before and after its encounter with Arrokoth (Porter et al. 2022).  \nThe search area was variegated enough to include some areas of 'dark lanes' where foreground dust in our galaxy obscures the stars that would normally have been seen. We did not specifically target these areas, but recognized that they would enhance our probability of finding objects due to a reduced amount of confusion. Unfortunately, nothing useful was found in those limited areas.", '2.8. Lessons learned': 'The results of this effort showed us that we clearly needed a much deeper search to have a reasonable chance to find an extended mission target. This effort also clearly showed the critical need for a secondary reference star catalog that could be directly tied to the search data. It likewise demonstrated the influence of high quality seeing ( < 0 . 7 arcsec) on our object discovery rate; quality seeing with lower transparency was better than poor seeing with clear skies. All of these lessons provided essential information that guided the implementation of the next search effort.  \nFinally, nearly all the software tools and the lessons learned therein were then available as a starting point for follow-on searches that were characterized by ever shortening timescales for data processing. This initial search effort did not result in finding a targetable object, but was an absolutely essential step along the way.', '3. DEEP GROUND-BASED SEARCH - 2011-2014': "A deeper search was undertaken starting in 2011. By this time, the search area for 80% of the candidate pool had dropped down to ∼ 2 deg 2 and we needed only 7 discrete SuprimeCam fields to cover the search region. The reduction in the area made it feasible to consider much longer integrations and thus deeper image stacks to reach a fainter limiting magnitude while still covering most, if not all, of the search area. Once again we used the Subaru SuprimeCam imager, but there was a second option available with Megacam (McLeod et al. (2015))and IMACS (Dressler et al. (2006)) on the Magellan telescopes. A list of observing runs for this search effort which includes a run identification, the instrument used, the UT start and end time for each run, the number of distinct fields imaged and the number of individual frames collected can be found in Table 4.  \nTable 3. Summary of 2004 Search Discoveries  \nNote -The astrometry for these objects from this search is provided in the supplemental data file for this table. All values given here are from the discovery data only. Blank ∆V values indicate that close approach conflicted with the Pluto encounter (on or before). ∆V values listed as ' > ' indicate a non-physical value.", '3.1. Astrometric support catalog': "Driven by the lessons learned in the initial 2004-2005 pilot search, we developed a special astrometric support catalog for the calibration of the new search data. Targeted observations with the Canada-France-Hawaii Telescope (CFHT) were obtained in queue mode. The goal each year was to cover the sky for the following year so that the catalog was ready for each new apparition of search data. All observations were taken with MegaCam (Boulade et al. 2003) and an r ' filter with 30 second exposures. Good image quality was essential for these crowded fields and the queue observing system ensured we got the necessary quality data. Most data were taken under photometric conditions even though this wasn't absolutely required. The data were processed at CADC using their standard tools for the MegaCam mosaic (Gwyn 2014). With its 1 deg 2 FOV camera it was able to easily cover the search region. Toward the end of the search phase this catalog was extended to cover the entire field through the expected target position at the time of the New Horizons encounter. The catalog reaches a useful depth of r=23-24 which provided ample overlap with the brightness range of the ground-based data and the available astrometric star catalogs. Initially, the UCAC4 catalog (Zacharias et al. (2012)) was used as the underlying reference upon which the support catalog was based. As the project progressed we moved to hybrid catalogs, URAT1 from Zacharias et al. (2015), and the Gaia DR1 catalog (Gaia Collaboration et al. 2016a,b; Lindegren et al. 2016), and then ultimately to the Gaia DR2 catalog (Brown et al. 2018; Lindegren et al. 2018). Mean proper motion corrections were applied so that the support catalog had a mean  \nepoch contemporaneous with the search. This catalog ensured that there were at least a couple of hundred reference stars for each detector in the search cameras.  \nFigure 3. Sample SuprimeCam image data from 2011. Two epochs are shown of the discovery images for 2011 HL 103 . The red arrow (or image) shows the location at the first epoch and the cyan arrow (or image) shows the location at the second epoch. The stacked images are in the left column. Template subtracted images are in the middle column. The right column shows the color stack of the two epochs. See text for more details.  \n<!-- image -->", '3.2. Data Processing': "The data analysis for this phase of the search was very similar to that used for the prior search, but there were some improvements and other changes needed. First, we had multiple independent pipelines for searching the data. Algorithmically, these pipelines performed similar operations and are not separately described. The Subaru data were again processed as a single multi-chip astrometric device. The Magellan data were processed with each image in the detector mosaics as independent systems.  \nThis search used much longer dwell times per visit to extend the search to fainter objects, but still took short exposures to be stacked. An important change was in the creation of template images for each pointing. This time, we used all available images to build a master template image optimized for each pointing for each telescope. Given the range of image quality in the data, we could not use all images for the template but needed to ensure a sufficient range of time so that any KBOs were filtered out of the template images.  \nTo make a template image, all of the candidate images were sorted by photometric depth. We tried using a list sorted by image quality, but that does not de-weight images that might be taken through some clouds. The photometric depth is dependent on both transparency and image quality and worked best as the sorting metric. A reference image was defined based on the mean pointing of all images, the mean orientation on the sky, and a linear image scale that matched the nominal scale of the images in the center of the mosaics. All images were resampled to this reference.  \nAs the list of images is traversed from good to bad, each new image was assessed for its ability to improve the stack. We used a simple noise model based on the sky signal, sky noise, and photometric depth of the image to estimate  \nKEM target search  \nTable 4. Summary of 2011-2014 Search Runs  \nNote -Nf is the number of distinct fields imaged and Ni is the number of individual images collected. The UT times are the start times of the first and last images collected in that run.  \nTable 5. Logs Notes from 2011-2014 Observing Runs (Sample 1 )  \nNote -(1) This table is available in its entirety in machine-readable and Virtual Observatory (VO) forms in the online journal. A portion is shown here for guidance regarding its form and content.  \nthe change in the SNR of the proposed new stack from adding each new image. If the new image would improve the SNR, it was then added to the stack and the SNR was updated. Any image not improving the SNR was dropped from inclusion in the template. In all cases, every search field had enough data and enough distinct epochs to ensure a good template where moving objects and cosmic ray strikes were filtered out. Once the templates were complete, all individual short-exposure images were differenced against the appropriate template image.  \nThe consequence of moving to longer dwell times was that trailing losses were no longer negligible, even for KBOs. A simple stack as done for the previous search would leave a smeared image of a KBO that would adversely affect our limiting magnitude. To generate the shift needed we used a population of synthetic KBOs that matched the orbital element distribution of cold-classical KBOs. The mean velocity of this population was computed for each visit, field, and device within the detector mosaic. This single shift value was applied to the subtracted images requiring another resampling of the images. Once these images were shifted, they could then be stacked. At this stage, cosmic ray strikes and some of the residuals of the stellar background were removed. The resulting image products could be searched using our standard color-composite and blinking tools.  \nA very small sample image that is typical of our dataset is shown in Fig. 3. This sample image shows two of the twenty-two epochs of data on object 2011 HL 103 . Each of these images is a stack of seven 90-second images for a total of 630 seconds of integration. The first epoch had a FWHM of 3.4 pixels and the second epoch was 3.5 pixels. All displayed images are stretched identically from -3 σ to +8 σ about the sky mean where σ is the measured noise in the sky at the first epoch. In the figure, the left column shows the stacked images from the same region in the two images with the first epoch on the top and the second epoch on the bottom. The arrows point to the position of the object and the red arrow is for the first epoch and the cyan arrow is for the second epoch. The middle column shows the same data but after the template image is subtracted. The object is clearly visible at both epochs in the each type of presentation, but is clearly easier to see once the background stars are subtracted.  \nThe images in the right column show what we examined to discover the objects. On the top is a color composite of both epochs where the red channel contains the first epoch and the green and blue channels contain the second epoch. The static sky shows as shades of grey to white. The transient component shows as either red or cyan depending on its location or properties in the two epochs. Once again the object can be seen in both epochs, but is much easier to spot in the lower right where the background template has been subtracted. It was this background subtracted color image composite that was used for all object detections. Note that additional transient features are apparent in this image product. Most notable are three apparent variable stars, two of which were brighter in the first epoch and thus show as red, and one that was brighter in the second epoch and shows as cyan. These transients are discarded during the search when there is no pairing at a sensible KBOs sky-plane velocity.  \nIn addition to the manual search process described above, the 2011 data were also searched by automatic computer searching methods as well as by the 'Ice Hunters.' As with the earlier results, Ice Hunters did no better than our internal team search efforts and was not pursued for later epochs of search data. The computer search tools did provide some complementary coverage to the visual searching that increased the aggregate completeness of the search effort. From 2012 onward, the search continued with just in-house processing and search efforts by manual and computer searches.  \n3.3. Results  \nTable 6 . Summary of 2011-2014 Search Discoveries  \nKEM target search  \nTable 6 (continued)  \nNote -(1) R magnitude is the apparent brightness at the time of discovery.  \nThe search depth of the survey was quite variable due to weather and site conditions. The best dataset from the entire effort was at Magellan where we reached a 50% efficiency detection limit of R=26.0 in multi-hour stares with exceptional seeing of ∼ 0 . 45 arcsec. While seeing on individual images ranging from 0.4 greater than 2.0 arcsec, the  \naverage was 0.84, 1.06, 0.75 and 0.71 arcsec per year respectively, from 2011-2014 (Table 5). Likewise, each run was given a unique set of field names, but with significant overlapping sky area between the runs in different months. This allowed for orbit extensions on any new objects to determine their observability/accessability by the spacecraft. In effect, all of the project search area was observed each year with multi-month overlap of the fields as conditions allowed. 2011 was our most productive year accumulating 12.5 degrees 2 of sky area (covering the search area multiple times at different epochs). The following years accumulated 4.4, 2.17 and 0.48 degrees 2 , respectively (Table 4). All of our observations were submitted to the MPC using the standard 'OBS80' formatting at that time. Although it is now possible to include detailed error analysis in ADES-format submissions, it was not standard at the time, and we did not re-evaluate this information given the level of effort that would have been required post-facto. Astrometry for objects that we later observed with New Horizons were, however, re-analyzed with those uncertainties and that process is reported in detail in Porter et al. (2018, 2022).  \nTable 6 provides a summary of the objects discovered during this three year search effort, sorted by increasing ∆V. The information provided also includes the object orbital elements ( a,e,i ), the estimated H V -magnitude, measured R-magnitude, accumulated object arc, discovery date, number of observations, MPC designation if one exists, the dynamical classification for the object, the encounter date and spacecraft encounter distance as well as the delta velocity required for such an encounter. The dynamical classification shown is meant to be a rough guide where 'C' is for some type of classical KBO, 'S' is a scattered object, and those with pairs of numbers denote the type of mean-motion resonance with Neptune. A close-approach date prior to the planned date of the post-Pluto-encounter trajectory correction maneuver (TCM) has no ∆V listed. One object had a close approach just days after the TCM planning date and would have required an extremely large ∆V. In this case only ' > ' is shown to denote a large value that is of no interest. The final object in the table is believed to be a real detection of a KBO, but was only seen twice and had insufficient information to constrain its orbit well enough to compute useful close approach circumstances.  \nThe result of the efforts through the end of 2013 made a significant improvement by finding objects that New Horizons would pass close to but the closest, 2011JY31, was not close enough, requiring a factor of two more ∆V than available. Without any course adjustment we expected to get to within 0.15 au of 2011JY31 on 7 October 2018. If we used all of the available fuel, this would still leave us with a flyby distance of ∼ 10 million km and no fuel to support the encounter. As with the earlier search, many of the objects we found became useful distant encounter targets for New Horizons .  \nThe results to this point were such that most of the search area was covered to a brightness limit that could be reasonably reached by the best ground-based facilities available. There was a still a slim chance that we would find something suitable in 2014, but it became clear that another tool was required to ensure finding a target. The tool of choice to extend the search was the Hubble Space Telescope (HST). The use of HST was by no means assured so a parallel full search effort was also carried out from ground-based facilities in 2014. As we began to consider the use of HST, the Hyper Suprime-camera (HSC) was commissioned for use on Subaru (Miyazaki et al. (2018)). This new camera had a field of view of 1.5 deg 2 , considerably larger than the search area for that year. The time allocated for our search on HSC was granted as shared-risk and the weather conditions and instrument problems prevented any new and useful search data. In the end, the ground-based data from 2014 and the following years were used only for recoveries and orbit refinement work for the distant encounter targets and will not be discussed further in this work. A subsequent new search for additional fly-by targets begun on Subaru in 2020 and onward will be discussed in other papers such as Fraser et al. (2024).", '4. HST SEARCH - 2014': "The immutable properties of the New Horizons flyby of Pluto gave rise to strong constraints on the search efforts. Re-targeting the spacecraft to the extended mission destination would use the least amount of fuel if it was done as soon as possible after the Pluto encounter. Due to the slow data transmission time from New Horizons to the ground, the project wanted a significant amount of time to downlink critical data before re-targeting. The time for the trajectory correction maneuver (TCM) was set for October 2015. The analysis and commanding development for the TCM usually takes a couple of months so when that work began, a good orbit had to be known for the new target. It was not enough to have just discovered the object prior to the TCM, we had to find it soon enough to get sufficient data to get an accurate orbit for the object as well. Our analysis showed that we needed astrometry spanning at least a year from which an accurate enough orbit estimate could be obtained to support the TCM. These constraints meant we absolutely had to find an object in 2014.  \nThe obvious advantages of HST are the stable PSF and the very fine image scale. The main difficulty with using HST is that the best survey camera for this project, the Wide Field Camera 3 (WFC3)/UVIS (with a image scale of 0.0395 arcsec/pixel) has a rather small FOV at just 0.0021 deg 2 (Dressel (2019)). Fortunately, 2014 was the year where the search area reached its minimum size of 0.3 deg 2 . Covering 80% of this area thus would require 83 separate pointings of the telescope.  \nIn the end, we secured two allocations of time on HST - 40 orbits for a pilot study which if successful would lead to another 154 orbits for the full survey. The pilot study was supported by Director's Discretionary (DD) time while the full survey came from the Cycle 23 GO allocation (GO-13633, PI Spencer). The agreed upon success criteria for the pilot study was to find two cold-classical KBOs as a rudimentary check on the predicted size distribution down to this small size. The pilot was successful in finding two objects and we were given clearance to conduct the full survey (83 search fields - 7.29 sq arcmin each, ∼ 605 sq arcmin total). 28 orbits from the full survey allocation were available for additional followup observations on newly discovered objects.", '4.1. Survey design': "The goal of the search was to find an object that New Horizons could reach for a close flyby investigation. This goal allowed a much narrower scope for the search than a generic search for KBOs. In this case, we could fine tune the data collection and processing to maximize our ability to find the object we sought. It was entirely permissible to be blind to objects that New Horizons could not reach.  \nThe area to search was defined by two primary constraints, (1) the region of space that New Horizons could reach after October 2015 and (2) the properties of the Kuiper Belt. The predicted trajectory of New Horizons defined a time-dependent path through space and the amount of fuel on board set the size of a cone centered on that path. Any KBO in that cone of 130 m/sec ∆V at the right time is encounterable. The requirement of being in the cone at the right time provided a constraint on the orbit an object would need to have. Our knowledge of the spatial distribution of the Kuiper Belt and the size distribution of KBOs was used with that positional constraint to construct a virtual population of KBOs that were a subset of all objects having the right orbital properties to be at the right place at the right time.  \nThere were only two dynamical groups that had a useful probability of being in the accessible cone: cold-classical and hot-classical objects. Of these, the cold-classicals were strongly favored over hot-classicals with a ratio of 4:1. We adopted a KBO population model which used the orbital distribution derived from the Canada-France Ecliptic Plane Survey (CFEPS) model of the Kuiper Belt Petit et al. (2011). The size distribution was constrained by the absolute and apparent magnitude distribution of KBOs from recent models (Fuentes et al. (2010), Fraser et al. (2014)). At the time of the survey proposal we deemed this 'FEPS/Fraser' model to be somewhat conservative; the true number to be found might actually be higher. We then used this orbit population to derive a set of semi-major axis, eccentricity, and inclination values that were consistent with this model of the Kuiper Belt. From there, we selected random values of the angular elements. If the randomly drawn elements placed the virtual object in the accessible cone, the elements were saved. At the end of this process we had a selection of orbital elements representing virtual objects that were all encounterable.  \nGiven the set of encounterable virtual objects, we could then project their positions back to some prior time and thus know the region of sky that needed to be searched to find a real object. Figure 4 shows a plot of the sky plane centered on the cloud of virtual objects in our sample. Each dot represents one virtual object and the cloud is roughly described by a two-dimensional Gaussian distribution. Each square represents one field of the WFC3 camera on HST and the fields are laid out to cover 80% of the distribution. We chose this coverage as a compromise between completeness and a practical request of time for the search. Covering a larger region would have taken much more time with quickly diminishing returns for that time investment. The properties of this cloud of virtual objects was time-dependent, but we ignored this in designing the search, provided the data could be obtained quickly enough. However, the orientation of the imaging camera is governed by the roll angle of the spacecraft as a function of time. The figure depicts the image orientation at the time of the pilot survey, shown in green. The remainder of the search fields, shown in blue, have the orientation expected for that part of the full search. In reality, the true search grid could not avoid small gaps between adjacent fields, but the missed search area was small enough to ignore.  \nA special virtual orbit was chosen from the entire cloud. This orbit was the closest to the center of the cloud and was used to set the pointing and tracking for HST. This strategy works well because in this case, the differential rate of motion of the cloud of virtual objects was 0.5 pixels (20 mas) during a single exposure for fields on the outer edge  \nof the distribution thus avoiding any trailing losses for objects of interest. Other KBOs, such as resonant or scattered disk objects, would have different rates of motion and were not expected to be detected at all. All search observations were based on this same guiding orbit throughout the survey.  \nFigure 4. A plot of the 83 HST search fields. Each box marks the approximate area of a combined WFC3/UVIS field using both detectors. 5 fields ended up not being useful due to guide star failures (red) or moonlight (orange). The green boxes represent the 40 orbit HST pilot fields and the blue show the additional full program fields excluding the repeated fields to follow up on the detected objects. The 5 discoveries from the search are marked with a red X and labeled with the eventual object designation. The black dots are a simulated KBO population from the CFEPS model with hot and cold Classical KBOs included.  \n<!-- image -->  \nThe search observations were designed as a pair of orbits looking at a given search field from the set shown in Fig. 4. There was a small amount of motion during one visit and sufficient motion between two visits to get an initial estimate of the target distance and orbital inclination. This two orbit time span was sufficient to produce an initial orbit constraint that reduced the list of objects requiring followup and further orbit refinement, but was not in itself enough to know if an object was encounterable. However, by the very nature of the survey design, all detected objects were highly likely to be interesting.  \nThe optimal time spacing was for back-to-back orbits. In some cases, this was impossible so 4.2% of the fields had a one-orbit gap between the two epochs. The increased time base was marginally useful for constraining the orbit of a detected object, but came at the price of a much larger search effort. Even longer time gaps were categorically excluded from consideration in scheduling observations. No other constraints were placed on scheduling the search observations other than to pick a set of fields near the center of the distribution for the pilot project. In retrospect,  \nTable 7. Summary of New Horizons Search Fields  \nNote -(This table is available in its entirety in machine-readable form in the online journal. A portion is shown here for guidance regarding its form and content.)  \nwe should have increased the lunar avoidance distance constraint rather than using the default settings. Two of our search fields (marked in orange) were taken close enough to a full moon that the extra scattered light made the images useless for the survey. Table 7 gives a list of each HST field and the results of our searches and Fig. 5 shows a sample image pair.  \nFigure 5. Sample image from field KBO1-R7R6 in each of 2 visits separated by 1.05 hours. This field was part of the pilot search and 2014 MU69 is in the top chip on the right side in the area marked by the white box. The visit numbers of HST program 13633 are shown on the lower right hand corner of each image.  \n<!-- image -->  \nWe used the WFC3/UVIS instrument with the F350LP filter for all observations. The data volume for the full FOV readout led to a particular exposure time choice. Five 370-second exposures fully fill a single visit without losing time to data management overhead. Some visits needed a small tweak (no more than a few seconds) to get the sequence of  \nobservations to fit in a single visibility window. The pair of visits was further constrained to be taken with the same guide star pair and same roll angle. Due to the guiding orbit used, all background stars were trailed in the image. This factor in the images required the development of new techniques to obtain accurate astrometry of the discovered KBOs.", '4.2. Data analysis': 'Each set of five images from a single visit was processed to generate a stellar template image. Referenced and stacked on the star positions, the stack is a median of the images. This template had only stars with cosmic-ray strikes and KBOs naturally filtered out. Each image was then processed to subtract the template image. This step suppressed the stars, including the wings of the stellar images. In the cores, there is too much residual noise to permit searching those areas for KBOs. In the wings, the KBOs can be completely hidden until the faint outer part of the star image is subtracted. For most visits, the trailing of the stellar images is not exactly the same in all the orbits. As a result, the stacked template does not completely remove the core of the star trails and leaves a characteristic residual image.  \nThe next step was stacking the set of five images. This process was done twice, once with and once without subtracting the stars. Two things happen at the same time: the images are re-sampled to eliminate the distortion and register on the same linear 40 mas/pixel grid while shifting the image based on a given velocity vector. Since the images are not Nyquist sampled nor cleaned of cosmic ray strikes and other blemishes, we chose to use bilinear interpolation for the re-sampling. This choice was cosmetically much better than using a proper sinc interpolation.', '4.3. Rate Corrections': "In order to sensitize the search to KBOs on targetable orbits which differed from the nominal guiding orbit, we conducted a 'shift-and-stack' digital tracking approach. Frames were shifted in software to counter the differential sky-plane motion between the nominal guide orbit and other potential targetable KBO orbits. As mentioned previously, the source population of targetable KBO orbits was drawn from the 'CFEPS/Fraser' model and they were downselected to a small set of guide orbits using the methods of Parker & Kavelaars (2010). Unlike previous KBO searches using the digital tracking approach, frame-by-frame shifts were not parameterized as a grid of linear rates of motion in the sky plane; instead, the frame-by-frame motion of each guide orbit was computed and a frame-wise sky-plane offset was computed from the exact motion of a KBO on that orbit as seen from the HST reference frame. The down-selected orbit grid was chosen with a density such that no KBO in the large source sample of synthetic targetable KBOs would have a single frame-wise correction that differed by more than a single HST WFC3 pixel ( ∼ 0.039 arcsec).  \nThe sample of virtual objects used to design the search provided a means to calculate the bounds for the search in velocity space. The size for the velocity bins was set such that the worst-case smear within a bin was no more than one pixel across the five images in a single visit. There are two important consequences of this design. First, an object is unlikely to be seen in more than one velocity bin, especially at the detection limit of the survey. Second, the number of discrete stacks to search varied with the time-span of the set. For observation pairs of the same field that occurred on back-to-back HST orbits (1.05 hours), typically of order 20 guide orbits were required to achieve the required sampling density leading to 20 image stacks to search. For observation pairs that spanned larger periods of time, a rapidly-increasing density ( ∝ t 2 ) of guide orbits was required and we ended up with as many as 80 stacks to search. Due to the increased workload and need for rapid turn-around, we adjusted our observation constraints as the program progressed to strongly prefer the shorter time spacing.  \nFigure 6 illustrates one resulting grid of linear sky-plane offsets - for the entire population of targetable orbits as well as the chosen guide orbits - from the first image in a single orbit to the last, computed for visit icii11 (where (486958) Arrokoth was discovered).  \nFigure 6. Sky-plane offsets from the nominal guide orbit over the course of a single HST orbit for visit icii11. Black points are the motion of all targetable KBOs from the synthetic orbit sample. Red points are the sky plane motions of the chosen guide orbit grid, and blue circles are the buffers of sky motion each covers to within a single HST WFC3 pixel. The point at (0,0) is the nominal guide orbit.  \n<!-- image -->  \nWe had five nearly independent teams pursuing their own search tools all built more or less around the same concepts of the shift-stacking. This discussion largely concentrates on the pipeline written by Buie since that was responsible for all five of the discoveries within the survey project. We will touch on aspects of the other systems as needed.  \nIn one of the final stages of development for the software, we implanted fake objects into our data. This was designed primarily to measure the sensitivity limits of our data. In the end, it was far more important as a training dataset to test the software. During the pilot phase, significant portions of the data were scanned multiple times with no results. The test data convinced us that we had the expected sensitivity and that the software was only finding objects when the shift was zero. This was quickly determined to indicate a sign error in the velocity stacking and Buie's detections in the pilot survey quickly followed after correcting the error.", '4.3.1. Searching for Candidates': 'The search methodology employed a mix of automated and manual steps. Each of the shift-stacked images was scanned for sources. The scanning was done with findsrc with a background window (WINDOW) of 20 pixels and a separation between the window and the target pixel (GAP) of 6 pixels. The detection threshold was set at 3 σ . With this setting, the number of false detections per stacked images was quite large - typically in the range of 3000 to 6000 per velocity stacked image and well beyond the ability to visually scan all candidates. To reduce the number of false positives, the pairs of visits were combined. A candidate worth examining was required to be seen on both visits with the same velocity and at a relative position within the stacks that was consistent with that velocity. A threshold on this test was set to 1.7 pixels multiplied by the ratio of the largest velocity to the mean velocity in the bin. The value  \nof 1.7 was empirically tuned to set the number of false positives to a manageable level. No threshold on consistent brightness was imposed at this step. For good data, there were typically ∼ 30 candidates per image per velocity, but the worst visit pairs had 10 times more.  \nThe visit-combined detection lists fed a tool that built a composite graphic image summarizing the detection candidate. This graphic contained small thumbnail images of the detections in each stack as well as the relevant portion of the original image. In all cases, the thumbnails were centered to the nearest pixel of the position of the candidate object. Also included was an aperture based magnitude estimate. Figure 7 shows the vetting graphic for (486958) Arrokoth. In this figure, the circle in the top left image is drawn at the measured location for the first visit stack. The top right image shows the stack for the second visit and the circle is drawn at the location that is predicted based on the stacking velocity. The size of the circles are comparable to the size of the positional coincidence imposed to find the candidates. The bottom of the figure contains two rows of five thumbnail images each. The top set of five is from the first visit and the bottom set is from the second visit. Each set of five is ordered in increasing time from left to right. The only processing done on the bottom images is the suppression of background stars. In this example, stars are seen as white streaks with some black pixels around the edge. Other than the object itself, the rest of the bright points are cosmic ray strikes that do not correspond to any object. This example shows the success of our spurious noise suppression algorithm. In the stacks there are a few low-level artifacts that are easily dismissed. For example, in the first stack there is a dark smudge to the lower right of the object. This is due to the residual of a star that was not perfectly subtracted. In both visits, there are other faint bright smudges but these are far enough away and at uncorrelated velocities so they are never considered. This graphic was designed to allow quick examination of all candidates for a determination of being real or not. In this case, there is clear evidence in all ten images for a source that leads to the two stacks. One of these ten (first image in second row) shows the object as being very close to a cosmic ray strike. The estimated apparent magnitude in both stacks is superimposed on the graphic below each detection.  \nFigure 7. Confirmation graphic for discovery of (486958) Arrokoth.  \n<!-- image -->  \nFigure 8 shows an example of a false positive that was discarded during the vetting process. The image stacks show an apparent source and is seen to agree with the velocity. The first clue against it being real is the lack of agreement in  \nthe magnitudes. This is apparent both in the numbers and visual appearance. The individual images provide the rest of the evidence against. In this case, the first visit alone is enough to discard the candidate. Image 1 has no source. Image 2 and 3 are on the edge of a subtracted star. Image 4 has a bright cosmic ray strike. Image 5 looks more like a background noise fluctuation. Similar logic applies to the second visit, Image 1 and 5 show nothing. In this case, there was clearly a confluence of spurious detections in each set of 5 that triggered a false positive detection. This combined graphic made it very easy to make a decision on the candidate.  \nFigure 8. Confirmation graphics for a typical false positive detection.  \n<!-- image -->  \nIt is perhaps possible to encode all of the logic used to grade candidate sources to avoid the manual step. No attempt was made to do this. The vetting process was already very fast and took much less than 10% of the total time to process an entire visit pair. Also, the time required to develop a fully tested automated tool was completely impractical on the timescale in which we needed to complete our search efforts.', '4.3.2. Cross-Checking Searches': "In order to cross-check our discoveries and to make sure no potentially real object was missed, our team compared the results of two independent search algorithms. The second algorithm worked with the HST 'drizzled' images which are both flat-fielded and distortion-corrected through the HST calibration pipeline. The intensity unit for these images is in electrons/s. Figure 9 records the steps of the second pipeline. Basically the full visit list of images for a pair of visits was read into a buffer. A star template was created by co-adding the images weighted by their exposure times, excluding high sigma pixels, and stacking on the positions of a selection of 6 user identified stars, ideally spread across the array. This process could be iterative if necessary. Once the star template for each visit was created it was subtracted from each individual image and a composite image for each visit was created for the family of KBO orbits provided by the rate correction previously discussed (i.e. if the family of orbit rates had 10 entries one would have 10 composite images for each HST visit to evaluate). The composite images for each pair of HST visits to the same position and at the same shift rate was then automatically searched for KBOs with image stamps of the discovered objects saved as .png files for visual confirmation. A display tool allowed one to blink-compare the results and the individual shifted and differenced images could also be brought up for searching or confirmation.  \n<!-- image -->  \n<!-- image -->  \nFigure 9. Sample images from our second reduction pipeline for the 2014 MU69 discovery image. (Top) Images from a single visit are shifted both by eye, then through an optimization routine to remove the stars. (Bottom) The star subtracted images are then searchable for KBOs in the field. The inlay on the bottom right is an extracted .png used for evaluation of reality.  \n<!-- image -->  \nIn addition to identifying potential targets via an algorithm, an independent manual-visual search was conducted in parallel. To systematically search each HST visit, a custom IDL program ('HSTlooker') was created to present side by side views of stacked and original images from corresponding visits, matching the predicted locations of an object from one visit to the next, and connecting the pixels of the stacked image to the data that created them. This tool drew from the prior work on searching the ground-based data. The HSTlooker program allowed users to cycle through the stacked and original images for each velocity shift and systematically cover the full frame to ensure no part of the image was missed. Each portion of the shifted and stacked images was visually scanned to identify any artifacts that persisted through median combination. Artifacts found in the stacks were verified via visual examination of the five original images to identify the source of the bright median pixel, be it a KBO or incidental alignment of cosmic rays, star subtraction remnants, scattered light or bad column. All five KBOs were independently 'discovered' via this process, albeit well after the objects were first located with the computer driven detection processes.", '4.4. Measuring KBO Positions': 'The next step in the processing was to obtain accurate positions for the validated objects. These measurements were always based on the original images, not the stacks. To start, a manually supervised aperture-based centroid was used. At this step, if the centroid looked reasonable for the position of the object, it was kept. If not, a manually derived position based on judgement of the position of the object was used. This aperture measurement is called a navigation position and fixes the region of interest for the next phase.  \nGiven a navigation location, a region of interest with a half-width of 9 pixels was defined. Within the region, an PSF-fit was performed on the object with the navigation position as its starting point. The PSF used was generated  \nby TinyTim (Krist et al. 2011). This location is saved separately from the navigation location. We used bad-pixel masks for the PSF fitting step to avoid contamination from cosmic ray strikes or residual blemishes from imperfectly subtract background stars. Most of the time the measurements were derived from the star-subtracted images but occasionally a better result was possible on the untouched data.', '4.5. Astrometric support catalog': 'The support catalog developed for the ground-based search turned out to be ideally suited for the HST data as well. There was useful overlap between the catalog and well-measured stars in the field such that there were generally ∼ 200 stars for a single detector in WFC3. The much improved image scale and smaller PSF enabled even more precise astrometry and the CFHT-based catalog held up amazingly well. Despite the quality of the catalog, the HST data clearly showed the limitations of the underlying astrometric catalogs. This led to the development of more and more intricate corrections to minimize the obvious systematic errors. None of these difficulties affected our ability to determine initial orbits or predict future positions for followup but it became clear that final targeting of an object would benefit from a better catalog. When the Gaia catalog was released in 2016 (Gaia Collaboration et al. 2016a,b; Lindegren et al. 2016), we switched to using that as the reference catalog for the CFHT data. We will return to a more in-depth discussion of our support catalogs in section 4.7.', '4.6. Measuring Star Positions and Updating the WCS': 'The world-coordinate system (WCS) header information (Gennaro et al. 2018) provides the means to transform pixel coordinates from the image to sky-plane coordinates. The WFC3 camera suffers from significant distortions and must be accurately known for precision astrometry. Kozhurina-Platais et al. (2016) did extensive work to map out the image scale and focal plane distortion. This map provides relative calibration down to the level of 0.01 pix (0.15 µ m) or better between any pair of pixels in the image. The absolute pointing information comes from the guide stars observed with the fine-guidance sensors (FGS) that are providing pointing and guiding control. These stars are known to have typical errors of 0.1-0.2 arcsec but sometimes are significantly worse than that. As a consequence, the absolute pointing information is of a much lower quality and must be improved.  \nTo start with we measured the stars that are in the support catalog with PSF fitting. Again, we used the TinyTim PSF appropriate for the location of each star on the detector. Still, the star fitting poses one additional challenge due to tracking of the KBO cloud. All star images show significant trailing which changes on every exposure. The dominant apparent motion of the stars is due to parallax as HST moves around the Earth in its orbit. Even though the images look like linear trails most of the time, the trailing is actually a time-dependent quantity and the residence time of a star at any given point on the trail varies throughout the exposure.  \nTo address this time-variable trail, we modified the PSF fitting process by introducing a convolution kernel into the model image production. The routine jitterk.pro takes the time-dependent position and computes the necessary kernel. More details about this tool can be found in Appendix B. The time-dependent position is derived from the HST-centric position of the guiding orbit that represents the center of the cloud of candidate objects. Also required was the definitive post-observation HST orbit provided by JPL 2 .  \nOne minor limitation of our tools is that the image distortion was not taken into account for the star smear kernel. The kernel is computed for the center of the array and then used over the entire image. There are very subtle shifting features in the residuals from the fit across the field that look like the model image is trailed at a slightly different angle from the image. We looked for patterns in the astrometric analysis that could be traced back to this approximation but no convincing signatures were found. Once the new Gaia catalog (Brown et al. 2018) is included the pipeline WCS can be dramatically improved but this will only ever be constrained by the position of two stars. By using more stars in the field, the solutions can always be improved.  \nOur final WCS update consisted of solving for a translation of the reference position in the WCS stored as CRVAL1 and CRVAL2 (Gonzaga et al. 2012) and a small angle rotation to the CD matrix. All other components of the WCS remain untouched. In practice, the rotation corrections are extremely small and amount to worst-case corrections in the image corners that are much smaller than a pixel (a few milliarcsecs at most). The apparent offsets in our data are larger than one would expect from just the errors in the guide star catalog but this is due to the bookkeeping of the  \npointing information for data taken with moving target tracking. Our pointing solution update removes all of these effects.  \nWe received a report of errors in the reconstructed trajectory of HST deduced from astrometric observations of C/2013 A1 (Siding Springs) relevant to its close passage to Mars (Farnocchia et al. 2016). In that work, the internal errors from the HST observations did not agree with the post-fit errors when combined with other sources. In that work the HST data were de-weighted to 0.3 arcsec due to presumed errors in the HST position at the level of 10 km. If this interpretation is correct, we should presume this same error is present in our data. However, our targets are nearly 30 times further away and perhaps that reduces the problem. Taken at face value, this issue indicates HST positional errors could contribute as much as 10 mas to the total error budget. However, in light of our successful occultation work on Arrokoth (Buie et al. 2020), this level of error in our astrometry is just not seen and we have no evidence in our data and reductions for any problems with the reconstructed positions of HST that supported our astrometric work.', '4.7. Final Astrometry on KBOs': 'The KBO position for each image is derived using the methodology described in Porter et al. (2018) (sections 2.3-2.4) and further discussed in Porter et al. (2022). In short, a MCMC algorithim is used in combination with Tiny Tim PSFs to determine the optimum KBO pixel position PDF which then has the WCS applied to extract the R.A./Dec sky-plane position. All measurements are provided in J2000 (EME200). Table 8 lists the positions for all objects measured in the HST survey.  \nTable 8 . HST Astrometry  \nTable 8 (continued)', '4.8. Targetability Assessment': 'Our HST discoveries are listed in Table 9. The apparent brightness measured is shown along with the positional data. The photometry tabulated has a very large uncertainty and the variations in the measurements is a consequence of the low signal-to-noise ratio and not any intrinsic variability in the objects. In the case of Arrokoth, a much more  \nBuie et al.  \nTable 9. HST List of Potential TargetsTable 10. Orbital Elements for HST Discoveries  \nNote -The values in italics are just estimates based on limited data. No additional data were ever acquired to improve these estimates.  \nextensive dataset was collected and reported by Benecchi et al. (2019) but no rotational modulation could be identified. The post-flyby results from New Horizons showed that our vantage point from the Earth is essentially looking down the rotation pole thus explaining the lack of lightcurve. In fact, we suspect that this geometric circumstance helped us find Arrokoth since it maintained a consistent brightness just above the detection limit of HST. A strong lightcurve would have likely led to missing followup observations and complicating the orbit determination work.  \nThe astrometric measurements in Table 8 were used as time progressed to produce orbit fits and to evaluate the likely ∆V needed to re-target New Horizons for a close encounter. On this basis alone, the candidate object 4510067S (PT5) was immediately discarded as out of reach and no further followup observations were attempted. For the remaining candidates, we obtained additional observations to improve the orbit estimates (GO-13663, GO-13311). Objects 0720090F (PT4, 2014 MT69) and E31007AI (PT2, 2014 OS393) were eventually revealed to also be beyond reach, the latter became a distant encounter object. The objects G12000JZ (PT3, 2014 PN70) and 1110113Y (PT1, 2014 MU69) survived as potential candidates into 2015. Of the two objects 2014PN70 was slightly brighter, implying a larger object, but the required ∆V was at the very limit if the remaining resources and would leave very limited margins. In the end, our choice was (486958) Arrokoth (formerly 2014 MU69) which required the smallest trajectory maneuver and left ample fuel reserves for other contingencies. Table 10 summarizes the final orbital elements for the three objects that were carried forward from this search based on all data collected through the end of 2016 with estimates of their uncertainties.', '5. DISCUSSION/LESSONS LEARNED': "In summary, our combined ground-based efforts from Subaru and Magellan resulted in the discovery of 80 new KBOs. An additional 5 KBOs came out of our HST search including the fly-by target (486958) Arrokoth. Eleven of these objects (excluding Arrokoth) became New Horizons long-range observation targets which received measurements at phase angles far greater than can ever be observed from earth-orbit locations. These observations are allowing us to map surface scattering properties of small KBOs for the first time (Verbiscer et al. (2019),Verbiscer et al. (2022)). The  \ntools we developed using the excellent Gaia astrometric catalog have also allowed us to make occultation predictions for many of these objects at high precision (Porter et al. (2022)). The successful Arrokoth occultations were an excellent validation of this approach (Buie et al. 2020).  \nWe also gained significant insight into searching for KBOs in general, but in the galactic plane specifically. When we started we expected that a target search would be a relatively straight forward process resulting in multiple potential mission objects to select from or prioritize among. Many on the search team had previously been part of dedicated KBO search projects from the time of the Kuiper Belt's observational discovery and we thought we had a pretty good handle on the Kuiper Belt population density and distribution.  \nOur initial 2004-2005 survey didn't find a target, but it served a critical role in learning just how challenging it was to work in extremely crowded galactic plane fields. It also provided sample datasets for developing tools and techniques to deal with identifying faint objects in high density star fields and was thus very successful in moving us forward in the search. These lessons learned were carried forward into the second 2011-2014 survey epoch and helped us to optimize both data collection and level of effort or image analysis. In particular we learned that image quality, both observing conditions and seeing stability, was critical. Observations at high ( > 1.6 airmass) or under non-photometric conditions were insufficient to meet the search goals. We also quickly learned the importance of having an astrometic star catalog that went deeper than the standard databases of the time (USNO B1.0). Our solution was to use MegaCam on the CFHT from 2012-2018 to build a deeper (fainter magnitude) high precision astrometric star sub-catalog specific to our search fields. Stephen Gwyn built these catalogues from a variety of references, building off various Gaia (DR2) and 2MASS references sets for the global frame and then internally building up to measure high precision relative positions. When the official Gaia DR2 version became available we re-calibrated the final catalogue. This took a tremendous amount of effort, but using this secondary reference sub-net made a huge difference in the earlier efforts and was essential for the second and third epoch searches.  \nOur search effort was also fairly early in the timeline of large data volumes among astronomical projects. Due to the volume of information, we quickly learned the critical importance of having effective pipelines that could process the data inclusive of data quality metrics on image characteristics, performance etc. so that one could quickly identify if quality data were being collected and also how prioritize data analysis. These lessons have directly informed searches for additional objects following the New Horizons fly-by of Arrokoth. The new search epoch, 2019 and following, as a result, are far better designed and characterized (Fraser et al. 2024). Likewise, it became clear by the end of 2013 that the likelihood of the data we were collecting at the time for identifying a fly-by object was minimal at best and a new strategy needed to be employed, hence a search with HST.  \nDesign of the HST search was itself a critical point in the discovery process as it required a full review of the state of knowledge of the Kuiper Belt at the time, in 2014. A 10-year improvement was made in our understanding of the absolute density and extent of the region in addition to search technique and analysis gains. We also accurately estimated what we would be able to find using HST. Finding 2 objects that we could consider targeting out of the sample of 5 detected was in fact an expected outcome rather than an under or overestimate. It is safe to say that by the time the survey was completed we knew the types of targets we were looking for, and had a good handle on the systematics, characteristcs and processes required for success.  \nLooking to the future, searching for additional KBOs, this work has led us to push our techniques even further into building deeper images stacks, applying machine learning, and considering additional (non-classical) KBO populations. These details are left to a future publication. However, it can't be overstated that the work presented here was a critical step along the pathway for identifying a fly-by object for New Horizons and in designing, characterizing and executing these sorts of search surveys while also sharing the methods and tools along the way.", '6. ACKNOWLEDGMENTS': "This work is dedicated to the memory of Peter Collins (1949-2017), instrumental to much of the software for the 2004 search that began this entire project.  \nBased on observations made with the NASA/ESA Hubble Space Telescope, obtained at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555. These observations are associated with programs 13633, 14053, 14092, 14485, and 14629. We are greatly indebted to Mike Ryschkewitsch (John's Hopkins Applied Physics Lab), Jim Green (Director of NASA's Solar System Exploration Division) and Ken Sembach (Space Telescope Science Institute) for their efforts in helping us to acquire the time for these observations.  \nThis paper includes data gathered with the 6.5 meter Magellan Telescopes located at Las Campanas Observatory, Chile. It is also based on data collected at the Subaru Telescope and retrieved from the HSC data archive system, which is operated by Subaru Telescope and Astronomy Data Center at National Astronomical Observatory of Japan. Access to the Subaru telescope was supported through NOAO Keck-Subaru exchange time, NASA HQ support and direct support from Taiwan through the Academia Sinica Institute of Astronomy and Astrophysics. The Subaru Telescope is operated by the National Astronomical Observatory of Japan. Likewise, the authors wish to recognize and acknowledge the very significant cultural role and reverence that the summit of Maunakea has always had within the indigenous Hawaiian community. We are most fortunate to have the opportunity to conduct observations from this mountain.  \nThis research used the Canadian Astronomy Data Centre facilities operated by the National Research Council of Canada with the support of the Canadian Space Agency. This work also made use of the Gemini Observatory Archive, NASA's Astrophysics Data System Bibliographic Services, and the JPL HORIZONS web interface (https: //ssd.jpl.nasa.gov/horizons.cgi). Some of the follow-up astrometry to refine orbits of newly discovered objects to access object targetability by New Horizons were obtained at the international Gemini Observatory, a program of NSF's NOIRLab, which is managed by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation on behalf of the Gemini Observatory partnership: the National Science Foundation (United States), National Research Council (Canada), Agencia Nacional de Investigaci'on y Desarrollo (Chile), Ministerio de Ciencia, Tecnolog'ıa e Innovaci'on (Argentina), Minist'erio da Ciˆencia, Tecnologia, Inova¸c˜oes e Comunica¸c˜oes (Brazil), and Korea Astronomy and Space Science Institute (Republic of Korea).  \nThe Zooniverse.org platform, development of which is funded by generous support, including a Global Impact Award from Google, and by a grant from the Alfred P. Sloan Foundation, hosted and promoted the 'Ice Hunters' project with all of the data handling and analysis being prepared by the authors of this paper.  \nThis work was also supported by a Harvey Mudd College Clinic Program team in 2009, consisting of faculty member Patrick Little, and team members Austin Lee, Claire Robinson, Florian Scheulen, Steven Berry, Chris Sauro, and Cullen McMahon. We also thank James Miles of Boulder, Colorado for his volunteer work in support of software development, observing support, and data processing.", 'A. ASTROMETRIC DESCRIPTION OF AN IMAGE': 'The astrometric tools for this project pre-dated any development of a FITS standard for the astrometric coordinate description of an image that required a non-linear transformation. The current standard system is described in Greisen & Calabretta (2002) and Calabretta & Greisen (2002). As a result, Buie developed a system for this purpose. This appendix describes the system and how it is encoded into a FITS header. As the FITS standard was developed, the system described here was adjusted to be complementary to any existing system by using a set of non-overlapping keywords so that both can co-exist in the same header. This alternate system was designed to be used with software that can modify an existing FITS header and place these keywords in the header (an operation that is now very commonplace).  \nThe task at hand is to describe a transformation from raw coordinates internal to some image or digital detector to coordinates on the celestial sphere, that is, right ascension and declination. This process involves the four coordinate systems indicated below,  \n( x, y ) ⇔ ( x m , y m ) ⇔ ( ξ, η ) ⇔ ( α, δ ) (A1)  \nwhere ( x, y ) denote the native coordinate system implied by the image array or detector pixels, ( x m , y m ) are the monolith coordinates, ( ξ, η ) are the tangent plane coordinates, and ( α, β ) are the right ascension and declination on the sky. For this problem one needs to be able to convert between any coordinate system as indicated. Note that these are considered to be linked in the order shown. Thus, to compute tangent plane coordinates from raw pixel positions you first convert to monolith coordinates and then to the tangent plane. In the case of a single detector focal plane, the monolith and native coordinate systems are one and the same. A monolith coordinate system is a virtual pixel scale that mimics the native coordinate system of a chip but this single system extends to include the area of all of  \nthe detectors. Typically, detector coordinates are never negative and the implementation of this system uses (0,0) for the coordinates of the first pixel in the image and the coefficients and IDL software reflect this design choice.', 'A.1. Tangent plane and celestial sphere': "The transformation between the tangent plane ( ξ, η ) and celestial sphere ( α, δ ) can be found on p. 283 from Smart & Green (1977) and is reproduced below:  \nη = tan δ -tan δ 0 cos( α -α 0 ) cos( α -α 0 ) + tan δ tan δ 0 and (A2)  \nξ = sin( α -α 0 ) cos δ (cos( α -α 0 ) + tan δ tan δ 0 ) , (A3)  \nwhere α 0 and δ 0 are the coordinates of the tangent point. These expressions have been re-written from those found in the reference to minimize the impact of mathematical singularities in the equations and to speed their execution on a computer. Despite the re-write, take note that equation A3 still contains a singularity at the celestial poles ( ± π/ 2). The inverse transformation is written,  \ntan( α -α 0 ) = ξ cos δ 0 -η sin δ 0 , and (A4)  \ntan δ = cos( α -α 0 )( η +tan δ 0 ) 1 -η tan δ 0 . (A5)  \nNote that care must be taken when evaluating equation A4 so that all 4 quadrants of the arctangent are preserved. Equations A2 and A3 are evaluated by the IDL procedure astrd2sn.pro . Equations A4 and A5 are evaluated by IDL procedure astsn2rd.pro . Both transformations are included in astcvt.pro . By virtue of these transformations, one can treat tangent plane and sky coordinates as completely interchangeable and the discussion below will often make little distinction between the two. Tangent plane coordinates are also referred to as 'standard coordinates.' The natural units of calculating the tangent plane coordinates would be radians. However, it is usually more convenient to work in units of arcseconds instead.", 'A.2. Image and tangent plane': "A plane tangent to a sphere at a given point is also known as the tangent plane. This plane is very similar to the physical mapping that takes place when a camera takes a picture of the sky. For small enough angular range images the mapping is linear. Departures from non-linearity are usually a consequence of optical distortions in the camera optics but also inevitable with a large field-of-view. The mapping from detector (image) pixels to the tangent plane is thus an approximate process that needs to be adjusted to match the characteristics of the image so that a sufficiently accurate mapping is achieved.  \nDigital image data collected with a single device is a common type of data and represent a simple initial case for consideration. The mapping between these two coordinate systems can be accomplished by any set of basis functions capable of reaching the desired accuracy. A common and sufficient basis set is an n -degree polynomial in x and y. In its simplest, linear form the transformation would look like this:  \nξ = a 0 + a 1 x + a 2 y (A6)  \nη = b 0 + b 1 x + b 2 y (A7)  \nwhere the a n and b n coefficients determine the transformation. This approach is just a power series expansion in x and y . Typically, the coefficients are determined from some least-squares fit between the image and a stellar catalog. For the sake of this discussion the method of determination is not important. A linear expansion works well for sufficiently small detectors but the systems we used for this search required higher order terms and the IDL-based tools support up to a fifth order expansion. Some of the original tools that were used as a starting point for this development used an r term ( r = √ x 2 + y 2 ), but this term is not part of the polynomial basis set and is incorrect to use in a cartesian framework. Support for this term is included in the IDL tools, but is present for historical reasons and should not be used. Other basis sets are possible, such as Chebyshev and Legendre polynomials. Thus, to specify this transformation it is necessary to indicate which basis set is to be used and then which terms are provided along with the coefficients.  \nThe IDL routine astcvt.pro supports all three of these basis sets. Note that this methodology does not require or enforce orthogonality between the ξ and η coordinates.  \nIf the transformation is non-linear it becomes important to specify the location of the optical axis. This point is special in that it will represent a point of symmetry in the mapping. If the device is well fit by a linear transformation it is common to use the center of the device as the optical center. In this case the center ( x c , y c ) is defined to be ( n x / 2 , n y / 2) where n x and n y are the x and y size of the image, in pixels. In non-linear cases it works best to use the optical axis position, wherever that might be. While it is important to know the optical axis position, the determination of this via a least-squares fit often leads to a poorly constrained answer. Knowledge of this position will usually give a better answer than relying on a fit but a precise determination is not needed.  \nIt is often necessary to rescale x and y prior to solving for coefficients. For a linear polynomial fit this rescaling is optional. For non-linear polynomials, the rescaling also helps make the solutions more robust even for a simple polynomial basis set. For Chebyshev and Legendre polynomials, the rescaling is required since the input range of values must be between -1 and 1. Normally this is done by dividing by a renormalization factor, d , such that x ' = x/d and y ' = y/d . Renormalization is usually done by using the size of the image array or detector. One form of normalization is d = √ n 2 x + n 2 y . However, this can lead to an incomplete filling of the allowed range of x ' and y ' between ± 1. A better choice is d = max(n x / 2 , n y / 2) if the center is used for the optical axis position. Otherwise, use the maximum of the absolute value of the corners of the array after subtracting the optical axis position. In the IDL software, the renormalization is carried only in the monolith coordinates.", 'A.3. Image mosaic conversions': "A closely related problem to a single image is dealing with a mosaic of images from an array of detectors within a single camera. They are also related though a presumed constant physical bond between detectors. Processing and fitting individual images from a mosaic can, in principle, be handled as if they were each independent images subject to the treatment in the previous section. Unfortunately, such handling is subject to numerical instabilities in the fitting process leading to often erroneous coefficients. For the case of a mosaic it is still true that there is a native raw coordinate system for each device that is identical to the ( x, y ) system in the previous section. This case benefits from a new coordinate system intermediate between ( x, y ) and ( ξ, η ) that we will call the monolith coordinates denoted by the subscript m . A well-constructed mosaic camera is one where all of the imaging devices lie in the focal plane of the instrument. In such cases, the raw coordinates of each device is related to the other devices by a translation and a rotation. The following treatment explicitly ignores any out-of-plane tilt a detector might have. The choice of the monolithic coordinate system is arbitrary, so for convenience we choose one of the devices to be co-joined with the monolithic system. It is also useful to choose this one device to be the one that contains or is at least closest to the optical axis. In this coordinate system the position of the optical axis is defined to coincide with the origin of the monolithic system. The transformation from the i th raw device to monolith coordinates is thus,  \nx i m =(( A i 1 , 1 ( x -x i c ) + A i 1 , 2 ( y -y i c )) /d (A8)  \ny i m =(( A i 2 , 1 ( x -x i c ) + A i 2 , 2 ( y -y i c )) /d (A9)  \nwhere A is a rotation matrix related to a rotation angle θ by  \nA = ( cosθ sin θ -sin θ cos θ ) (A10)  \nNote that the elements of A are often referred to as the 'CD matrix' found in FITS headers for WCS (world coordinate system) conversions (Greisen & Calabretta 2002; Calabretta & Greisen 2002). The values for ( x 0 , y 0 ) are often found in the CRPIXn keywords. The choice of the reference chip leads naturally to θ = 0 for that image. With this construction, the x i 0 values are the raw coordinates of the optical axis for each image. While it is useful and convenient to impose no rotation between the monolith system and the anchor chip, there are cases were a global rotation is useful. The rotation for the anchor chip will impart that global rotation if applied at the right step. In our IDL software, astcvt.pro , this global rotation is supported.", 'A.4. Inverting an astrometric solution': 'The equations above make it very easy to go from position ( x, y ) to sky-plane coordinates, either ( α, δ ) or ( ξ, η ). However, it is often very handy to be able to go in the other direction. That is what astcvt.pro does. Note that  \nthe convertion between ( α, δ ) and ( ξ, η ) in either direction is very simple and is not discussed here. If the solution is linear the conversion is deterministic and exact via an analytic expression. If the solution is non-linear the conversion requires an iterative numerical solution. The routine astcvt will do this inversion for any input solution.', 'A.5. FITS keywords': 'All of the keywords in the Buie system start with AST . The variables expressions shown refer back to the quantities defined above. The units of the variables are provided in brackets.  \nASTINFO (boolean) - Flag to indicate this header contains the astrometric information in this format.  \nASTXCEN (float, x 0 ,[pixels]) - X pixel location of the optical axis.  \nASTYCEN (float, x 0 ,[pixels]) - Y pixel location of the optical axis.  \nASTPROT (float,[radians]) - Chip pre-rotation angle.  \nASTRNORM (double, d ,[pixels]) - Normalization factor.  \nASTTRA (double, α 0 ,[HMS string]) - Right ascension of the tangent point. It is not required, but highly desireable that this point be coincident with ( x 0 , y 0 ).  \nASTTDEC (double, δ 0 ,[DMS string]) - Declination of the tangent point. It is not required, but highly desireable that this point be coincident with ( x 0 , y 0 ).  \nThe following three-keyword set is given for every term in the basis set. The end of the keyword is a numeric digit, starting at 1 and increasing for as many terms as are required for the mapping. The order of the terms is arbitrary.  \nASTTNi (string) - Name of the term. The names indicate the term and the basis set to be used. There are three basis sets supported by astterms.pro : simple polynomial, Chebyshev degree 1 or degree 2 polynomials, and Legendre polynomials (M=0). The recognized term names are CONST , X , Y , XX , YY , XY , XXX , YYY , XYY , XXY , XXXX , YYYY , XYYY , XXYY , XXXY , XXXXX , YYYYY , XYYYY , XXYYY , XXXYY , XXXXY , or all of the preceeding with a T prepended for degree 1 Chebychev polynomials, a U prepended for degree 2 Chebychev polynomials, or P prepended for Legendre polynomials. The order of the term is indicated by the repetition of X or Y . For example, the x 3 term is denoted with XXX .  \nASTCXi (double, a i , arcsec) - The transformation coefficients from ( x, y ) to ξ . The units are in arcsec raised to the degree of the term.  \nASTCEi (double, b i , arcsec) - The transformation coefficients from ( x, y ) to η . The units are in arcsec raised to the degree of the term.', 'A.6. Useful Buie IDL library tools': 'The preceeding sections describe the mathematical basis and algorithmic approach needed to extract astrometry from an image. This approach can, of course, be implemented in most programming languages. The rest of this section describes our specific implementation of tools in the IDL programming language. The following software tools are in direct support of handling the operations related to astrometry and can be used directly or to guide development and testing of other versions. The top level tools are the ones most often used to build applications to extract or use astrometry. The low-level routines are usually just used internally in the library. The mid-level routines are somewhere in between the other two categories.', 'A.6.1. Top-level routines': 'astcvt - Convert from one astrometric coordinate system to another given an astrometric description of the image. The coordinate systems supported are raw pixel positions in the image (xy), detector mosaic monolithic coordinates (XY), tangent plane coordinates in radians (sn), tangent plane coordinates in arcsec (SN), and equatorial celestial coordinates (rd).  \nastinfo - Read (or write) astrometric information from (to) the header.  \nmkastinfo - Create a synthetic astrometric description. This is intended only for creating linear descriptions.', 'A.6.2. Mid-level routines': 'astchi1 - Compute the goodness of fit for a solution based on some rotation and offset. This is used internally to help solve for inter-chip rotations in a multi-detector camera.', 'A.6.3. Low-level routines': 'asteval - Evaluate an astrometric polynomial function given the independent variables ( x, y ), a list of terms, and the polynomial coefficients.  \n- astrd2sn - Convert from ( α, δ ) to ( ξ, η ). This routine is called by astcvt but is useful on its own.\n- astrd2xy - Convert from ( α, δ ) to ( x, y ). This routine works only for linear transformations. astcvt is almost always preferred.\n- astsn2rd - Convert from ( ξ, η ) to ( α, δ ). This routine is called by astcvt but is useful on its own.\n- astsn2xy - Convert from ( ξ, η ) to ( x, y ). This routine is called by astcvt and is rarely called directly.  \nastsolve.pro - This routine performs the least-squares fit for an astrometric solution. The inputs are correlated lists of pixel coordinates ( x, y ) and tangent plane coordinates ( ξ, η ) along with a list of fitting terms to use. This tool performs a robust fit, automatically culling discrepant positions based on the noise statistics in the fit. Considerable work is required before-hand to set up the correlated lists but this tool does the basic work of the fit. One example of an interactive tool for solving for astrometry that feeds this solver is astrom.pro . This one example is a heavily used tool for other projects but the scale of the data reduction process for the New Horizons search was such that it was more efficient to write dedicated and very non-general solving tools.  \n- astterms - Evaluate the independent vectors to match an astrometric polynomial function.\n- astxy2rd - Convert from ( x, y ) to ( α, δ ). astcvt is usually used rather than this tool.\n- astxy2sn - Convert from ( x, y ) to ( ξ, η ). This routine is called by astcvt and is rarely called on its own.', 'B. DETAILS OF KEY BUIE IDL LIBRARY ROUTINES': 'The following routines were particularly useful in the data processing and results presented in this paper. Many more were also used but these contain either the essential computational engines or interesting and perhaps non-obvious algorithmic concepts. This discussion is provided to highlight the algorithms used rather than to be a guide for the use of this specific implementation.', 'B.1. basphote.pro': 'Aperture photometry summation and centroiding is provided by basphote.pro . This routine and its use has been discussed in Buie & Bus (1992) and Buie (1996). Alternate versions are available in C and were also ported to IRAF many decades ago. The object aperture integral is the sum of all pixels weighted by the exact area of overlap between the pixel and the circular aperture. Sky signal is determined from a sample of pixels whose centers lie between the inner and outer radii of the sky annulus. The robust mean of the sample ( robomean.pro ) is computed and the sky mean and the sky noise is used to remove sky signal from the object aperture summation and to propagate into the final noise estimation. Detector gain, exposure time, and detector read-noise can also be provided for more accurate and consistent flux and noise computations. During the measurement, the input position is treated as a rough location. The image is scanned for the maximum within the distance of the object aperture radius. From this position, an initial centroid is calculated without subtracting sky. This first centroid is then used to place the sky annulus, determine sky, then compute the object aperture flux and final centroid.', 'B.2. dewarp.pro': 'This tool was created in a collaborative project with a Harvey Mudd College clinic team. The basic function of this tool is to map image data from one grid to another. The input is an astronomical image and an image descriptor that contains an astrometric solution for the image. Another input is a second astrometric solution for an image to be constructed. The two solutions are used to map flux in the input image to the correct location in the output image. The images are resampled, guided by the astrometric solution, using an interpolator that is equivalent to a sinc interpolation.  \nThis tool had many uses in our project. It can be used to map one image onto another eliminating any distortion, sampling, or rotation from the resampled image and then replicating any distortion in the destination image. This type of operation is particularly valuable when building a deep stack of the sky that matches a single image in preparation for subtracting that template without sampling the original image.  \nFor the ground-based searches it was more advantageous to define a synthetic, linear image that all search data mapped to for a particular field. In this case, all images are resampled but the result is strictly linear and registered. Given that the image transformation is always done with full astrometric knowledge, the resampled images are just as valid for astrometric measurements of discovered objects as the originals.', 'B.3. findsrc.pro': 'Automatic source detection is important element of getting astrometry and photometry from a digital sky image. This tool is similar in function to the DAOPHOT tool, SExtractor (Stetson 1987; Bertin & Arnouts 1996). However, the internal algorithmic workings of this program are quite a bit different and its output can have very different characteristics. In particular, the number of false positives is usually much lower.  \nThe first step in this tool is a scan of the image to find pixels that exceed a given threshold SIGTHRESH relative to the noise in the background of the image. For each pixel, four windows on the background are computed. These windows extend away from the pixel in the cardinal directions of the image ( ± x and ± y ). The length of each window is controlled by the WINDOW keyword and defaults to 6 pixels. The windows are displaced by GAP pixels from the target pixel with a default of 2 pixels. This value is normally set to a value similar to the FWHM of the PSF in the image. From these windows, four sky signal values are computed with an associated standard deviation. If the pixel signal exceeds the threshold above the sky in 3 or 4 directions, that pixel is flagged. The result of this step is a detection array with a binary value for the pixel exceeding the threshold or not. Note that an additional check is made to restrict the detected pixels to have a signal level less than the saturation limit imposed by the MAXPHOTSIG keyword.  \nThe next step is to collapse the detection array. This is needed since many pixels from the same source can be flagged as detections. This step is performed by an internal rountine, collapse . The detection image is eroded so that only one pixel is left flagged per source and that pixel is at the local maximum in the image. The local maximum are further constrained to be at least GAP pixels apart. The erosion logic is designed to also collapse large regions from saturated stars and minimize the number of pixels flagged as detected. This step is reasonably effective but local variations along bleed trail from a saturated source can have local maxima that will confuse this tool and lead to a string of detections along the trail.  \nThe eroded detection image is then a map of source locations. Aperture photometry is computed with basphote.pro for each flagged pixel in this map. The radius of the object aperture defaults to GAP but can be altered with the OBJRAD keyword if desired. The sky annulus for a local sky computation is set to objrad+1 and objrad+5. Those sources with sensible values for FWHM and instrumental magnitude are kept. Also, the SNR of the source must also be above the value of SIGTHRESH . During the photometry calculation, a centroid location is computed as well.  \nAt the end, the centroid position, FWHM, flux, magnitude, magnitude error, SNR, sky and sky noise are returned to the caller or possible saved to a file for later processing. This routine is reasonably fast with the most expensive part being the computation of the aperture photometry and centroids.', 'B.4. frmdxdy.pro': 'One of the most important tasks during the determination of an astrometric solution is the correlation of two lists of positions. In this case, we need to take a list of star positions, usually in ( x, y ) coordinates and determine the correspondence into a list of star catalog positions, usually in ( α, δ ). The algorithm used here was inspired by a method developed by Bruce Koehn at Lowell Observatory for the LONEOS project. Approximate information about the image scale and orientation is used to convert the catalog positions to estimated pixels positions. This tool is extremely resilient against differences in photometric depth between the two lists. Rough knowledge of the plate scale is needed. The position shift between the two lists can be large with no penalty on speed of execution or quality of result. The rotation of the field, however, must be known very well. The requirement on this knowledge gets more stringent as the field-of-view increases. This limitation can be ameliorated somewhat by restricting the initial determination of shift to a subset of a large image. The presence of uncharacterized non-linearities in the image also limits the ability to find a good shift but images subsets help here too.', 'B.5. frmdxyr.pro': 'This routine is an attempt to fix the problem of needing to know an exact rotation angle for frmdxdy . Rotation is scanned over a range and a binary search is run to find the best rotation angle that gives a sensible list correlations. This tool can be very time-intensive to run if the rotation range is large and is really only a stop-gap measure when there is small rotation error in the datasets. This tool is completely useless for the class of problem where the rotation is not known at all. For most astronomical data, one usually has a pretty good idea of the information for these routines.', 'B.6. jitterk.pro': "PSF-fitting to an image to get the position and flux of a source requires a good understanding or characterization of the PSF. For stellar data taken with stellar guiding, knowing the PSF is sufficient and is also easy to determine from the data (see psfstack ). For solar system observations, the sources of interest are moving. The situation is even more complex for Hubble Space Telescope data where parallax in introduced.  \nThis tool was specifically written to support HST data reductions though it can be usefully used on ground-based data as well. The input is a time-history of the motion of an object relative to the pixel coordinates during an exposure. In the case of the HST search data for this project, we track the KBO. As a consequence, the stars are trailed by the target motion and parallax induced by the telescope's motion around the Earth. In principle, one can also add in the measured jitter of HST pointing as measured with its fine-guidance sensors but this level of detail was never needed.  \nAt first glance, most of the trailed HST images look like linear trails. In practice, this visual appearance is misleading. The parallactic motion can still be along a line but sometimes, the rate of motion is slower at the start of the exposure than at the end or vice versa. If the images are fitted as lines and the line center used, there can be an effective temporal shift in the derived astrometry. This tool builds a time-history of the pointing that is used to create a convolution kernel that can be applied to the PSF prior to the fitting process. This ensures that the temporal registration between smeared and un-smeared targets is preserved.", 'B.7. lowess.pro': 'This tool is a LOcally-WEighted Scatterplot Smoother. It provides a means to draw a smooth curve through data to follow the trend in the data even in the presence of noise without imposing a smoothing or fitting function.', 'B.8. ois.pro': 'This tool was originally written by Miller et al. (2008) as a port of the Alard & Lupton (1998) optimal image subtraction (OIS) software into IDL. That tool was very slow and not very practical to use. Our version was optimized during a Harvey Mudd College clinic project and runs at essentially the same speed as the original Ccode. An additional enhancement was added by Miller et al. (2008) to encode a δ -function basis set for the PSF manipulations rather than the nested gaussian approach from Alard & Lupton (1998). During our extensive work on the ground-based data for this project, we found the δ -function basis set, coupled with very precise astrometric registration eliminated all need for spatially variable convolution kernels in the OIS process and that feature of the software was never used.', 'B.9. robomean.pro': 'Numerous operations in our data require the computation of a mean from a sample of data. Most real data unfortunately are affected by samples that do not fit the underlying noise distribution of the quantity being measured. This tool is designed to provide a robust mean from a sample is is probably the most heavily used routine from the Buie library. The goal is to iteratively process the data until all of the surviving points are consistent with the noise distribution in the surviving points. A threshold can be set on the filtering and is normally set to 3 σ . Convergence depends on getting convergence in the value of the mean, standard deviation, skew, and kurtosis consistent with the errors on these four quantities as points are (possibly) removed. This is used for measuring sky brightness around a source for position and flux measurements as well as just about any other place where a mean is desired.', 'B.10. psfstack.pro': 'Many operations require a PSF for an image. This routine determines a numerical PSF by stacking many sources from the images to get an average PSF. An input list of sources is provided and they are stacked in descending order of SNR down to a limit that defaults to 25. Sources that are not good PSFs (such as a double-star blend) are excluded.', 'B.11. skyfit.pro': 'The degree of image crowding in the image made it difficult to measure the signal level due to the light from the sky. A typical method is to collect a large sample of pixels across the image and then compute a robust mean, removing unusual outliers caused by such things as pixels illuminated by stars. This method requires the sample to be dominated by actual sky measurements. Levels of contamination up to 10-20% are easily removed but when that fraction rises too high our usual tool failed and returned an estimate of the sky that is systematically and significantly too bright.  \nThis new routine includes robust estimation and simultaneous fitting of a low-order 2-D polynomial to the sky background with one important enhancement. An examination of a histogram of the random image pixels shows an asymmetric distribution. The appearance of this histogram in Fig. 10 reveals what looks like a normal distribution for signal at or below the mode of the distribution. Signals higher than the mode are contaminated by faint star signals that are just barely brighter than the sky and thus add flux to the apparent distribution of signal from the sky. If there were no stars, the distribution will always look like a symmetric Gaussian but as the number of star-illuminated pixels increases the distribution becomes more and more distended on the upper half. Pixels below the mode are minimally affected by background star light.  \nFigure 10. Sky signal example. The plotted histogram shows the residuals from a 2-d fit to the sky in a single image. The thick orange curve is a lowess curve of the histogram and the peak is used as an estimate of the mode. The vertical green lines mark the mode and the -1 σ point on the low side of the distribution. The gray shaded region shows the approximate Gaussian representation of the sky signal and noise.  \n<!-- image -->  \nIn detail, the program first selects 20,000 random pixels from the image (2k by 4k array size) representing 0.2% of the pixels. A first-order 2-D polynomial is then fitted to these pixels. Those pixels whose residuals are more than 3 σ from the mean residual are rejected. If there are new rejected pixels, the fitting process is repeated. When there are no additional pixels removed with up to 10 passes, the 2-D fitting is deemed complete. The mode of the surviving pixels is computed and then the 1/e half-width is measured from a smoothed curve based on the data using lowess.pro applied to data falling below the mode. The true sky signal is assumed to be normally distributed and we use the mode and 1/e half-width to approximate the sky signal and its uncertainty. In the example demonstrated in Fig. 10, 23% of the samples were excluded by the initial robust sky fit and the plotted histogram is of the surviving samples. The mode from this histogram is then added to the constant term of the 2-d polynomial fit and the width of the grey shaded region is taken to represent the sky noise in the image. In the histogram shown, the portion above the gaussian represents pixels that are contaminated by faint starlight. The fraction thus contaminated are too numerous to be excluded by typical outlier rejection algorithms. Without this correction to the sky fit, we would systematically fail to detect the faintest sources which happen to be the most likely to be of interest in our search. This approach still contains a small systematic bias to a higher sky level but this bias is significantly smaller than the sky noise level and was thus neglected. Note that zero in this plot is effectively the naive mean of the sky computed without sufficient  \nregard to the residual star contamination. In this case, the correction to the sky level is nearly 200 DN which is roughly a 3σ shift.', 'B.12. stacker.pro': 'Building image templates or stacks is done with this program. A set of images can be added up and a robust mean is estimated from the stack at each pixel. This tool also supports shifting the images to get a coherent stack on a object moving with respect to the pixel coordinates.', 'C. ICE HUNTERS - CITIZEN SCIENCE PARTICIPATION': "The rise of citizen science projects has been a consequence of attempts to harness the interest and attention of a large number of initially untrained people to tackle a large and either tedious or difficult to code task, generally in the form of sorting, classification, or searching through data. The search for KBOs in these crowded stellar fields was seen to be an interesting case.  \nFor this project, all the data were pre-processed as described in the main section of this paper. The basic product of these calculations were difference images where the stars are (mostly) subtracted and moving objects are seen as white and black pairs. We provided annotated training images to acquaint the participants with what to look for. After successful training, access to live data was granted.  \nThe images being scanned were much larger (2k by 4k) than could be displayed given the technology at the time. Rather than rely on image navigation driven by the user, we chose to present sub-images from the larger images. The size of the sub-array was dictated by the least capable graphical system to be used - in this case we were limited by the capabilities of the first version of the Apple iPad and a sub-array size of 424 x 424. To eliminate missing objects on the edge, we built-in a 10% overlap between adjacent sub-arrays. At the time of the sub-array creation, a fixed stretch of the image was used - scaled to -3 to +8 σ of the sky signal and the result was saved as a PNG-format image.  \nBreaking up the images resulted in a very large number of files. These files had to be loaded on a cloud server which was hosted on Amazon Web Services (AWS). The tool for uploading data to the server was GUI-based and designed for small batches of files. The transfer process was extremely time-consuming and would often fail if too many files were attempted at once. This problem proved to be a nearly fatal bottleneck for this project.  \nThe system was designed to present images to a participant for examination. Performing a complete search meant that each sub-image was viewed and processed by up to 16 different people with a small fraction seen as many as 40 times. The list of images was stepped through as each load request came in. Our design allowed for very fast processing of the data given enough participants. Once an image was viewed enough times, it was removed from the candidate list for processing. The goal during processing was to ensure that there were always images in the queue waiting for processing. In practice it was difficult to load images to AWS faster than they could be processed by the citizen scientists.  \nThis method of processing created a very large number of false positive detections. The total number of flagged transients was 728,954 and we had to apply some automated filtering steps to reduce the work load. Any transient seen ≥ 3 times was considered worthy of a confirmation effort. This threshold excluded 61% of the flagged transients. An additional 13% were excluded on the basis of visual examination of the candidate. The surviving candidates were visually graded into 3 different classes: a likely KBO (124, 0.02%), a likely main-belt asteroid (16,011, 2.2%), and a likely variable star (79,383, 11%). An additional 93,734 (13%) candidates could not be excluded or placed into any of these three classes and simply left marked as 'unknown' and not considered further.  \nThe final performance of this system versus our normal processing methods was the same for KBO detections. The citizen science approach was far more effective at finding and recording all of the other category of time variable elements contained in the data for essentially the same amount of effort. The citizen science approach would also be very effective if the data were populated with simulated objects for the purposes of characterizing the observational biases.  \nThe following is a list of all those that participated as an 'Ice Hunter:' T. Abbate, T. Abrahamson, M. Adam, J. Adamson, A. Agbedor, B. Alexander, U. Allan, X. Alldredge, A. M. H. Allen, D. Alton, M. P. Alves, V. S. Alves, O. Amerongen, W. L. Anderson, J. Andrews, S. Andrews, M. Anstett, K. Applebaum, J. Arrington, E. T. F. Ashton, W. Astl, L. G. Atkins, K. Badker, B. Badnaruk, E. Baeten, R. C. Bailey, R. Balick, J. L. Ball, B. E. Barnett, J. S. Bartel, C. Bartholomew, P. Baumgaertner, A. Becker, S. Becker, ' E. Bedard, I. Beer, J. Begg, J. Begonja, P. Bel, S. Berry, J. Bevc, V. A. Bezugly, G. Bianchini, M. Biel, B. Bigelow, M. Billiani, J. Bishop, M. C. Blanaru, S. Blystone,  \nG. Boden, T. Boesch, R. Bohnes, C. Boland, F. Bomeny, J. Borck, D. Borncamp, M. Boschat, G. P. A. Boscher, M. A. Boutilier, S. Boyd, N. B. Bradford, S. Bratcher, P. Brayshaw, A. Brinkmann, C. Brittain, D. Brogan, R. Brooks, D. Browne, T. Brydon, D. Burluraux, C. Butler, P. Caado, N. G. Cairns, F. Callebout, D. Cameron, P. Campbell, J. Campos, J. Candeias, A. Caracciolo, S. Carlino, A. Carpenter, C. Carrascal, J. Carruthers, U. Carsenty, F. Castro, L. D. Chadwick, E. Chaghafi, B. Chamberlain, K. Champion, C. Chandler, J. K. H. Cheng, D. Chestnov, C. Chinery, J. Chisholm, A. CIvinskas, A. Clark, D. Cleaver, M. Cline, S. Clingenpeel, B. Collins, E. Conseil, S. Corfield, R. Costello, R. Costin, M. Cotton, C. Craddock, M. Crosio, J. A. Cross, B. Crossland, C. Crowe, S. Cruze, R. Damerell-Moss, J. Danker, M. Darcy, G. David, J. M. A. Dawey, V. P. De A Barreto, C. De Grandis, D. de Vries, S. de Witt, K. Dekker, S. Delille, J. Demers, A. Dether, K. Devine, W. A. Dias, S. Diecks, E. Dimitrova, K. M. Doll, P. Dowden, P. N. Duggan, A. Dumitriu, L. M. Duncan, G. Dungworth, P. Dunlea, D. Dzierwa, A. Edwards, C. Elidoro, P. Ellinger, L. J. Elsnau, M. Elustondo, Y. A. Emery, H. F. Evans, J. Evans, S. Evdokimova, A. Eve, A. Fakahany, A. Farmer, P. K. Farmer, K. Farrelly, P. Fazekas, A. Felsenstein, N. Fequet, S. Finney, T. Fisher, M. S. Fletcher, R. Fletcher, K. Flippo, J. P. B. Fonseca, C. Foote, K. Fox, D. C. Fromm, S. Fronczak, R. Fuchs, H. Fukushi, M. Furskog, M. Gaebel, R. Gagliano, D. Gash, S. Ghegoiu, B. Gilbert, J. Gilbertson, R. Gill, A. Gillis, R. Glover, V. Gonano, S. Gonzaga, S. J. Goodman, A. Gorton, L. Goznell, H. Grady, T. Green, P. Greutmann, C. Griffing, P. Gronowski, R. Guerreiro, E. Guti'errez G'omez, D. Gutler, C. Guttau, P. Gwenn, L. Gyarmati, L. Hackl, T. Hager, E. Hall, J. Haller, T. Hantel, J. Harkin, A. Harlander, C. Hartel, E. Hartmeijer, R. F. Harvey, A. Hatwood, R. Heikes, H. Heilman, F. Helk, E. Helmersen, R. Hendricks, S. Hennessey, B. Hensley, B. Herentrey, M. Herr, J. Herridge, J. F. Higgins, L. Higgins, S. Hignell, T. Hodge, T. Hoffman, T. Hollerung, T. Horvath, D. Hrdlicka, M. Huefner, A. E. Huerta-Martin, P. Hughes, M.-T. Hui, C. Hunt, C. Hutchinson, N. Ihle, B. Imp, S. Ivanchenko, K. Izakovic, P. Jackson, W. Jammer, J. D. Jamulla, N. Janssen, A. Jaureguiberry, D. Jaworski, P. Jennings, R. Jewell, W. Johnson, L. F. Jones, S. L. Kaczmarek, A. Kaczmarowski, M. Kaluzny, N. Kaminski, M. Kamp, K. Kanegae, S. Karge, A. Karl, D. Karthikeyan, A. Kastengren, S. Kavalli, N. Kaylor, B. Kelly, C. Kendall, A. Kennedy, B. Kerr, S. Kidner, M. Kidulich, R. King, B. Kirschenmann, J. H. Klouwen, K. Klug, K. Koester, M. Kolling, R. Korber, P. Kostrzewa, R. Kracht, H. Krawczyk, J. W. Krawczyk, R. Krikken, P. Kr'ol, H. Kumar, M. Kumrucu-Lohmiller, P. Kyriakides, D. J. Lactin, A. Lamperti, J. Lander, M. Larson, F. Laurentiu, C. Le Garff, P. Le Gendre, K. Le Tourneau, P. Ledin, D. V. Leegwater, N. Lenke, P. A. Leon, A. V. Levin, D. Lindberg, J. Lipinski, C. Lloyd, V. P. Loeffler III, R. Lopez-Fabrega, E. Luchinat, S. Luers, J. Lummus, P. J. Lunn, A. Lysiak, C. Macmillan, A. Macumber, R. Madala, K. L. P. Madsen, S. Magee, R. Maherjinqiu, C. Mangili, A. Mankevich, B. Manning, D. Marion, A. Martin, D. Martinez, J. A. Martins, M. H. Massuda, B. Matter, C. Mays, I. Mazouzi, K. Mazowiecka, M. T. Mazzucato, R. McClure, B. McDaniel, P. McGarry, G. D. McKee, D. McMillan, R. Mellor, P. Mellors, J. Merc, S. Mercer, K. Miller, A. Mimeev, G. Mitchell, J. L. Moe, G. Mohr, B. B. Molders, R. M. Moldoveanu, E. Moody, C. Mooney, C. Moreno-Martinez, C. Morford, A. Morgan, D. Morgan, K. Morrill, V. Mottino, B. Moyant, W. Muehlisch, E. Muller, T. Myllymaki, S. Napier, J. Naujok, C. Neely, M. Negus, D. Nelson, S. Nembrini, J. Nethercott, B. Newgrosh, C. Newman, O. Nica, J. Nicholson, M. Nicholson, P. Nicholson, M. Nissinen, J. Nordnes, A. G. Norris, O. North, P. Northrop, K. Nuber, M.-C. Nuta, G. O'Callaghan, A. Oikkonen, R. Olejarski, R. Orsval, J. Ostler, K. Ozga, N. N. Paklin, A. Pandey, M. Pani, P. Pantoni, E. R. Pearsall, J. K. Pedersen, J. Peixoto, S. Pemberton, C. Perello, P. R. G. Perlmutter, E. Perozzi, B. Perry, D. Pettitt, M. Phillips, K. Pidgley, N. Piekarska, R. Pierce, F. Piroddu, T. Playle, M. Plewinska, L. Poffel, S. Pogrebenko, T. Pollard, C. W. Poole, A. Post, D. L. M. Preston, T. Pulokas, J. Purontakanen, E. Pusch, S. Pyne, J. Quinn, E. Quiroga-Rodr ' Iguez, H. Raab, C. A. Radford, V. M. Rahimi, C. E. Ravasi, I. Regan, F. Reisch, M. Renz, B. Replogle, J. Richmond, E. Rike, Q. Ringsak, F. Ringwald, A. E. M. Rojas, J. P. Ronsen, K. Roovers, J. Ros, J. Rossman, S. L. Roubroeks, T. Rounds, G. Roynon, C. Ruedinger, M. Ruh, N. Rutledge, G. F. Salinas, O. Sanislav, H. Sankala, L. Saracino, M. Savels, J. A. Sayle, T. Schafer, U. Scheuss, M. Schindler, S. Schmalz, K. Schmida, T. Scholze, P. Schulz, L. Schupper, J. Scott, M. Scott, J. Sejpka, G. Selig, B. M. Shaw, G. Shearwater, L. Sheldon, J. Shelton, A. Sheppard, A. Shirvanian, M. Shockley, K. L. Siddall, M. J. Sijens, P. Silva, D. J. Simister, M. Simmons, D. Skillman, R. Slivka, I. A. Smith, M. Smith, R. Smith, J. Snyder, B. Sobczuk, L. L. Somsikova, M. Souders, M. Spanjaard, A. Springer, T. Stallone, M. Stano, D. Stelter, R. W. Stephan, R. Stevens, P. D. Stewart, A. D. Stoca, J. Stuart, J. Sue, P. S. Summerlin, S. R. Taylor, M. Tepper, R. Tertnes, J. Thebarge, P. Thomaidis, J. Thornton, D. Tomic, J. Tompkins, T. Tormanen, A. Tors˚a, G. Triltzsch, D. Truog II, D. Tukendorf, T. Turaev, I. Uhl, C. Unger, C. van Boven, R. Van Der Hilst, P. Van Grijfland, M. Veillette, R. H. B. Velthuis, P. Verdelis, M. Verga, G. Verhoeven, A. Vertinskis, P. Vigil, D. Volkel, T. Vorley, T. Vorobjov, V. Vorochilov, L. L. W. Wah, E. Walravens, D. Ward, S. Ward, J. D. Warner, J. Wartenberg, G. Washbourn, P. Waterfield, S. Wegert, A. Wegner-Kaminski, D. Weilant, D. Westbrook, J. Wheat, N. Wheeldon, D.  \nWhitfield, A. Wille, M. V. Winkle, G. Wollenhaupt, D. Woodhams, A. Wrobel, C. Yandle, J. Zeddam, D. R. Zeigler, L. Zinke, B. Zinn, and I. Zormpa.  \nFacility: Subaru, Magellan, HST, CFHT, Gemini", 'REFERENCES': 'Monet, D. G., Levine, S. E., Canzian, B., et al. 2003, AJ, 125, 984, doi: 10.1086/345888 Parker, A. H., & Kavelaars, J. J. 2010, Icarus, 209, 766, doi: 10.1016/j.icarus.2010.04.018 Petit, J. M., Kavelaars, J. J., Gladman, B. J., et al. 2011, AJ, 142, 131, doi: 10.1088/0004-6256/142/4/131 Porter, S. B., Buie, M. W., Parker, A. H., et al. 2018, AJ, 156, 20, doi: 10.3847/1538-3881/aac2e1 Porter, S. B., Spencer, J. R., Verbiscer, A., et al. 2022, PSJ, 3, 23, doi: 10.3847/PSJ/ac3491 Smart, W. M., & Green, E. b. R. M. 1977, Textbook on Spherical Astronomy (Cambridge University Press) Spencer, J., Buie, M., Young, L., Guo, Y., & Stern, A. 2003, Earth Moon and Planets, 92, 483, doi: 10.1023/B:MOON.0000031963.58573.97  \nStern, S. A. 2008, Space Sci. Rev., 140, 3, doi: 10.1007/s11214-007-9295-y  \nStetson, P. B. 1987, PASP, 99, 191, doi: 10.1086/131977  \nVerbiscer, A. J., Porter, S., Benecchi, S. D., et al. 2019, AJ, 158, 123, doi: 10.3847/1538-3881/ab3211  \nVerbiscer, A. J., Helfenstein, P., Porter, S. B., et al. 2022, PSJ, 3, 95, doi: 10.3847/PSJ/ac63a6  \nZacharias, N., Finch, C. T., Girard, T. M., et al. 2012, VizieR Online Data Catalog, I/322A  \nZacharias, N., Finch, C., Subasavage, J., et al. 2015, VizieR Online Data Catalog, I/329', 'All Authors and Affiliations': 'Marc W. Buie , 1 John R. Spencer , 1 Simon B. Porter , 1 Susan D. Benecchi , 2 Alex H. Parker , 1 S. Alan Stern , 1 Michael Belton , 3 Richard P. Binzel , 4 David Borncamp , 5 Francesca DeMeo , 4 S. Fabbro , 6 Cesar Fuentes , 7 Hisanori Furusawa , 8 Tetsuharu Fuse , 9 Pamela L. Gay, 10 Stephen Gwyn, 11 Matthew J. Holman , 12 H. Karoji, 13 J. J. Kavelaars , 11 Daisuke Kinoshita, 14 Satoshi Miyazaki , 15 Matt Mountain, 5 Keith S. Noll , 16 David J. Osip , 17 Jean-Marc Petit , 18 Neill I. Reid, 5 Scott S. Sheppard , 19 Mark Showalter , 20 Andrew J. Steffl , 1 Ray E. Sterner, 21 Akito Tajitsu , 22 David J. Tholen , 23 David E. Trilling , 24 Harold A. Weaver , 21 Anne J. Verbiscer , 25 Lawrence H. Wasserman, 26 Takuji Yamashita, 27 Toshifumi Yanagisawa, 28 Fumi Yoshida , 29 and Amanda M. Zangari 1  \n1 Southwest Research Institute, 1301 Walnut St., Suite 400, Boulder, CO 80302  \n2 Planetary Science Institute, 1700 East Fort Lowell, Suite 106, Tucson, AZ 85719  \n- 3 Belton Space Exploration Initiatives, 430 S. Randolph Way, Tucson, AZ 85716; Emeritus Astronomer, Kitt Peak National Observatory,  \nTucson, AZ 85719  \n4 Department of Earth, Atmospheric, and Planetary Sciences, Massachusetts Institute of Technology, 77 Massachusetts Avenue,  \nCambridge, MA 02139  \n5 Space Telescope Science Institute, INS, Baltimore, MD, 21218  \n- 6 National Research Council of Canada, Herzberg Astronomy & Astrophysics Program, 5071 West Saanich Road, Victoria, BC, V9E 2E7, Canada  \n7 Northern Arizona University, Flagstaff, AZ  \n- 8 Astronomy Data Center, National Astronomical Observatory of Japan, Osawa 2-21-1, Mitaka, Tokyo 181-8588, Japan  \n9 Public Relations Center, National Astronomical Observatory of Japan, Osawa 2-21-1, Mitaka, Tokyo 181-8588, Japan  \n10  \nSouthern Illinois University  \n- 11 National Research Council of Canada, Herzberg Astronomy & Astrophysics Research Centre, 5071 West Saanich Road, Victoria, BC V9E 2E7, Canada  \n12 Center for Astrophysics, Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138  \n13 National Astronomical Observatory of Japan, Osawa 2-21-1, Mitaka, Tokyo 181-8588, Japan  \n14 Institute of Astronomy, National Central University, Jhongli, Taiwan  \n15 Subaru Telescope, National Astronomical Observatory of Japan, 650 N Aohoku Place Hilo Hawaii 96720 USA  \n16 NASA Goddard Spaceflight Center, 8800 Greenbelt Rd., Code 693.0, 20771, Greenbelt, MD  \n17 Carnegie Observatories, Las Campanas Observatory (Chile)  \n18 Institut UTINAM UMR6213, CNRS, OSU Theta F-25000 Besancon, France  \n19 Earth and Planets Laboratory, Carnegie Institution for Science, 5241 Broad Branch Rd. NW, Washington, DC 20015  \n20 SETI Institute, 189 Bernardo Ave., Mountain View, CA 94043  \n21 Johns Hopkins University Applied Physics Laboratory, 11100 Johns Hopkins Road, Laurel, MD 20723  \n- 22 Subaru Telescope Okayama Branch, National Astronomical Observatory of Japan, 3037-5, Honjou, Kamogata, Asakuchi, Okayama 719-0232, Japan  \n23 Institute for Astronomy, 2680 Woodlawn Drive, Honolulu, HI 96822  \n24 Department of Astronomy and Planetary Science, Northern Arizona University, P.O. Box 6010 Flagstaff, AZ, 86011; School of Informatics, Computing, and Cyber Systems, Northern Arizona University, P.O. Box 5693, Flagstaff, AZ, 86011  \n25 Department of Astronomy, University of Virginia, Charlottesville, VA  \n26 Lowell Observatory, 1400 W Mars Hill Road, Flagstaff, AZ 86001  \n27 Thirty Meter Telescope (TMT) Project, National Astronomical Observatory of Japan, Osawa 2-21-1, Mitaka, Tokyo 181-8588, Japan 28 Chofu headquarters, Japan Aerospace Exploration Agency, Jindaiji Higashimachi 7-44-1, Chofu, Tokyo 182-0012, Japan 29 School of Medicine, Department of Basic Sciences, University of Occupational and Environmental Health, 1-1 Iseigaoka, Yahata, Kitakyusyu, Fukuoka 807-8555, Japan; Planetary Exploration Research Center, Chiba Institute of Technology, 2-17-1 Tsudanuma, Narashino, Chiba 275-0016, Japan'}