Issue 
A&A
Volume 615, July 2018



Article Number  A159  
Number of page(s)  8  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/201732254  
Published online  01 August 2018 
A new approach to distant solar system object detection in large survey data sets
Hamburger Sternwarte, Universität Hamburg, Gojenbergsweg 112, 21029 Hamburg, Germany
email: vperdelwitz@hs.unihamburg.de
Received:
7
November
2017
Accepted:
19
April
2018
Context. The recently postulated existence of a giant ninth planet in our solar system has sparked search efforts for distant solar system objects (SSOs) both via new observations and archival data analysis. Due to the likely faintness of the object in the optical and infrared regime, it has so far eluded detection.
Aims. We set out to reanalyze data acquired by the WideField Infrared Survey Explorer (WISE), an allsky survey well suited for the detection of SSOs.
Methods. We present a new approach to SSO detection via parallactic fitting. Using the heliocentric distance as a fit parameter, our code transforms groups of three or more singleobservation point sources to heliocentric coordinates under the assumption that all data stem from an object. The fact that the orbit of a distant SSO is approximately linear in heliocentric coordinates over long timescales can be utilized to produce candidates, which can then be confirmed with followup observations.
Results. We demonstrate the feasibility of the approach by a posteriori detecting the outer SSO Makemake within WISE data. An allsky search for Planet Nine yielded no detection.
Conclusions. While the postulated Planet Nine eluded detection by our algorithm, we tentatively predict that this new approach to movingobject analysis will enable the discovery of new distant SSOs that cannot be discovered by other algorithms. Especially in cases of sparse data observed over long time spans, our approach is unique and robust due to the use of only one fit parameter.
Key words: planets and satellites: detection / surveys / minor planets, asteroids: general / minor planets, asteroids: individual: Makemake
© ESO 2018
1. Introduction
The last few decades have seen great advances in the understanding of the outer solar system. Seventyfour years after the discovery of Pluto (Tombaugh 1946), Brown et al. (2004) report the discovery of Sedna, the first of a number of objects of the inner Oort cloud to be detected in subsequent years. The discovery of 2012 VP_{113} (Trujillo & Sheppard 2014) leads the authors to conclude that a perturber at 250 au could explain the apparent clustering in the arguments of perihelion of the distant scattered disk population around ω ≈ 0^{°}, an effect also observed in subsequently discovered objects (Sheppard & Trujillo 2016; Bannister et al. 2017). The underlying simulations were, however, limited, and the authors stress that the properties of the potential perturber are not unique. Batygin & Brown (2016a) not only show that this clustering extends to physical space (i.e., orbital planes), but that it can be explained by the existence of a ninth planet (referred to here as Planet Nine) in our solar system with a mass of ≈10 M_{⊕} on an orbit with a ≈700 au semimajor axis. This publication has sparked a vivid debate on the existence of Planet Nine, leading to publications which support the hypothesis (e.g., Bailey et al. 2016; Lai 2016; Batygin & Brown 2016b; Brown 2017; Gomes et al. 2017; Becker et al. 2017; de la Fuente Marcos & de la Fuente Marcos 2017; Batygin & Morbidelli 2017), as well as those that argue against it (e.g., Lawler et al. 2017; Nesvorný et al. 2017; Shankman et al. 2017). Furthermore, other authors derive low probabilities for the production of the proposed orbit (Li & Adams 2016). The fact that Planet Nine has not been detected to date must then be due to one of the following premises: (i) It does not exist; (ii) previous surveys lack the necessary sensitivity; or (iii) the object’s signature is present in survey data but has hitherto eluded all applied detection algorithms.
Several authors provide restrictions on the location of Planet Nine based on Cassini data (Fienga et al. 2016; Holman & Payne 2016b), motion of comets (Medvedev et al. 2017), a Monte Carlo approach (de la Fuente Marcos & de la Fuente Marcos 2016), astrometry of transNeptunian objects (TNOs; Holman & Payne 2016a), mean motion resonance (Malhotra et al. 2016; Millholland & Laughlin 2017) and the sky coverage of sufficiently deep surveys (Brown & Batygin 2016). Despite these spatial constraints, the search for an object as distant as the proposed Planet Nine illustrates the tradeoff between sky coverage and limiting magnitude. Brightness estimates (Fortney et al. 2016; Toth 2016; Linder & Mordasini 2016) place it at the edge of the detection limit of most large surveys, while a dedicated search with large telescopes is relatively expensive in terms of observing time.
However, Planet Nine may be hidden within existing survey data. A variety of methods have been proposed to discover faint, moving objects within these vast data sets. For example, Kuchner et al. (2017) have developed a citizen science project in which volunteers visually check images acquired by the WideField Infrared Survey Explorer (WISE, Wright et al. 2010), and Meisner et al. (2018) use timeresolved coadds to rule out the existence of Planet Nine within the WISE data up to a W1 magnitude of 16.7. Despite differing in the precise implementation, most other methods for the detection of distant TNOs (e.g., Petit et al. 2004; Denneau et al. 2013; Brown et al. 2015; Weryk et al. 2016) are based on the same principle: data acquired during a time span of several days is utilized to identify tracklets, that is moving sources. Tracklets from different epochs are then crosschecked via an algorithm based on or similar to the work of Bernstein & Khushalani (2000). All of these algorithms seem to require multiple tracklets within a time frame of ≈60 days. In a recent publication, Holman et al. (2018) describe the discovery of a TNO with the PanSTARRS1 Outer Solar System Pipeline (Holman et al. 2011), which transforms single detections from topo to heliocentric coordinates and subsequently searches for linearity. However, in this approach the heliocentric distance is not a fit parameter, but rather looped through stepwise.
While having successfully detected TNOs, all previous approaches based on tracklets rely on the existence of several exposures acquired within a short span of time. In order to query data in areas with low coverage depth and/or contamination through various effects, it is possible to implement an algorithm which corrects for the objects’ parallactic motion and which is then able to correlate sparsely sampled data acquired over long time intervals.
In this publication, we present our new approach to moving object detection, which we then employ to reexamine the WISE/NEOWISE singleexposure source catalogs. Throughout this publication, we use the TNO (136472) Makemake (Brown et al. 2005) to illustrate the challenge of movingobject detection and as a proofofconcept of our method. Our paper is structured in the following way: In Sect. 2, we describe the concept of the algorithm along with a mathematical description of the steps involved. In order to demonstrate functionality, we test the algorithm on the field containing Makemake in Sect. 3, followed by an allsky search for Planet Nine and injected artificial test planets. We conclude with a summary and outlook in Sect. 4.
2. A new approach to movingobject detection
Typical photometric catalogs contain millions or even billions of entries (e.g., the NEOWISE singleexposure source table alone comprises a total of ≈6 × 10^{10} entries). The goal of our new movingobject detection approach is to correlate those singular pointsource detections belonging to one physical object, without any boundary conditions on scan frequency and temporal spacing, using a specifically optimized algorithm in order to detect distant moving objects in such large data sets.
Our algorithm combines routines that filter, correlate, and select movingobject candidates with capabilities to handle solar system object (SSO) trajectories that are contaminated by random background stars. The approach can be divided into four steps, of which we give a summary here and which we explain in detail in Sects. 2.1–2.4.
Data selection and filtering: Choose the field of interest and remove all data with temporal duplicity (i.e., background stars).
Clustering: Form groups of singleobservation point sources which could originate in real moving objects.
Selecting: Compute all possible permutations of three point sources out of the clustered groups and iteratively try to linearize these subgroups via transformation to heliocentric coordinates; assume that in the time span of several years the orbit is linear in these coordinates; compile a set of subgroups with linearity in heliocentric coordinates at a corresponding heliocentric distance.
Backtracing and position prediction: All resulting subgroups are reclustered and selected, this time requiring a larger subgroup size. The linearized orbit of sufficiently large subgroups can be extrapolated to future dates so as to provide a position prediction for followup observations.
In the following we give a detailed description of the different steps performed by our algorithm (Fig. 1).
Fig. 1.
Flowchart of the search algorithm. All process are displayed as square boxes, while all (intermediate) catalogs are marked by rounded boxes. The corresponding section in the text is given in the right column. 
2.1. Data selection and filtering
Due to its spectral coverage and observing strategy, WISE (Wright et al. 2010) is wellsuited for the search for distant SSOs. Covering each area of the sky multiple times during halfyear intervals, the WISE mission provides allsky data spanning a total of 4 yr. We therefore used WISE data for the first tests of our algorithm, but we stress that the method works with all large area/large temporal coverage survey data sets. We downloaded singleexposure source lists from the WISE catalog via the NASA/IPAC Infrared Science Archive, requiring a separation from the South Atlantic Anomaly of more than 1^{°}, no direct neighbors in the 2MASS catalog (Skrutskie et al. 2006; Cutri et al. 2013), a detection in the W2 band (w2mpro ≠ null), and all other quality and saturation flags with appropriate values. Furthermore, the SSO Flag (sso_flg) in the WISE catalogs was required to be zero to avoid rediscovery of known objects.
In order to remove stationary background sources (e.g., stars and galaxies), we first checked for temporal duplicity and then filtered the data in the following manner: For each observed event we found all other data within a specific search radius. This radius is determined by two factors: (i) the maximum orbital motion of an SSO at a given distance and (ii) its maximum parallactic motion. In order to approximate these two contributions, we followed the reasoning of Cowan et al. (2016) and determined the hourly orbital movement assuming a semimajor axis of 700 au as(1)
and the hourly parallactic movement as(2)
where d denotes the distance of the object. Throughout most of the sky, WISE covers a given region ≈12 times per half year, and all visits occur within ≈1 day, so the overall apparent movement becomes μ ≈ 10 arcsec d^{−1}, assuming a distance of d = 180 au. We note that, since singular points (i.e., those without neighbors within the search radius) can be carried along for further processing (as explained below), the choice of a smaller radius does not result in the loss of true events. It does, however, yield a higher number of points for further processing. The choice of a larger search radius, on the other hand, eliminates real SSO observations along with the background
For those groups with two or more points, the average right ascension, declination, and modified Julian date (MJD) was calculated. Finally, the angular velocity of each group was computed and used as a further filtering tool. Since the angular velocity caused by Earth’s orbit is dependent on distance of the object, time, and direction, we used the AstroPy package to compute an upper and lower limit for the angular motion within one day by defining a minimum and maximum distance search window. Allowing for a margin of 1 arcsec d^{−1} to account for the orbital motion of objects, we selected only those objects exhibiting a velocity in the determined window for further processing.
2.2. Clustering
Using the compiled input catalog, the algorithm computes clusters of data which could physically belong to an orbit. This is done by comparing every pair of object vectors(3)
in terms of their temporal distance Δt and angular great circle distance(4)
A search radius is defined by assuming a proper motion and a parallax, inside of which a pair of vectors X_{i} is accepted to belong to a cluster. The number of vectors X_{i} clustered together is then called the cluster size.
For a nearly geocentric observer, all TNOs show a characteristic loop whose shape is controlled by the Earth’s orbit as well as the proper motion of the object and its coordinates. In order to illustrate this we used the HORIZONS WebInterface^{1} to compute an artificial sample of observations of Makemake in the period 2016 Jan 1 and 2017 Sep 1 as displayed in Fig. 3.
Fig. 2.
Schematic of the vectors introduced in Sect. 2.3.1. The Sun is denoted by ⊙, Earth by ⊕ and the TNO by ◯. 
Fig. 3.
Makemake geocentric path between 2010 May 13 and 2015 April 7 generated using the HORIZONS WebInterface (blue line). All WISE observations during which Makemake was in the field of view are marked with a dot. The WISE detections of Makemake are marked with a red cross. 
An upper limit for the semimajor axis of the parallactic loop is given by(5)
where d_{min} denotes the minimal geocentric distance. For a given temporal separation Δt, we can estimate an upper limit for the parallactic component of the search radius via(6)
where ω_{⊕} = 2π/yr is the Earth’s mean orbital angular velocity. For distant SSOs, the parallactic loop is mildly distorted by a proper motion component. Imposing a minimum semimajor axis and a minimum distance d_{min}, we can calculate a maximum perihelion orbital velocity(7)
where we approximate the object’s mass to be zero. With this peak velocity, we can calculate a maximum daily proper motion component via(8)
Finally, the total search radius is now determined via(9)
The clustering algorithm was implemented in the FORTRAN programming language and fully parallelized for sharedmemory machines using OpenMP. Even large sets of data containing millions of objects can be clustered within hours or a few days. An example of a typical cluster size distribution is given in Sect. 3.2. For this example, the required computing time on a single Intel^{Ⓡ} Xeon^{Ⓡ} CPU E52680 v3 @ 2.50 GHz was ≈20 h.
2.3. Selecting
Even for nonSSOs, we have a nonzero probability of two objects having just the right temporal and angular separation to be regarded as part of the path of one virtual new SSO. Thus, we established a filtering pipeline that selects the physically plausible clusters. To do so, all clusters are transformed from their initially nearly geocentric reference system into a heliocentric system assuming a distance d. As a result, paths of actual distant SSOs are to a good approximation linear in the (α_{⊙}, δ_{⊙}) plane over long timescales if d is chosen correctly. Figure 4 shows the same observations as Fig. 3, but transformed into heliocentric coordinates. Therefore, the core idea of our algorithm is to assume that a given cluster represents the observed path of an SSO within some search distance range [d_{min}, d_{max}]. Subsequently, the algorithm tries to linearize the cluster by applying subsequent transformations to the heliocentric system and finding the bestfit distance.
Fig. 4.
Illustration of the Makemake orbit in heliocentric coordinates. Transformation to heliocentric coordinates has removed the ellipsoidal component of the orbit and the epochal sets of observations (chronologically from upper right to lower left) are linearized to good approximation. All WISE observations during which Makemake was in the field of view are marked with a dot. The WISE detections of Makemake are marked with a red cross. 
2.3.1. Linearization
For the sake of simplicity and to reduce the computational effort, we apply two simplifications. First, we perform no barycentric corrections leading to peak errors on the order of ≤0.01 au. Furthermore, we assume that the SSO’s heliocentric distance does not change significantly throughout the measurements. Depending on the heliocentric distance and the time range of the cluster, additional uncertainties up to ≈0.2 au can occur. Both simplifications are easily met by distant Kuiper belt objects and lead to total positional errors of just a few percent (assuming heliocentric ranges greater than 30 au).
We consider a heliocentric coordinate system and the xaxis pointing towards the vernal equinox.
We take SE to be the vector pointing from the Sun towards Earth, EP the vector pointing from Earth towards the SSO candidate and SP = (x_{p}, y_{p}, z_{p})^{T} the vector pointing from the Sun towards the SSO candidate (see Fig. 2 for a visualization). We make use of the IAU SOFA library^{2} for the calculation of the SunEarth position vector. Further, we take α and δ to be the geocentric equatorial coordinates of the SSO candidate and γ to be the angle between SE and EP. The EarthSSO vector can be written as EP = Δ · e_{EP}, denoting the EarthSSO distance as Δ and e_{EP} = (cos δ cos α, cos δ cos α, sin α)^{T} being the unit vector of EP. Finally, we take d to be the distance between the Sun and the SSO. From this, we calculate γ as(10)
with SunEarth unit vector e_{SE} and the negative sign to make sure that we calculate the inner angle between the two vectors.
Now, according to the law of cosines, the identity(11)
holds for the triangle spanned by the three vectors, which is equivalent to(12)
with the SunEarth distance r_{⊕}. Solutions for Δ are given by(13)
but as we are interested in objects beyond Earth’s orbit, only the positive solution is relevant. This equation implies that for a given guess of the heliocentric distance d, the geocentric distance can be calculated, and thus the full SunSSO vector via SP = SE + Δe_{EP}. Given that vector, we can calculate heliocentric equatorial coordinates α_{⊙} and δ_{⊙} using(14)
In these coordinates, any set of vectors of equatorial positions for some given time t_{i}(15)
can be transformed into a set of vectors of heliocentric positions(16)
for a given heliocentric distance guess d. Random sets of vectors X_{i} are not expected to show any special behavior in the (α_{⊙}, δ_{⊙}) plane, while those sets of vectors X_{i} belonging to an SSO are transformed into a line in the (α_{⊙}, δ_{⊙}) plane if the provided distance guess d proves to be consistent with the object’s true heliocentric distance. The quality of the bestmatch heliocentric range is estimated by calculating the linear Pearson correlation coefficient for the set of n vectors X_{⊙,i} through(17)
Actual SSOs with a large semimajor axis are transformed into lines in the (α_{⊙}, δ_{⊙}) plane with correlation coefficients close to unity. Large samples of data covering sufficiently large temporal intervals lead to welldefined maxima of the squared Pearson correlation coefficient which can be calculated via numerical maximization algorithms. Figure 5 shows as a function of heliocentric distance for a sample of Makemake data sets. Notably, the bestfit heliocentric distance coincides with the mean heliocentric distance of the object for the given time span.
Fig. 5.
Squared Pearson correlation coefficient ρ^{2} as a function of the heliocentric distance guess d for a set of 37 calculated geocentric positions of Makemake between 01012017 and 01012018. The function peaks prominently at 52.5 au which is the average heliocentric distance of Makemake for the given time interval. 
Very small samples (≈3) covering small temporal intervals (less than a month) feature less prominent ρ^{2} peaks and show a strong dependence on the observation time. For instance, during quadrature the apparent motion takes a minimum. As a result, we obtain smaller geocentric arcs, less well defined ρ^{2} peaks and more complicated ρ^{2} functions in general, which may result in incorrect bestfit distances if the provided fitting range is too large (see Fig. 6). Therefore, we apply a number of postfit sanity checks to detect failed linearizations (see Sect. 2.3.2).
Fig. 6.
Squared Pearson correlation coefficients as functions of the heliocentric distance guess for different samples (1day cadence) of Makemake ephemeris, all covering 1 month in 2017. Despite almost identical time spans, numbers of observations and observation frequency, the ρ^{2} functions differ greatly. For such data, the fitting range [d_{min}, d_{max}] must closely bracket the actual physical peak. Otherwise, the maximization routine may end up stuck at the borders of the fitting range or at one of the secondary peaks. While this could be avoided by checking linearity for small steps in heliocentric distance, this would void the speed of the fitting routine; so instead, the analysis is carried out with steps of d_{max} − d_{min} = 10 au. 
We conclude that our linearization algorithm works best for larger group sizes with a minimum of three data points covering a sufficiently long temporal range of ten days or more. Larger time spans and subgroup sizes increase the quality of the fit. For very small time spans, higher sampling rates are necessary in order to reproduce at least a small fraction of the parallactic loop and obtain the correct distance. However, tools such as MOPS cover that parameter range and are a natural complement to our approach. Compared to full orbital fitting routines which typically minimize parameter vectors with six fit parameters (see e.g., Bernstein & Khushalani 2000), our simplified approach deals with only one single fit parameter (i.e., the heliocentric distance d) to identify SSO candidates, and allows us to process much larger sets of data.
2.3.2. Subcluster analysis
Given the nonzero probability of cluster contaminations through nearby stars, we perform a subcluster decomposition for every SSO candidate cluster identified by our pipeline. Explicitly, for every cluster of size n we consider every possible permutation of subsets with subgroup size k leading to a total number of subclusters per cluster of(18)
ensuring a piecewise direct selection of uncontaminated planetary trails. For all of our first runs, we choose a subgroup size of k = 3 or k = 4. All subclusters are transformed into the bestfit heliocentric distance and tested for linearity as described in Sect. 2.3.1. Because of the large number of subclusters per cluster (e.g., 1140 for n = 20 and k = 3), we apply additional constraints to reduce the number of false positives. First, we require a proper time sequence of the heliocentric coordinates. Secondly, we check if the subgroup path length is consistent with the proper motion of an SSO. Finally, recalling that actual planetary paths show very high linearity in heliocentric coordinates, we require a squared Pearson correlation coefficient of . In Fig. 7 we applied a subcluster decomposition to the Makemake data cluster shown in Fig. 5 with a subgroup size of k = 4, resulting in subclusters spanning between 40 and 360 days, and calculated the distribution of the squared Pearson correlation coefficients. The distribution peaks at values larger than 0.9999.
Fig. 7.
Distribution of the squared Pearson correlation coefficients of all subgroups of size k = 4 for the Makemake data cluster presented in Fig. 5. The vast majority of subgroups are wellfitted and yield good linearity values larger than 0.9999. 
2.4. Backtracing and position prediction
After the subcluster analysis, we are typically (for the whole WISE catalog) left with ≈100 subgroups deg^{2}, every one of which consists of the coordinates and time of observation as well as the bestfit distance. As the probability of producing subgroups of size three with a high degree of linearity from random data (noise) is nonzero, all data from the previous step (Sect. 2.3.2) are then reclustered and selected, this time requiring the same squared Pearson correlation, but with a larger subgroup size. Since every group of four must consist of groups of three, all valid candidates are preserved while eliminating random matches. This process can be repeated until the required subgroup size deemed a “real” candidate is reached. In order to be able to carry out followup observations, the linearized (i.e., heliocentric) orbit is then extrapolated to the observation date and transformed to geocentric coordinates.
3. Results
3.1. Makemake
As a proofofconcept we tested our algorithm by performing a blind detection of Makemake (Brown et al. 2005), that is with no prior assumptions. The Known Solar System Association Table of the NEOWISE mission denotes a total of 32 detections of Makemake, which were attributed to it a posteriori. However, the SSO was in the fieldofview a total of 124 times, meaning that it was not detected in the majority of exposures. This can be explained by the fact that, with a mean brightness in bands W1 and W2 of and , Makemake is at the detection limit of WISE single exposures. Furthermore, it exhibits a large variability in brightness in both bands, with a range of 3 mag in W1 and 2 mag in W2, caused by effects such as rotational modulation, changes in background and noise. Figure 8 shows the light curve of Makemake in the W1 and W2 passbands for the entire WISE data set released so far, along with nondetections. It therefore presents a good example to illustrate the difficulty of discovering such objects within catalog data via standard methods, such as the search for “tracklets”, as well as a test for the method described in this publication.
Fig. 8.
Light curve of Makemake in the WISE bands W1 (red crosses) and W2 (blue dots). The upper row shows the entire light curve, while the second row shows a zoom of the observation epochs in 1.5 d windows. Nondetections are marked with an arrow at the S/N = 5 level. We note that most detections have a S/N below 5. 
We downloaded WISE data covering the Makemake orbit during the period 2010–2016 (i.e., α = 187.5 to 194, δ = +25 to +29) with the settings described in Sect. 2.1, producing a total of 280 515 point sources^{3}. Since the sso_flg was set to zero in the standard data acquisition, all WISE detections of Makemake with a detection and error estimate in both bands (a total of 16 individual detections) were added manually at this stage, that is, before any data reduction, thereby ensuring a representative trial. The selection and velocity filtering yielded a total of 4472 groups consisting of two single detections, and 918 groups consisting of more than two singular detections. Therefore, after computing the mean coordinates and date, we were left with 5390 points in the input catalog. After applying the velocity filter described in Sect. 2.1, this number was reduced to 4836, which constituted our input catalog for the clustering and selection processes.
Running the processes described in Sects. 2.2 and 2.3 with a subgroup size of three, we used the selected subgroups as input for a second run with subgroup size 4, which left a singular candidate with a Pearson coefficient of ρ^{2} = 0.999999999 at a bestfit heliocentric distance of 52.6 au. A lookup in the original data and comparison to the Known Solar System Object Possible Association List shows that the subgroup constitutes seven singular detections of Makemake. As a traceback reveals, the remaining detections were discarded during the data selection process due to proximity to background sources. However, transforming the entire input catalog to heliocentric coordinates assuming the derived bestfit distance and identifying all sources which fit the linearity in α_{⊙} and δ_{⊙} as a function of time to within 2^{″}, we retrieved all detections of Makemake from the data set. We therefore conclude that our approach allows for the discovery of Makemake based on WISE data alone.
3.2. Allsky search for Planet Nine and blind test
As described in Sect. 1, there are several, sometimes conflicting, position predictions for Planet Nine. For this reason, we did not limit our search to a specific region of the sky, but rather perform an allsky analysis. We obtained WISE data of this region with the flag settings described in Sect. 2.1 and without any a priori assumptions regarding colors or brightness. In this case, however, we limited our input data to the NEOWISE mission phase, since the presence of the AllWISE flag within this catalog allowed us to discard a large fraction of the data from the beginning. The query yielded a total of ≈4.2 × 10^{7} individual point sources. In Fig. 9 we illustrate the given source density per square degree. Using the same settings as described in Sect. 2.1, the removal of stationary sources resulted in ≈ 25 × 10^{6} data points.
Fig. 9.
Density of detected sources (sources/deg^{2}) of the input catalog. For clarity, the color code has been truncated at densities above 2000, since the regions at the ecliptic poles exhibit source densities almost two orders of magnitude above the average. The injected simulated Planet Nine observations are overplotted as black dots. 
To ensure that our algorithm is capable of finding SSOs at distances proposed for Planet Nine, we carried out blind tests very similar to that described in Sect. 3.1. The basic idea was to generate artificial Planet Nine data points which in principle could have been observed by WISE. To this end, we randomly generated 10^{3} sets of Keplerian orbital elements which included values for the semimajor axis a, the numerical eccentricity e, the orbital inclination i, the longitude of the ascending node Ω, the argument of the perihelion ω, and the mean anomaly M. The limits for these values were chosen according to the parameters predicted by Brown & Batygin (2016), and the random values of the mentioned parameters are uniformly distributed within these limits. To determine the ephemerides of these artificial Planet Ninelike objects we used the publically available software tool PyEphem^{4}. We evaluated the ephemerides for alltime stamps available in the NEOWISE survey frame metadata table, resulting in 7.7 × 10^{6} data points for each simulated object, and crossmatched them using the framecorner coordinates, thus excluding all data not covered by simultaneous WISE pointings, yielding approximately 290 data points per object as observed by WISE. From these, we decided to randomly remove further data points to account for observations at the detection limit, blends, or spikes caused by bright background stars, as well as cosmic ray contamination. The number of deleted data points was randomly determined with the goal to achieve artificial object data sets ranging almost uniformly from 3 to 20 data points. Our artificial data set now constitutes 11 315 data points of 1000 Planet Ninelike objects in total. In Fig. 9 we show these data points in comparison to the density of data from the WISE catalog.
Finally, we injected the candidates into the real NEOWISE data. Based on the experience with the recovery of Makemake described in Sect. 3.1, we used a minimum required Pearson coefficient of ρ^{2} > 0.999999 (the minimum linearity of all subgroups of Makemake) and a minimum heliocentric distance of 180 au, resulting in a total of 4.7 × 10^{5} subgroups. In order to demonstrate the influence of pointsource density on the total computing time, we show the number of clusters as a function of cluster size as a solid line in Fig. 10. The total number of linearizations, and therefore computing time, scales with the product of the binomial coefficient and the number of clusters of a given size (dotted line). We note that the number of iterations converges toward a constant value, indicating that the processing time is manageable.
Fig. 10.
Number of clusters (solid line) and number of linearizations necessary to process the data (dotted line). The number of iterations converges toward ~10^{6} per cluster, an achievable value with reasonable computing power. 
Tests have shown that the number of subgroups of size 3 and 4 formed from random data (i.e., noise) is too large to regard each of them as candidates, so a successive analysis with subgroup sizes of 4 and 5 was carried out. We stress that it is also possible to start the analysis with a higher subgroup size, although this does require larger computing times. The resulting subgroups of size 5 were merged based on the occurrence of data within multiple groups, and our algorithm retrieved a total of 836 candidates of size 5 or more from our artificial Planet Nine candidates. Considering that 872 out of the 1000 artificial planets consisted of five or more points, our algorithm successfully retrieved 95.8% of the candidates, including many with only one observation per epoch. We could not find any systematics within the orbital parameters of the injected planets not recovered by the code.
Finally, since no candidates of size 5 or more were detected within the real NEOWISE data, and we were able to recover most of the injected test planets, we exclude the presence of five or more data points of Planet Nine within the NEOWISE catalog with a 2σ confidence level.
4. Summary and outlook
The detection of TNOs is a notoriously challenging task and requires great observational efforts. The visual discovery of such objects via systematic observations of welldefined parts of the sky has already been performed by many other groups. We present the development of a new, highly efficient algorithm, which can deal with the unprecedented amounts of archival data that have become available within the last years. With the additional capability of handling data covering large time spans, our code complements other existing movingobject detection algorithms.
We used WISE/NEOWISE data to test our algorithm by serendipitously detecting the known TNO Makemake. Our code proved to be fully capable of locating objects such as Makemake in crowded fields and is accurate enough to provide position predictions for followup observations.
Furthermore, we demonstrated that our algorithm would have been able to detect Planet Nine in WISE data, by injecting test planets into our data, of which we were able to detect 95.8%.
Our search for Planet Nine within WISE data was unsuccessful. It is, however, more likely that the object, should it exist, is not bright enough to be detected with the WISE satellite, which is in agreement with the nondetection of Planet Nine by Meisner et al. (2018).
Having demonstrated the capabilities of the approach, we continue the search for Planet Nine by analyzing a combined data set spanning decades and comprising all large public catalogs covering the visual and infrared range. Furthermore, the fact that the code is not reliant on a tight temporal spacing of data may allow for the search for Planet Nine via stellar occultations (however rare) by compiling candidate occultation catalogs from existing photometric surveys or dedicated projects (Lehner et al. 2012; Pass et al. 2018).
Acknowledgments
This publication makes use of data products from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation. This publication makes use of data products from the Widefield Infrared Survey Explorer, which is a joint project of the University of California, Los Angeles, and the Jet Propulsion Laboratory/California Institute of Technology, funded by the National Aeronautics and Space Administration. This research has made use of the NASA/IPAC Infrared Science Archive, which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration. The authors thank J.H.M.M. Schmitt for suggestions and proofreading.
References
 Bailey, E., Batygin, K., & Brown, M. E. 2016, AJ, 152, 126 [NASA ADS] [CrossRef] [Google Scholar]
 Bannister, M. T., Shankman, C., Volk, K., et al. 2017, AJ, 153, 262 [NASA ADS] [CrossRef] [Google Scholar]
 Batygin, K., & Brown, M. E. 2016a, AJ, 151, 22 [NASA ADS] [CrossRef] [Google Scholar]
 Batygin, K., & Brown, M. E. 2016b, ApJ, 833, L3 [NASA ADS] [CrossRef] [Google Scholar]
 Batygin, K., & Morbidelli, A. 2017, AJ, 154, 229 [NASA ADS] [CrossRef] [Google Scholar]
 Becker, J. C., Adams, F. C., Khain, T., Hamilton, S. J., & Gerdes, D. 2017, AJ, 154, 61 [Google Scholar]
 Bernstein, G., & Khushalani, B. 2000, AJ, 120, 3323 [NASA ADS] [CrossRef] [Google Scholar]
 Brown, M. E. 2017, AJ, 154, 65 [NASA ADS] [CrossRef] [Google Scholar]
 Brown, M. E., & Batygin, K. 2016, ApJ, 824, L23 [NASA ADS] [CrossRef] [Google Scholar]
 Brown, M. E., Trujillo, C., & Rabinowitz, D. 2004, ApJ, 617, 645 [NASA ADS] [CrossRef] [Google Scholar]
 Brown, M. E., Trujillo, C. A., & Rabinowitz, D. 2005, IAU Circ., 8577 [Google Scholar]
 Brown, M. E., Bannister, M. T., Schmidt, B. P., et al. 2015, AJ, 149, 69 [NASA ADS] [CrossRef] [Google Scholar]
 Cowan, N. B., Holder, G., & Kaib, N. A. 2016, ApJ, 822, L2 [NASA ADS] [CrossRef] [Google Scholar]
 Cutri, R.M., et al. 2013, VizieR Online Data Catalog, II/328 [Google Scholar]
 de la Fuente Marcos, C., & de la Fuente Marcos, R. 2016, MNRAS, 459, 66 [NASA ADS] [CrossRef] [Google Scholar]
 de la Fuente Marcos, C., & de la Fuente Marcos, R. 2017, MNRAS, 471, 61 [NASA ADS] [CrossRef] [Google Scholar]
 Denneau, L., Jedicke, R., Grav, T., et al. 2013, PASP, 125, 357 [NASA ADS] [CrossRef] [Google Scholar]
 Fienga, A., Laskar, J., Manche, H., & Gastineau, M. 2016, A&A, 587, L8 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fortney, J. J., Marley, M. S., Laughlin, G., et al. 2016, ApJ, 824, L25 [NASA ADS] [CrossRef] [Google Scholar]
 Gomes, R., Deienno, R., & Morbidelli, A. 2017, AJ, 153, 27 [NASA ADS] [CrossRef] [Google Scholar]
 Holman, M. J., & Payne, M. J. 2016a, AJ, 152, 80 [NASA ADS] [CrossRef] [Google Scholar]
 Holman, M. J., & Payne„ M. J. 2016b, AJ, 152, 94 [NASA ADS] [CrossRef] [Google Scholar]
 Holman, M. J., Protopapas, P., & Chen, Y., et al. 2011, in AAS Meeting Abstract, 43, 435.02 [Google Scholar]
 Holman, M. J., Payne, M. J., Fraser, W., et al. 2018, ApJ, 855, L6 [NASA ADS] [CrossRef] [Google Scholar]
 Kuchner, M. J., Faherty, J. K., Schneider, A. C., et al. 2017, ApJ, 841, L19 [NASA ADS] [CrossRef] [Google Scholar]
 Lai, D. 2016, AJ, 152, 215 [NASA ADS] [CrossRef] [Google Scholar]
 Lawler, S. M., Shankman, C., Kaib, N., et al. 2017, AJ, 153, 33 [NASA ADS] [CrossRef] [Google Scholar]
 Lehner, M. J., Wang, S.Y., Alcock, C. A., et al. 2012, Proc. SPIE, 8444, 84440D [CrossRef] [Google Scholar]
 Li, G., & Adams, F. C. 2016, ApJ, 823, L3 [NASA ADS] [CrossRef] [Google Scholar]
 Linder, E. F., & Mordasini, C. 2016, A&A, 589, A134 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Malhotra, R., Volk, K., & Wang, X. 2016, ApJ, 824, L22 [NASA ADS] [CrossRef] [Google Scholar]
 Medvedev, Y. D., Vavilov, D. E., Bondarenko, Y. S., Bulekbaev, D. A., & Kunturova, N. B. 2017, Astron. Lett., 43, 120 [NASA ADS] [CrossRef] [Google Scholar]
 Meisner, A. M., Bromley, B. C., Kenyon, S. J., & Anderson, T. E. 2018, AJ, 155, 166 [NASA ADS] [CrossRef] [Google Scholar]
 Millholland, S., & Laughlin, G. 2017, AJ, 153, 91 [NASA ADS] [CrossRef] [Google Scholar]
 Nesvorný, D., Vokrouhlický, D., Dones, L., et al. 2017, ApJ, 845, 27 [NASA ADS] [CrossRef] [Google Scholar]
 Pass, E., Metchev, S., Brown, P., & Beauchemin, S. 2018, PASP, 130, 014502 [NASA ADS] [CrossRef] [Google Scholar]
 Petit, J.M., Holman, M., Scholl, H., Kavelaars, J., & Gladman, B., 2004, MNRAS, 347, 471 [NASA ADS] [CrossRef] [Google Scholar]
 Shankman, C., Kavelaars, J. J., Bannister, M. T., et al. 2017, AJ, 154, 50 [NASA ADS] [CrossRef] [Google Scholar]
 Sheppard, S. S., & Trujillo, C. 2016, AJ, 152, 221 [NASA ADS] [CrossRef] [Google Scholar]
 Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2006, AJ, 131, 1163 [NASA ADS] [CrossRef] [Google Scholar]
 Tombaugh, C. W. 1946, ASP Leaflets, 5, 73 [NASA ADS] [Google Scholar]
 Toth, I. 2016, A&A, 592, A86 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Trujillo, C. A., & Sheppard, S. S. 2014, Nature, 507, 471 [NASA ADS] [CrossRef] [Google Scholar]
 Weryk, R. J., Lilly, E., Chastel, S., et al. 2016, Icarus, submitted [arXiv:1607.04895] [Google Scholar]
 Wright, E. L., Eisenhardt, P. R. M., Mainzer, A. K., et al. 2010, AJ, 140, 1868 [NASA ADS] [CrossRef] [Google Scholar]
All Figures
Fig. 1.
Flowchart of the search algorithm. All process are displayed as square boxes, while all (intermediate) catalogs are marked by rounded boxes. The corresponding section in the text is given in the right column. 

In the text 
Fig. 2.
Schematic of the vectors introduced in Sect. 2.3.1. The Sun is denoted by ⊙, Earth by ⊕ and the TNO by ◯. 

In the text 
Fig. 3.
Makemake geocentric path between 2010 May 13 and 2015 April 7 generated using the HORIZONS WebInterface (blue line). All WISE observations during which Makemake was in the field of view are marked with a dot. The WISE detections of Makemake are marked with a red cross. 

In the text 
Fig. 4.
Illustration of the Makemake orbit in heliocentric coordinates. Transformation to heliocentric coordinates has removed the ellipsoidal component of the orbit and the epochal sets of observations (chronologically from upper right to lower left) are linearized to good approximation. All WISE observations during which Makemake was in the field of view are marked with a dot. The WISE detections of Makemake are marked with a red cross. 

In the text 
Fig. 5.
Squared Pearson correlation coefficient ρ^{2} as a function of the heliocentric distance guess d for a set of 37 calculated geocentric positions of Makemake between 01012017 and 01012018. The function peaks prominently at 52.5 au which is the average heliocentric distance of Makemake for the given time interval. 

In the text 
Fig. 6.
Squared Pearson correlation coefficients as functions of the heliocentric distance guess for different samples (1day cadence) of Makemake ephemeris, all covering 1 month in 2017. Despite almost identical time spans, numbers of observations and observation frequency, the ρ^{2} functions differ greatly. For such data, the fitting range [d_{min}, d_{max}] must closely bracket the actual physical peak. Otherwise, the maximization routine may end up stuck at the borders of the fitting range or at one of the secondary peaks. While this could be avoided by checking linearity for small steps in heliocentric distance, this would void the speed of the fitting routine; so instead, the analysis is carried out with steps of d_{max} − d_{min} = 10 au. 

In the text 
Fig. 7.
Distribution of the squared Pearson correlation coefficients of all subgroups of size k = 4 for the Makemake data cluster presented in Fig. 5. The vast majority of subgroups are wellfitted and yield good linearity values larger than 0.9999. 

In the text 
Fig. 8.
Light curve of Makemake in the WISE bands W1 (red crosses) and W2 (blue dots). The upper row shows the entire light curve, while the second row shows a zoom of the observation epochs in 1.5 d windows. Nondetections are marked with an arrow at the S/N = 5 level. We note that most detections have a S/N below 5. 

In the text 
Fig. 9.
Density of detected sources (sources/deg^{2}) of the input catalog. For clarity, the color code has been truncated at densities above 2000, since the regions at the ecliptic poles exhibit source densities almost two orders of magnitude above the average. The injected simulated Planet Nine observations are overplotted as black dots. 

In the text 
Fig. 10.
Number of clusters (solid line) and number of linearizations necessary to process the data (dotted line). The number of iterations converges toward ~10^{6} per cluster, an achievable value with reasonable computing power. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.