A 3π Search for Planet Nine at 3.4 μm with WISE and NEOWISE

Over the past fifteen years, observations have emerged that challenge the conventional understanding of our solar system's outer regions. Chief among these are the discoveries of Sedna (Brown et al. 2004) and 2012 VP₁₁₃ (Trujillo & Sheppard 2014), trans-Neptunian objects (TNOs) with orbits that cannot be explained by perturbations from the solar system's known planets (Gladman et al. 2002; Brown et al. 2004; Morbidelli & Levison 2004). These two objects are often grouped together with other high perihelion (q > 30 au), large semimajor axis (a > 250 au) bodies into a class of objects referred to as extreme TNOs (ETNOs). Thirteen such ETNOs are now known (Brown et al. 2004; Trujillo & Sheppard 2014; Sheppard & Trujillo 2016; Bannister et al. 2017). These objects have been found to cluster in argument of perihelion (Trujillo & Sheppard 2014) or longitude of perihelion (Batygin & Brown 2016a; Brown 2017). To explain the clustering in argument of perihelion, Trujillo & Sheppard (2014) proposed the existence of a super-Earth mass perturber at ∼250 au. In a related but different analysis, Batygin & Brown (2016a) suggested that an unseen giant planet ("Planet Nine") of mass 5–20 M_⊕ can simultaneously explain many enigmatic features of the solar system. These include the longitude of perihelion clustering among ETNOs, the existence of Sedna and 2012 VP₁₁₃, the population of high-inclination TNOs (Batygin & Brown 2016b), and the solar obliquity (Bailey et al. 2016; Lai 2016; Gomes et al. 2017). Criticisms of this Planet Nine narrative have also been put forth. For instance, the ETNO clustering in longitude of perihelion has been argued to result from observational bias (Shankman et al. 2017), and the presence of such a massive unseen planet has been claimed inconsistent with Cassini ranging measurements (Folkner et al. 2016). Lawler et al. (2017) and Nesvorný et al. (2017) did not find evidence that ETNO clustering in ω or perihelion longitude would result from the existence of an body like Planet Nine. Nesvorný et al. (2017) also suggested that Planet Nine would create a distribution of ecliptic comet inclinations in tension with that observed.

Nevertheless, the Planet Nine hypothesis yields concrete predictions (Brown & Batygin 2016) for the proposed planet's orbital parameters (380 ≲ a/au ≲ 980, i ∼ 30°, 0.3 ≲ e ≲0.75), facilitating searches that can directly observe or else place limits on its presence. Planet Nine is expected to be faint in the optical, but not prohibitively so (22 < V < 25; Brown & Batygin 2016). Studies have also suggested the possibility of detecting Planet Nine at a range of other wavelengths (Cowan et al. 2016; Fortney et al. 2016; Linder & Mordasini 2016), depending on its detailed physical properties. Owing in part to these encouraging forecasts, many archival and/or observational searches for Planet Nine are already underway (Fortney et al. 2016; Sheppard & Trujillo 2016; Weryk et al. 2016; Kuchner et al. 2017; Meisner et al. 2017a).

In this work, we focus on searching for Planet Nine in the infrared. Evolutionary models of Planet Nine spanning 5–60 M_⊕ predict that it could be self-luminous (Fortney et al. 2016; Linder & Mordasini 2016). In particular, Fortney et al. (2016) suggested that Planet Nine may be much brighter than blackbody expectations at 3–5 μm, based on a detailed analysis of model atmospheres. Fortney et al. (2016) found that Planet Nine may be detectable by the Wide-field Infrared Survey Explorer (WISE) in its bluest channel (W1; Wright et al. 2010), although the models span an enormous range of luminosities at mid-infrared wavelengths.

The W1-brightest model from Fortney et al. (2016) still remains fainter than the single-exposure detection limit, but would be readily detectable in coadds (at a fiducial distance of Planet Nine of d₉ ∼ 600 au). In Meisner et al. (2017a), hereafter M17, we developed a methodology to convert the WISE W1 imaging into a serendipitous search for Planet Nine using a custom set of time-resolved coadd images. We applied this procedure to a ∼2000 square degree region, which, at the time, was thought to represent a likely Planet Nine location (Holman & Payne 2016). Here, we extend our W1 search to cover more than three-quarters of the sky and include an additional year of WISE observations. Although several predictions have been made for Planet Nine's most likely current sky location (Brown & Batygin 2016; de la Fuente Marcos & de la Fuente Marcos 2016; Holman & Payne 2016; Millholland & Laughlin 2017), our present strategy is to uniformly search the widest area possible. The very high and homogeneous quality of the space-based imaging delivered by WISE is conducive to deriving detailed constraints on Planet Nine's apparent mid-infrared brightness in the event of a non-detection.

In Section 2, we describe the archival WISE data used in the course of this search. In Section 3, we review the Planet Nine search methodology presented in M17. In Section 4, we provide details of the footprint searched during the course of this work. In Section 5, we explain our orbit linking procedure. In Section 6, we discuss the sensitivity of our search. In Section 7, we describe the results of our search. We conclude in Section 8.

We begin by very briefly reviewing key aspects of the WISE mission and its timeline of observations. WISE is a 0.4 m telescope onboard a satellite in low-Earth orbit (Wright et al. 2010). The satellite was launched in late 2009. WISE has surveyed the entire sky in four broad bandpasses, centered at 3.4 μm (W1), 4.6 μm (W2), 12 μm (W3), and 22 μm (W4). WISE performed a seven-month full-sky survey in all four of its channels starting in 2010 January. In 2010 August and September, WISE surveyed the sky in its three bluest bands only. It then carried out an asteroid hunting mission called NEOWISE until 2011 February, operating only in the W1 and W2 channels (Mainzer et al. 2011). It then entered a hibernation period of nearly three years and recommenced survey operations in 2013 December. The post-hibernation mission phase has been dubbed NEOWISE-Reactivation (NEOWISER). WISE remains operational in W1 and W2, and has already completed roughly eight full sky passes throughout its first seven years in orbit.

NEOWISER publicly releases its data in chunks that each span one year of observations. These data releases consist of single-exposure images and source extractions and are published on a roughly annual basis. Thus far, three such data releases have been published, in 2015 March, 2016 March, and 2017 June.

A major difference between the present work and the initial search presented in M17 is our inclusion of an additional year of NEOWISER exposures obtained from the third-year NEOWISER data release, which became publicly available in 2017 June. In this work, we employ all publicly available W1 Level 1b (L1b) single-exposure images that fall within our search footprint described in Section 4. Our search footprint is defined primarily by cuts on ecliptic latitude (denoted by β, Section 4.1) and background source density (Section 4.2). In particular, we restrict our analysis to $| \beta | \lt 55^\circ$. The L1b acquisition dates range from 2010 January 7 to 2016 December 13 (UTC). In other words, we use all pre-hibernation W1 L1b frames, plus those frames included in the first three annual NEOWISER data releases. This amounts to ∼49 months of WISE observations. While the input data spans a roughly seven-year time period, the time baseline at a given sky location is typically 6.5 years.

We have downloaded⁵ local copies of all publicly available W1 L1b sky images (-int- files) along with their corresponding uncertainty maps (-unc- files) and bitmasks (-msk- files). We refer to one such group of three W1 files resulting from a single WISE pointing as a "frameset." The coadds presented in this work make use of ∼6 million such framesets, for a total of ∼40 TB and ∼18 terapixels of single-exposure inputs.

Our search strategy is modeled after Brown et al. (2015), although our analysis begins at the pixel level rather than the catalog level. Brown et al. (2015) developed a method for linking transients with arbitrary time sampling into Keplerian orbits for solar system objects beyond 25 au. We employ a similar orbit linking approach, which we describe fully in Section 5. The Brown et al. (2015) analysis begins with large catalogs of transient candidates provided by the Catalina Real-time Transient Survey (Drake et al. 2009). The starting point for our analysis is a set of ∼6 million single-exposure WISE/NEOWISE images. Therefore, much of our analysis is devoted to image processing steps that allow us to obtain a clean list of transients from the available WISE images. In brief, our method stacks W1 exposures into a set of epochal coadds, extracts transients via difference imaging, and finally attempts to link these transients into Keplerian orbits beyond 250 au. The following subsections briefly review each step of our Planet Nine search procedure. Full details are provided in M17.

3.1. Tiling

Our search region (Section 4) is composed of a series of "tiles", with each tile corresponding to a unique ∼2.5-square-degree astrometric footprint. These are the footprints onto which we reproject L1b exposures during coaddition. We adopt tile footprints matching those of the WISE team's Atlas stacks (Cutri 2011). The Atlas tiling covers the entire sky with a set of 18240 astrometric footprints. Each Atlas tile footprint is 156 on a side and is identified by a coadd_id value. The coadd_id is a string encoding the tile's central (R.A., decl.) coordinates. For instance, the Atlas tile centered at (R.A., decl.) = (3248, −197) has coadd_id = "3248m197." The Atlas tile centers trace out a series of isodeclination rings. For our time-resolved coadds, we adopt a pixel scale matching that of the L1b images, 275/pixel. Thus, each of our tiles has dimensions of 2048 pixels × 2048 pixels. In our terminology, "tiles" are distinct from "coadds": there are exactly 18240 unique tile footprints, whereas we have typically constructed 8–9 epochal coadds for each such tile footprint within our search region (see Table 1). Note that although the 18240 astrometric footprints of the Atlas tiles are unique, neighboring tiles overlap at their boundaries. The amount of overlap is ∼3' near the celestial equator.

Table 1.  Fraction of Pixels in Our Search Footprint Region as a Function of Number of Coadd Epochs

nepochs	$p({n}_{\mathrm{epochs}})$
0	3.08 × 10−6
1	1.13 × 10−6
2	2.90 × 10−6
3	3.98 × 10−6
4	6.30 × 10−6
5	7.98 × 10−6
6	8.13 × 10−5
7	5.35 × 10−2
8	3.08 × 10−1
9	3.85 × 10−1
10	2.21 × 10−1
11	3.23 × 10−2
12	4.62 × 10−6

Download table as:  ASCIITypeset image

3.2. Time Slicing

3.3. Coaddition

We employ an adaptation of the unWISE code (Lang 2014) to generate W1 coadds which incorporate the time slicing described in Section 3.2 and retain the native WISE angular resolution. Complete details of the unWISE code and the particular options we use in the context of our W1 Planet Nine searches can be found in Lang (2014), Meisner et al. (2017c) and M17. The number of coadd epochs generated per sky location is shown in the top panel of Figure 1. In total, our present search incorporates 129,483 time-resolved W1 coadds. For comparison, the bottom panel of Figure 1 shows the number of coadd epochs per sky location in the M17 search, which included only 6219 such coadds. The median, mode, and mean of the number of coadd epochs per sky location are 9, 9, and 8.87, respectively.

Zoom In Zoom Out Reset image size

Our analysis includes approximately eight full sky passes, plus an additional ∼1/6 of a sky pass from early 2011. ∼5% of the sky is missing one of the usual eight passes due to a WISE spacecraft "command timing anomaly", which halted science observations in 2014 April. As a result, under an idealized survey strategy in which WISE always points at exactly 90^◦ solar elongation, we would expect ∼80% of the sky to have 8 coadd epochs, ∼15% to have 9 coadd epochs, and ∼5% to have 7 coadd epochs. However, Moon avoidance maneuvers lead to extra coadd epochs, with each sky location typically affected by one such maneuver. This results in roughly one "bonus" coadd epoch everywhere on the sky. It is possible that the same sky region has been affected by more than one Moon avoidance maneuver, explaining the Table 1 tail of n_epochs, which extends up to 12 over a fractionally small area. The small fraction of area with ${n}_{\mathrm{epochs}}\leqslant 6$ is due to rare circumstances under which our coaddition fails, for example, in the heavily saturated cores of a few of the brightest stars or during a sighting of Mars, Jupiter, or Saturn.

3.4. Reference Templates

3.5. Image Differencing

3.6. Extracting Transients

3.7. Filtering Transients

We go to great lengths to eliminate identifiable artifacts from our transient catalogs while masking the smallest possible fraction of sky area. The transients we remove include those associated with:

• 1.  
  Static sources, as identified by 2MASS counterparts.
• 2.  
  Kernel matching artifacts. These are ringing features introduced when HOTPANTS attempts to establish consistency between the search and reference image PSFs.
• 3.  
  The cores and/or diffraction spikes of bright star profiles.
• 4.  
  Latents. These are persistence artifacts that appear at the detector position of saturated stars in subsequent exposures. Latents have the appearance of diffuse smudges and are detected as spurious transients because their positions vary as the WISE scan direction changes. Given a catalog of bright star locations, it is possible to very accurately predict the positions of latents. We use the AllWISE catalog to identify bright stars capable of causing latents.

In M17, we created custom masks for four "ultra-bright" stars whose profiles extend far beyond the PSF model used in our standard bright star masking approach. In this work, we opted not to continue the practice of making such custom masks for bright stars throughout our search footprint. This choice did not lead to a fractionally significant increase in spurious transients, as the transient catalog is dominated by artifacts associated with much more numerous moderately bright sources, especially in relatively low Galactic latitude regions of our search footprint.

At high $| {b}_{\mathrm{gal}}|$ ($| {b}_{\mathrm{gal}}| \gt 30^\circ$), our filtering masks 1.0% of the area. This fraction varies dramatically over the sky, mainly due to variations in the number of static background sources within the Galaxy along different sightlines (see Figure 2). In Section 6, we account for spatial variations of the fractional area masked when assessing our survey's sensitivity.

Zoom In Zoom Out Reset image size

3.8. Orbit Linking Overview

Our linking approach is similar to that of M17 but different enough in detail that we have fully re-described it in Section 5. Our orbit linking procedure employs the orbfit code of Bernstein & Khushalani (2000) to test whether groupings of transients ("tuples") can be well-fit by any bound Keplerian orbit in the outer solar system. This orbit linking methodology is modeled after that of Brown et al. (2015) and is capable of linking transients with arbitrary cadence. Notably, in this work we require that a Planet Nine "candidate" consist of at least five linked transients, versus only four required in M17. We found this choice to be important in reducing the number of false positive linkages, particularly in crowded regions where the number density of transients spikes dramatically. Because we have incorporated an additional year of W1 frames relative to M17, we suffer no loss in sensitivity despite requiring a minimum of five rather than four linked transients (see Figure 4).

Zoom In Zoom Out Reset image size

The astrometric uncertainties assigned to our transient centroids represent a critical input for orbit linking, as these dictate the χ² values returned by orbfit. We have sought to be very conservative/generous in our centroid uncertainty model, which is described in Section 7.1 of M17. The typical per-coordinate uncertainty we adopt is 08, nearly a third of a WISE pixel. As in M17, our choice of χ² threshold is guided by our desire for a false negative rate of less than 0.001. For quintuplet linkages, this translates to a requirement of ${\chi }^{2}\lt 21.5$. In practice, to be conservative, we have inspected finder charts for all candidates with five or more linked transients and χ² < 40. This can be viewed, in the case of quintuplet linkages, as effectively padding our centroid uncertainties by an additional ∼40%.

Although WISE observes the entire sky, some regions have been omitted from our search. The following subsections detail the circumstances under which we have trimmed area. Our search footprint includes 76.3% of the sky.

4.1. High Ecliptic Latitude

4.2. The Galactic Plane

We run all stages of our processing up to and inclusive of transient filtering for all 15058 tile footprints satisfying our ecliptic latitude cut. However, we find that an excessively high number density of transients remains in our filtered catalogs in low Galactic latitude regions toward the inner Galaxy. The enhanced transient number density near the Galactic center is due to the high density of static background sources in this region. Extreme, localized increases in transient density are particularly problematic because of the combinatorics associated with orbit linking—the number of transient n-tuples increases with the nth power of transient number density. We therefore sought to excise regions near the Galactic center containing huge numbers of static background sources. We do so by making a cut on the number of 2MASS (Skrutskie et al. 2006) sources per tile footprint, n_2mass. We remove tiles with ${n}_{2\mathrm{mass}}$ > 90000. This corresponds to one 2MASS source per 47 pixels. For comparison, the effective number of pixels corresponding to the W1 PSF is 15.8, meaning that our chosen source density threshold translates to roughly 0.34 2MASS sources per WISE beam. This is equivalent to 2.8 × 10⁻³ 2MASS sources per square arcsecond. We use 2MASS source density as a proxy for crowding, rather than WISE source density, because 2MASS has far superior angular resolution. In general, sources detected by 2MASS are also detected by WISE (WISE is deeper), though deblending issues may modulate/reduce the number of WISE sources actually cataloged in the densest regions.

1131 tiles are rejected due to our cut on number of 2MASS sources. These tiles trace out a roughly diamond-shaped area centered at (l_gal, b_gal) = (0°, 0°) when viewed in Galactic coordinates. This region is shown as the empty (light blue) stripe passing through the Galactic center at (λ, β) = (267^◦, −6°) in Figures 1, 2, 5, and 7. The region excised due to high source density reaches a maximum $| {b}_{\mathrm{gal}}|$ of ∼14^◦ near ${l}_{\mathrm{gal}}=0^\circ$, then tapers toward the Galactic plane, disappearing by l_gal = 90° toward the east and l_gal = 282° toward the west. We are able to perform our search all the way into the Galactic plane (through ${b}_{\mathrm{gal}}=0^\circ$) in the outer Galaxy—our background source density cut is lenient enough to retain most of the outer Galactic plane.

Zoom In Zoom Out Reset image size

We also found that a region of the Galactic plane near (l_gal, b_gal) = (109°, −2°) exhibited an unusually high number density of transients associated with diffuse emission. We presume this to result from time-variable illumination of dust clouds (e.g., Rest et al. 2008). As a result, we removed eight nearby tiles, with coadd_id values 3461p575, 3484p590, 3489p575, 3510p605, 3513p590, 3516p575, 3540p605, and 3542p590.

4.3. External Galaxies

Our orbit linking strategy follows the procedure we outlined in M17, with some modification to handle the increased number of epochs and transients. First, for each transient, we identify a pool of all nearby subsequent transients that might be associated with bound Keplerian orbits at distances of 250 au or more. Next we construct triplets that consist of the first transient and pairs of nearby transients from the pool around it. We then test if each triplet is consistent with a simple model of reflex motion plus linear drift in the sky plane, and we accept it if its χ² is less than 30. We repeat this process for quadruplets built from accepted triplets and a fourth transient from the pool. Quadruplets themselves are tested with the parallax-linear drift model, but with a χ² threshold of 150. These steps, which convert the full list of transients into potential linkages of quadruplets, constitute our pre-screening process.

We test the pre-screening process with mock data. By sampling Keplerian orbital parameters for hypothetical bound solar system objects beyond 250 au, we generate mock transients with cadence and astrometric uncertainties similar to the real data. Applying the same parallax and linear drift tests on triplets and quadruplets as we do for real data, we find that our pre-screen step falsely rejects no more than 1 in 10⁴ linkages.

In the next step of the linkage procedure, we evaluate each pre-screened quadruplet with orbfit, a code designed by Bernstein & Khushalani (2000) to fit observations of outer solar system objects with Keplerian orbital elements. We choose a χ² threshold of 40 in the decision to accept or reject a quadruplet, irrespective of order or degrees of fitting (orbfit uses several different models, depending on the details of the transients). This choice yields false negative rates that are well below the desired maximum value of 0.001 adopted in M17. We began the present search using χ² thresholds yielding false negative rates of exactly 0.001, to match M17, but expanded our analysis to higher χ² values in attempt to be conservative. The actual value of 40 is not special, and was chosen simply because it corresponds to a false negative rate below 0.001 while yielding an ambitious yet tractable number of candidates requiring finder chart inspection.

Finally, we check if the quadruplets that are admitted on the basis of orbfit have any transients in common. In cases where they do, we merge them together to form quintuplets, or even higher-order tuples, if the merged orbfit χ² is less than 40, as before. For all of the linkages that we find, the number of degrees of freedom is low enough that this χ² threshold preserves a false negative rate below the 0.001 threshold.

This procedure culls over 5 million transients into roughly 2500 quintuplets, many of which are subsets of nearly 200 higher-order n-tuples (n = 6–8; see Section 7). Along the way, the pre-screen process enumerates over 92 billion triplets and 275 million quadruplets. Our code ends up calling orbfit to check over 8 million potential linkages.

Our completeness estimates are based on an empirical detection efficiency measured from 300,000 fake sources we have injected into our real difference images. We measure the fraction of these fake sources recovered as a function of W1 magnitude and number of exposures (n_exp) contributing to the relevant epochal coadd pixels. Using the probability distribution functions for both n_exp and the number of coadd epochs per sky location (n_epoch), we then perform Monte Carlo simulations to assess the probability with which five or more transients will be recovered as a function of an object's W1 magnitude. This latter probability is what we refer to as our completeness, as our search demands that a Planet Nine candidate consist of a low-${\chi }^{2}$ linkage containing at least five transients. Our Monte Carlo simulations also account for the fraction of area masked (f_masked) by randomly zeroing out the appropriate (spatially varying) fraction of single-epoch detection probabilities.

The probability distribution function for n_epoch is provided in Table 1. We have also computed the probability distribution function of n_exp for the 129,483 tiles in our present search and show the results in Figure 3. Both the median and mode are ${n}_{\exp }=12$.

Zoom In Zoom Out Reset image size

Figure 4 shows the results of our Monte Carlo simulations at high Galactic latitude (${f}_{\mathrm{masked}}=0.01$), plotting the fraction of objects with ≥5 transient detections as a function of W1 magnitude. 90% (50%) completeness is achieved at W1 = 16.70 (16.86). Section 8 of M17 provides additional details of our fake source injection (Section 8.3) and Monte Carlo simulation procedures (Section 8.5).

Because our search spans a wide range of background source densities and includes sub-regions with different n_epoch values, we have also produced a map of the 90% completeness threshold's spatial variation. In doing so, we perform additional Monte Carlo simulations to simultaneously account for five factors:

• (i)  
  The local value of n_epoch as provided by the map shown in Figure 1. Regions with more coadd epochs will be more sensitive and therefore reach 90% completeness at higher W1 magnitudes.
• (ii)  
  The local value of f_masked as given by the map shown in Figure 2. Regions with larger f_masked have reduced sensitivity.
• (iii)  
  Additional background noise in regions near the Galactic plane, presumably due to large numbers of static sources, both resolved and unresolved, as well as diffuse Galactic emission. We find that this effect is well-approximated by a reduction in depth of:

  $$\begin{eqnarray}{{\rm{\Delta }}}_{\mathrm{bg}} & = & -1.08\times {10}^{-2}+(2.76\times {10}^{-6}){n}_{2\mathrm{mass}}\\ & & +(5.94\times {10}^{-12}){n}_{2\mathrm{mass}}^{2}\ [\mathrm{mag}].\end{eqnarray} \tag{ 1 }$$

  This effect is largest in the most crowded regions of our footprint, reaching a maximum value of 0.29 mag.
• (iv)  
  Variation in the amount of parallactic smearing as a function of ecliptic latitude, resulting from the fact that the time span of a single coadd epoch scales as 1/cos(β). This tends to reduce the sensitivity of our search at high $| \beta |$ relative to low $| \beta |$.
• (v)  
  The average increase of n_exp as 1/cos(β). This effect tends to increase the sensitivity at high ecliptic latitude relative to that at low ecliptic latitude.

Figure 5 shows the variation in 90% completeness threshold implied by taking these effects into account, for a fiducial distance of d₉ = 700 au. The map is a function of d₉ because of the parallax-dependent term described in item (iv) above. The vertical stripes in Figure 5 trace variations in n_epoch as a function of ecliptic longitude. Reduced sensitivity toward the Galactic plane is also apparent, due to the effects described in items (ii) and (iii) above. The ecliptic latitude dependent effects in items (iv) and (v) tend to cancel each other out and are not especially obvious in the Figure 5 map. The mean 90% completeness threshold within our search footprint is W1 = 16.69, and the 25th (75th) percentile value is W1 = 16.64 (W1 = 16.76). The full range of 90% completeness threshold values spanned is quite broad, 15.97 < W1 < 16.93. The map shown in Figure 5 is available online.¹⁰

6.1. Cautionary Notes on Completeness

5,603,879 transients remain after filtering our difference image detections as described in Section 3.7. Figure 6 shows the distribution of W1 magnitudes for our transients at high Galactic latitude. As expected, this distribution peaks near W1 = 16.6, which is ∼1.3 mag below the single-exposure detection limit at W1 = 15.3. The spatial distribution of these transients is shown in Figure 7.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Running our orbit linker on the full set of transients results in 2735 linkages with orbfit χ² < 40 and a minimum of five transients. 2560 of these linkages are quintuplets, 163 are sextuplets, 11 are septuplets, and one is an octuplet. The locations of all 2560 low-χ² quintuplet linkages are indicated by red dots in Figure 7. Because WISE scans are oriented vertically in this coordinate system, essentially all relevant quantities (such as the number of available coadd epochs) and systematics (such as periods of relatively high Moon glow) will imprint as vertical stripes. The transient density and therefore the density of false positive quintuplets also trace these vertical stripes. For instance, the vertical stripes in Figure 1 will cause vertical striping in the number of transients per unit area, simply because certain stripes have more or fewer coadd epochs available.

We have examined finder charts for all 2735 (χ² < 40) candidates. We note that this χ² threshold is quite conservative in terms of its corresponding false negative rate. For instance, had we implemented cuts to achieve a formal false negative rate of 0.001, only 1296 linkages would require finder chart inspection.

In scrutinizing the finder charts, we seek to determine whether any linkage exists for which ≥5 of its transients are legitimate, in the sense that they cannot be ruled out as associated with artifacts or previously known moving objects. Because we do not impose any cuts on the orbital elements returned by orbfit, it turns out that many of our linkages are known high proper motion objects in the solar neighborhood. We therefore downloaded the list of all SIMBAD (Wenger et al. 2000) proper motion sources with $| \mu | \gt 200$ mas yr⁻¹, and used this catalog to automatically flag linkages associated with these known moving objects, employing a match radius of 15''. 737 of our 2735 linkages are associated with known high proper motion objects cataloged by SIMBAD. We exercised caution by checking each such finder chart to ensure a legitimate association, and verified that the direction of proper motion seen in our finder chart matched that quoted by SIMBAD. A small number of our linkages displayed linear motion consistent with being stars in the solar neighborhood, and appeared to be new discoveries upon searching SIMBAD and Vizier (Ochsenbein et al. 2000). In these cases we checked 2MASS data to verify the linear nature of the motion on longer timescales.

For all remaining 1998 linkages, we were able to rule out a sufficient number of transients such that four or fewer legitimate transients remained. The spurious transients that we ruled out can be grouped into three dominant categories:

• 1.  
  Latents. It is difficult to automatically flag 100% of latents for several reasons. For instance, we found a number of cases in which the parent bright star is entirely absent from the AllWISE catalog from which we constructed a sample of potential latent parents. In a few similar cases, the parent bright star is present in AllWISE yet has a W1 magnitude listed as "null." In many cases the parent's AllWISE W1 magnitude is slightly fainter than the brightness threshold we employed in constructing our sample of latent parent stars. Latents of third order or higher often contaminated our linkages. Such latents occur ≥3 exposures after imaging of the parent star, whereas our automated masking of latents only accounts for first and second order latents.
• 2.  
  Diffraction spikes. Our bright star masks are highly detailed and do account for diffraction spikes, but these masks are based on a W1 PSF model of angular size 149. Even though this PSF model extends quite far into the wings, the diffraction spikes of exceptionally bright stars can remain detectable far beyond the angular extent of our masks.
• 3.  
  Known high proper motion objects that evaded our SIMBAD crossmatch. This includes cases where a known proper motion counterpart was found in Vizier but not SIMBAD and extremely fast-moving objects ($| \mu | \gtrsim 1\buildrel{\prime\prime}\over{.} 5\,{\mathrm{yr}}^{-1}$) for which our 15'' SIMBAD match radius was insufficient.

Other relatively common contaminants included spurious transients associated with bright star glints, real astrophysical variability of static background sources and kernel matching artifacts which evaded our filtering procedure.

One might view the higher-order linkages (sextuplets and higher) as the most interesting candidates, as in these cases the number of linked transients exceeds our minimum value of five. However, we found that all sextuplets, septuplets, and the one octuplet were especially clear cut cases in which sufficient numbers of transients could easily be ruled out as known high proper motion objects or bright star artifacts. In particular, the octuplet consisted of eight detections of the white dwarf 2MASS J14205484−0905086.

As we do not find any valid linkages of ≥5 transients, we are left with a non-detection of Planet Nine and the completeness constraints described in Section 6.

We have performed the deepest and widest-area WISE-based search for Planet Nine, covering over three-quarters of the sky. At high Galactic latitude, we rule out the presence of Planet Nine in the parameter space searched with W1 < 16.7 (16.86) at 90% (50%) completeness. The sensitivity of our survey is reduced toward the Galactic plane, as illustrated in Figure 5. Our search methodology is valid for $250\lt {d}_{9}/\mathrm{au}\lt 2250$, although the lower distance boundary is dependent on ecliptic latitude. For the most W1-luminous Fortney et al. (2016) model (which has a mass of 10 M_⊕), our search probes distances up to 800–900 au depending on ecliptic latitude. The optical-W1 color of Planet Nine is extremely model-dependent and is further complicated by the fact that intrinsic emission in W1 becomes fainter as ${d}_{9}^{-2}$ whereas reflected sunlight in the optical becomes fainter as ${d}_{9}^{-4}$. For the maximally W1-luminous Fortney et al. (2016) model, our W1 fainter than 16.7 constraint corresponds to a limit of VR fainter than 22.1 at a fiducial distance of 650 au.

Our methodology should be readily applicable to the WISE W2 channel. Typically, (W1–W2) ≈ 0.8 corresponds to equal signal-to-noise in W1 and W2, meaning that we should expect sensitivity to Planet Nine at W2 < 15.9. This survey limit would render the most W2-luminous model of Fortney et al. (2016) detectable at distances of up to ∼1700 au. One operational advantage of searching in W2 rather than W1 is that W2 has a factor of ∼2 fewer background sources relative to W1. A disadvantage is that scattered moonlight, one of the dominant time-dependent artifacts affecting WISE data, is more pronounced in W2 than in W1.

Whereas our survey has omitted some regions, we note that the ongoing Kuchner et al. (2017) search is examining the entire sky in both W1 and W2, based on the time-resolved coadds of Meisner et al. (2017b). Because the Kuchner et al. (2017) search employs deep W2 coadds, it can also serve as a fainter extension of the Luhman (2014) search for a roughly Jovian mass companion to the Sun. WISE continues to collect additional data in W1 and W2, so that there will eventually be at least one additional year of NEOWISER exposures in both bands relative to those currently available.

This work has been supported by grant NNH17AE75I from the NASA Astrophysics Data Analysis Program, the NASA Outer Planets Program through grant NNX11AM37G and Emerging Worlds program grant NNX17AE24G. Resources supporting this work were provided by the NASA High-End Computing (HEC) Program through the NASA Center for Climate Simulation (NCCS) at Goddard Space Flight Center. We thank Eddie Schlafly for helpful comments on this manuscript.

This research makes use of data products from the Wide-field 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 also makes use of data products from NEOWISE, which is a project of the Jet Propulsion Laboratory/California Institute of Technology, funded by the Planetary Science Division of 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 National Energy Research Scientific Computing Center, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231, provided staff, computational resources, and data storage for this project.