SEARCHING FOR PLANET NINE WITH COADDED WISE AND NEOWISE-REACTIVATION IMAGES

WISE surveys the sky by remaining pointed at very nearly 90° solar elongation, completing one 360^◦ scan per orbit while approximately following a great circle of constant ecliptic longitude. As a result, away from the ecliptic poles, a given sky region will intersect the WISE field of view at six monthly intervals, each time remaining in the 078 FOV for ∼1 day, during which one exposure is acquired every ∼95 minutes. We refer to such day-long periods of WISE observation as "visits."

Here we search for Planet Nine by creating per-visit inertial coadds of W1 single exposures. Because we use data from all phases of the WISE+NEOWISER mission, we typically have seven such epochal coadds available per sky location within our search footprint. This number of coadd epochs results from the fact that WISE has visited our search footprint during seven distinct "seasons," which we label with letters a–g, as defined in Table 1. A linkage of transients at four epochs along a Keplerian orbit is generally considered the minimum threshold to identify a TNO candidate. Therefore, our ∼7 coadded W1 epochs are sufficient to identify or place meaningful constraints on the presence of Planet Nine candidates.

Table 1 also includes a "parity" column, which essentially serves as a binary indicator of Earth's location relative to the Sun during each season. Because WISE always points at the very nearly solar elongation of 90^◦, observations of a particular sky location with the same parity have almost exactly zero relative parallax, whereas observations with opposite parity have relative parallax of approximately 2 au.

The foregoing depiction of WISE's survey strategy is overly simplistic because it ignores Moon-avoidance maneuvers, during which WISE can point many degrees away from elongation 90°. Our handling of these events is discussed in detail in Section 5.2. The remainder of this section's formulae and forecasts will ignore the complication of Moon-avoidance maneuvers, although in practice these actually provide us with the benefit of extra coadd epochs. As a result, we have on average 7.3 epochal coadds per sky location in the search region studied.

Apparent motion due to parallax during the course of a coadd epoch places restrictions on the range of distances at which a TNO will be detectable in our coadds, since moving sources become smeared out over an area larger than that of a point source, effectively decreasing their signal-to-noise ratio. To determine the extent of parallactic smearing, we must first calculate the expected timespans of our epochal coadds.

Consider an idealized version of the WISE survey strategy, assuming WISE always points at exactly 90° solar elongation. At the ecliptic latitude (for which we use the symbol β) of zero, a given sky location will be observed over a period of 0776 × (365.25 days/360°) = 0.787 days during a WISE visit. In reality, WISE does not point at precisely 90° elongation, but toggles back and forth by a few tenths of a degree on even/odd orbits.⁷ This effectively lengthens the β = 0° visit duration to 1.0 ± 0.05 days, based on examination of per-pixel maps of minimum and maximum MJD for real coadds.

WISE obtains a constant number of frames per unit ecliptic latitude, and therefore the number of exposures per unit area within a visit scales as 1/cos($| \beta |$), as does the time spanned by exposures overlapping a given point on the sky:

$$\begin{eqnarray}&&{\tau }_{\mathrm{coadd}}=\displaystyle \frac{1.0\,\mathrm{days}}{\cos (| \beta | )}.\end{eqnarray} \tag{ 1 }$$

Ignoring the curvature of the Earth's orbit, Planet Nine's parallax during a single visit of WISE observations is

$$\begin{eqnarray}&&{\pi }_{9}=\displaystyle \frac{{{dl}}_{\perp }}{{d}_{9}}=\displaystyle \frac{(2\pi \,\mathrm{au})\times \tfrac{{\tau }_{\mathrm{coadd}}}{365.25\,\mathrm{days}}\times \sin (| \beta | )}{{d}_{9}},\end{eqnarray} \tag{ 2 }$$

where d₉ is the distance to Planet Nine in astronomical units and the sin($| \beta |$) factor projects the Earth's motion onto the direction perpendicular to the line of sight. The combination of the previous two equations can be rewritten as

$$\begin{eqnarray}&&{\pi }_{9}=5\buildrel{\prime\prime}\over{.} 91\times \displaystyle \frac{600\,\mathrm{au}}{{d}_{9}}\times \tan (| \beta | ).\end{eqnarray} \tag{ 3 }$$

The smearing caused by parallactic motion means that there will be a minimum distance at which our coadd-based search is sensitive to Planet Nine. For instance, even in the absence of orbital drift, an object in the Kuiper Belt at $| \beta |$ = 15° and 35 au from the Earth would appear to move by ∼27'' within a coadd epoch, making it very highly streaked. On the other hand, at a fiducial distance of 600 au and the same $| \beta |$, Planet Nine would experience only 16 of parallactic motion, rendering its W1 profile very nearly pointlike.

To quantify and account for parallactic smearing, we compute the effective number of noise pixels of a TNO profile as a function of parallax within a coadd epoch. The effective number of noise pixels for a particular profile is given by

$$\begin{eqnarray}&&{n}_{\mathrm{eff}}=\displaystyle \frac{{\left(\sum {f}_{i}\right)}^{2}}{\sum {f}_{i}^{2}},\end{eqnarray} \tag{ 4 }$$

where f_i is the flux in each pixel i. We simulate Planet Nine profiles that translate by an amount ${\pi }_{9}$ during the course of a coadd epoch with 12 exposures, for 0'' < π₉ < 825. For these simulations, we use the W1 PSF constructed according to Meisner & Finkbeiner (2014), as shown in Figure 4 of Lang (2014). The resulting n_eff values are well fit by

$$\begin{eqnarray}\displaystyle \frac{{n}_{\mathrm{eff}}}{{n}_{\mathrm{eff}}({\pi }_{9}=0)} & = & 1-(1.42\times {10}^{-3}){\pi }_{9}\\ & & +\ (9.59\times {10}^{-3}){\pi }_{9}^{2}-(2.91\times {10}^{-4}){\pi }_{9}^{3},\end{eqnarray} \tag{ 5 }$$

with ${\pi }_{9}$ in arcseconds. 825 corresponds to 3× the native WISE pixel size, and 1.35× the W1 FWHM, but even so n_eff(${\pi }_{9}=8\buildrel{\prime\prime}\over{.} 25$) remains less than 1.5× that of a W1 point source. We choose ${\pi }_{9}=8\buildrel{\prime\prime}\over{.} 25$ as a threshold beyond which we will not claim any sensitivity of our search, since we do not aim to push into the regime of detecting/photometering highly elongated sources. We nevertheless believe that we have considerable sensitivity to less pointlike sources (see Section 6.3.7). We note that a typical galaxy (r_half = 045) in 1'' seeing has n_eff 1.99× larger than that of a Gaussian PSF—by restricting to ${\pi }_{9}$ ≤ 825, we remain strongly within the quasi-pointlike regime for which standard source detection and photometry techniques are applicable.

Our per-visit coadds very closely resemble the WISE Preliminary Release Atlas stacks in terms of depth of coverage and sensitivity. In unconfused regions of the ecliptic plane with typical integer coverage (12 exposures), the preliminary release 5σ detection limit for static point sources is ${\rm{W}}{1}_{\mathrm{lim},\mathrm{PSF}}=16.5$ (Cutri et al. 2011). We always quote WISE magnitudes in the Vega system.⁸ We can calculate W1_lim($| \beta |$, d₉) by scaling this static point source limiting magnitude to account for parallactic smearing via its influence on n_eff($| \beta |$, d₉), as well as the increasing number of exposures with increasing $| \beta |$:

$$\begin{eqnarray}{\rm{W}}{1}_{\mathrm{lim}}(| \beta | ,{d}_{9}) & = & {\rm{W}}{1}_{\mathrm{lim},\mathrm{PSF}}\,+\,{{\rm{\Delta }}}_{\mathrm{ref}}-2.5{\mathrm{log}}_{10}\\ & & \times \left[\sqrt{\displaystyle \frac{{n}_{\mathrm{eff}}(| \beta | ,{d}_{9})}{{n}_{\mathrm{eff}}({\pi }_{9}=0)}}\sqrt{\displaystyle \frac{{n}_{\exp }(\beta =0)}{{n}_{\exp }(| \beta | )}}\right],\end{eqnarray} \tag{ 6 }$$

where n_exp is the number of W1 exposures contributing to a single sky location within an epochal coadd. Allowing n_exp to represent an average capable of taking on fractional values, we adopt

$$\begin{eqnarray}&&{n}_{\exp }(| \beta | )=\displaystyle \frac{{n}_{\exp }(\beta =0)}{\cos (| \beta | )}.\end{eqnarray} \tag{ 7 }$$

We will search for transients in difference images, rather than the epochal coadds themselves (Section 6). The ${{\rm{\Delta }}}_{\mathrm{ref}}$ term represents a small loss in sensitivity to due to the noise added by subtracting the reference coadd, which typically has 6× more contributing exposures than the epochal coadd. In that case,

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{\mathrm{ref}}=-2.5{\mathrm{log}}_{10}(\sqrt{(1+1/6)})=-0.084\,\mathrm{mag}.\end{eqnarray} \tag{ 8 }$$

We can now evaluate ${\rm{W}}{1}_{\mathrm{lim}}(| \beta | ,{d}_{9})$ and compare against Planet Nine's minimum W1 apparent magnitude as a function of d₉, yielding a mask for accessible versus inaccessible regions of ($| \beta |$, d₉) parameter space. We calculate the minimum W1 magnitude as a function of d₉ by assuming a maximum W1 luminosity of H_W1 = 2.13 (W1 = 16.1 at d₉ = 622 au; Fortney et al. 2016, Table 2). The results are shown in Figure 1.

Zoom In Zoom Out Reset image size

Table 2.  Moon-avoidance Maneuvers

Season	MJDstart	MJDend	${\lambda }_{\mathrm{approx}}$
a	55219.5	55220.075	380
b	55381.28	55381.88	120
b	55410.91	55411.53	410
c	55573.63	55574.18	265
d	56665.1	56665.62	200
d	56694.815	56695.27	500
e	56856.71	56857.16	230
e	56886.10	56886.55	510
f	57019.74	57020.21	95
f	57049.49	57049.97	430
g	57211.67	57212.08	160
g	57240.87	57241.335	435

Download table as:  ASCIITypeset image

Per the discussion in Section 3.3, regions with ${\pi }_{9}$ > 825 are marked as inaccessible to our search in Figure 1. The d₉ < 250 au region of parameter space is also marked as inaccessible for multiple reasons. Planet Nine could only be so nearby if very close to perihelion, which is unlikely (Brown & Batygin 2016). Additionally, the total apparent motion over 5.5 years can be larger for smaller d₉, which requires fitting orbits to very large numbers of transient tuples. As a result, purely by chance, we expect to become disproportionately flooded with spurious linkages as we push toward lower d₉. Because the marginal benefit of pushing to d₉ < 250 au is outweighed by the increase in false positives, we perform our search only for d₉ $\geqslant \,250\,\mathrm{au}$.

In summary, we expect substantial sensitivity to the maximally W1-luminous Fortney et al. (2016) model, especially at relatively low $| \beta |$. For instance, at $| \beta |$ ≤30° the H_W1 = 2.13 model is detectable at all distances between 250 and ∼700 au. We also forecast sensitivity out to d₉ $\approx \,800\,\mathrm{au}$ at high $| \beta |$. We note that during orbit linking (Section 7), no upper limit is imposed on d₉. We are therefore able to detect Planet Nine out to thousands of astronomical units, provided it could be sufficiently more luminous than any of the Fortney et al. (2016) models. Our actual sensitivity turns out to be ∼0.2 mag better than the present forecasts, essentially because of the distinction between the 5σ detection threshold considered in this section and our 90% completeness threshold for acquiring a quadruplet of detections (Section 8).

Alternative search strategies based on linking single-exposure detections (for instance a W1 variant of Luhman 2014) are not expected to be sensitive beyond d₉ = 430 au, given the single-exposure limiting magnitude of W1 = 15.3 and adopting the maximum W1 luminosity from Fortney et al. (2016). Still, a search based on single-exposure transients would be complementary to our coadd-based analysis, filling in some of the low d₉ parameter space, which is problematic for our approach.

We have thus far ignored the complication of orbital drift within the timespan of a coadd epoch. If large enough, orbital drift could compromise our sensitivity by further smearing the appearance of Planet Nine in our coadds. To determine the maximum intra-coadd orbital drift, we calculate the maximum intra-coadd drift within the mask shown in Figure 1 for each (a₉, e₉) pair in a set of orbits tracing the two loci of a₉ versus e₉ provided in Section 2 of Brown & Batygin (2016). We find that the maximum intra-coadd drift across all orbits is 129, which occurs near $| \beta |$ = 30^◦, for orbits currently very close to perihelion and with perihelia very similar to the inner edge of our search space (q ≈ d₉ ≈ 250 au).

Thus, the intra-coadd orbital drift is always expected to be less than one-fourth of the W1 PSF FWHM. Because of the small amplitude of the intra-coadd drift, and because it does not in general add constructively with parallactic motion, we do not attempt to account in detail for its minimal effect on our sensitivity.

We note that the upper limit of d₉ for which our Section 8 completeness curve calculation applies is dictated by lack of sufficient orbital drift for very distant, eccentric objects near aphelion. Because a given coadd epoch's reference image (Section 5.3) will in general be constructed using other coadd epochs with zero relative parallax, an orbiting object drifting sufficiently slowly will behave somewhat like a static source, and at least partially cancel itself out in our difference images. Assuming $e\leqslant 0.75$, the annual drift⁹ does not become smaller than the W1 FWHM until ${d}_{9}\gtrsim 2250\,\mathrm{au}$. Our pixel-level simulations suggest a negligible loss in sensitivity even in the case of 61 yr⁻¹ drift, and indeed the linear drift model of Section 7.3 recovers faint high proper motion stars among our list of transients down to $\mu \approx 0\buildrel{\prime\prime}\over{.} 2\,{\mathrm{yr}}^{1}$. To be conservative, we claim our completeness curve of Section 8 to be applicable only for ${d}_{9}\leqslant 2250\,\mathrm{au}$, though we have considerable sensitivity at larger distances in the very unlikely event that Planet Nine could be so far away and detectably luminous.

One might reasonably ask whether the same outlier rejection that removes satellite streaks and main-belt asteroids could also erase Planet Nine from our epochal coadds. We have performed injection tests to address this concern. We consider the worst-case (fastest-moving) scenario by injecting point sources with parallax of 825 day⁻¹ into individual exposures, and then running the modified exposures through our coaddition process. We inject 100 such fakes per W1 magnitude, for dozens of W1 magnitudes spanning $14.0\,\leqslant$ W1 ≤ 16.8. The outlier rejection should be most problematic for the brightest Planet Nine injections. This motivated our choice to extend these fakes as bright as W1 = 14.0, since the maximally W1-luminous Fortney et al. (2016) model would have W1 ≈ 14.2 at the inner edge of our search space (d₉ = 250 au).

We find that the outlier rejection causes a very small, but nonetheless measurable, reduction in the average fluxes we recover relative to the true fluxes of the injections. Near our 90% completeness threshold of W1 = 16.66 (Section 8), on average 7 mmag of flux is lost. The fraction of flux lost ramps up slowly for brighter injections. At the bright end, W1 = 14.0, 19 mmag of flux is lost on average.

Because the impact of unWISE outlier rejection on a maximally fast-moving Planet Nine near our 90% completeness threshold is less than a hundredth of a magnitude, we neglect this effect in our completeness analysis. In short, Planet Nine is expected to be too slow moving and faint to trigger unWISE outlier rejection at a level that significantly impacts our sensitivity.

The largest departures of WISE's pointing relative to solar elongation of 90° occur during so-called "Moon-avoidance maneuvers." These maneuvers coincide with the roughly monthly occasions on which the Moon crosses the WISE 90° elongation sight line. At the start of such a maneuver, WISE will gradually shift several degrees ahead in ecliptic longitude (λ) relative to 90° elongation, then avoid the Moon by jumping ∼5°–10° backward in λ. This backward jump in ecliptic longitude has the effect of lengthening WISE visits near in time to a Moon-avoidance maneuver to last up to ∼5–7 days. Coadds spanning such long periods of time would dramatically dilute the signal from an object orbiting at hundreds of astronomical units, which could undergo many dozens of arcseconds of parallax during such a timespan.

In order to avoid long coadd timespans that are strongly inconsistent with Equation (1), we enforce the following rules:

• 1.  
  No epochal coadd may include exposures during a Moon-avoidance maneuver. A list of start/end MJD values for the maneuvers relevant to our search region is provided in Table 2. By discarding the frames acquired during moon-avoidance maneuvers, we lose 1.6% of exposures.
• 2.  
  No epochal coadd may include exposures acquired both before and after the same Moon-avoidance maneuver.

This time slicing splits the five to seven day visits created by Moon-avoidance maneuvers into two short coadd epochs that span lengths of time consistent with or shorter than Equation (1). As a result, even though there are only seven seasons of W1 observations available, we end up with an average of 7.7 coadd epochs per coadd_id, thanks to the "extra" epochs created by Moon-avoidance maneuvers. Once the above rules are enforced, the MJD ranges spanned by exposures contributing to actual pixels in our epochal coadds closely follow Equation (1). The largest measured value of ${\tau }_{\mathrm{coadd}}$ for any pixel in our analysis is 1.65 days, which occurs at β = −523 (the expectation from Equation (1) at this $| \beta |$ is 1.64 days). The mean measured value of ${\tau }_{\mathrm{coadd}}$ across all pixels in our analysis is 1.01 days.

Note that WISE Moon-avoidance maneuvers are spaced such that a single coadd_id can have at most two coadd epochs during a single season. For cases in which one coadd_id has two coadd epochs in one season, these epochs are separated by a small handful of days.

For each time-resolved coadd identified by its unique (coadd_id, epoch) pair, we generate a corresponding reference stack. To do so, we combine all epochal coadds for the relevant coadd_id astrometric footprint, excluding the epoch under consideration. We also exclude any other epoch during the same season as the epoch for which the reference is being created. This latter constraint is important because, in cases where a single coadd_id has two epochs during the same season, these epochs are typically separated by only a couple of days, and Planet Nine's apparent motion could be slow enough that its detections at these two epochs have some overlap. To combine the epochal coadds into a reference image, we apply a weighted median filter on a per-pixel basis. The mean integer coverage of our reference stacks is 71.5 exposures, ∼6× the typical integer coverage of our epochal coadds.

We begin by astrometrically calibrating the reference coadd to 2MASS using SCAMP (Bertin 2006). We next use SCAMP to astrometrically calibrate the epochal coadd to 2MASS by comparing its source positions against the 2MASS-calibrated reference stack source positions. We then use SWARP (Bertin et al. 2002) to warp the reference and its uncertainty image to match the astrometry of the epochal coadd.

We use the warped reference to create two difference images for each epochal coadd. First, we perform a subtraction using the HOTPANTS package,¹⁰ which applies the Alard (2000) kernel-matching technique. Second, we perform a direct subtraction without kernel matching. We find that for WISE, which delivers highly stable space-based imaging, kernel matching does not play the same crucial role it does for ground-based images, where the PSF shape/size can vary dramatically. When creating the direct subtraction, we do take into account small zero-point differences between the epochal and remapped reference images, using the KSUM00 value obtained by HOTPANTS.

To extract transient sources from the subtracted images, we use Source Extractor (Bertin & Arnouts 1996) in "dual image mode," which makes use of static source properties in the remapped reference image to inform source detection in the subtracted image.

In addition to legitimate transients, the source extraction performed on our difference images also yields large numbers of false positives. We need to cull the list of transients for several reasons. First, the number of quadruplets that must be fit with Keplerian orbits during linking scales as the fourth power of transient number density, meaning that our full list of difference detections would be computationally expensive to link. Aside from computational limitations, we also expect the number of chance occasions on which spurious transients happen to link into Keplerian quadruplets to scale as the fourth power of the transient number density. Thus, eliminating as many spurious transients as possible from the outset allows us to obtain a manageable number of candidate quadruplets well fit by Keplerian orbits.

To remove false positives, we would ideally employ an approach like that of Goldstein et al. (2015), analyzing a variety of pixel-level features for each transient in order to classify it as legitimate versus spurious. Unfortunately, that methodology has thus far only been developed for the case of pointlike transients, whereas we are explicitly searching for a moving object which may appear slightly streaked. We therefore resort to a number of heuristics. These heuristics are designed to have minimal impact on our completeness, and have effects that can be easily incorporated into our simulations of Section 8. The following subsections detail each of the transient filtering heuristics.

6.3.1. 2MASS Sources

6.3.2. Bright Static W1 Sources

6.3.3. Kernel-matching Artifacts

6.3.4. Reference Image Brightness

We reject any transient with total flux less than the reference coadd brightness at its centroid. Given the W1 PSF and our coadd pixel size, this is roughly equivalent to rejecting any source with peak brightness fainter than the reference image at the same location by a factor of more than 7. For the transient flux, we adopt the Source Extractor flux_best measurement. It is not entirely straightforward to translate our reference brightness cut into a fraction of area masked, but we note that the transient magnitude histogram (Figure 5) is sharply peaked near W1 = 16.6, and for this magnitude 0.48% of the total search area is masked.

Zoom In Zoom Out Reset image size

6.3.5. Negative Fluxes

6.3.6. Ultra-bright Stars

6.3.7. Latents

6.3.8. Merging Per-coadd Transient Lists

6.3.9. Filtering Summary

After filtering, the final transient catalog contains 206,891 detections, corresponding to a mean density of 113 per square degree. Figure 5 shows a histogram of the transient W1 magnitudes, which peaks near W1 = 16.6. The number of transients per N_side = 128 HEALPix pixel (0.21 square degrees, Górski et al. 2005) is shown in Figure 6.

Zoom In Zoom Out Reset image size

The regions masked due to our various filters are not mutually exclusive, and in total we have masked 1.0% of the search area. This small loss of area is accounted for in the Monte Carlo analysis of Section 8.

The goal of our astrometric calibration is to ensure that the most appropriate and accurate spatial coordinates are input to our orbit linker (Section 7). The orbit linker computes model orbits in ICRS coordinates, and therefore, as described in Section 6, we have calibrated the astrometry of each epochal coadd to ICRS via 2MASS.

It is crucial to have a full accounting of the uncertainties on our astrometric measurements of transients, since these uncertainties will strongly influence the ${\chi }^{2}$ values of our best-fit orbits. We include three terms in our astrometric uncertainty model.

7.1.1. 2MASS Calibration Residuals

7.1.2. Systematic Proper Motion Shifts

7.1.3. Statistical Scatter

Statistical scatter due to background noise is by far the dominant source of astrometric uncertainty for our typical transient. We use the 300,000 injections described in Section 8 to fit an analytic relation to the (point source) per-coordinate statistical scatter as a function of W1 magnitude and integer coverage n_exp. We find that the results are well fit by

$$\begin{eqnarray}&&{\sigma }_{\mathrm{stat},\mathrm{PSF}}=0\buildrel{\prime\prime}\over{.} 41\times \displaystyle \frac{\sqrt{12/{n}_{\exp }}}{{({10}^{(16.0-{\rm{W}}1)/2.5})}^{0.69}}.\end{eqnarray} \tag{ 9 }$$

To be conservative, we wish to scale this uncertainty up so that it is appropriate even in the case of maximum intra-coadd parallax (825). Typically, for a fixed amount of signal, the astrometric uncertainty scales linearly with both FWHM and the effective amount of noise. The FWHM (n_eff) for 825 of parallax is 1.30× (1.48×) larger than that of a W1 point source, and therefore we adopt ${\sigma }_{\mathrm{stat}}$ = 1.30 × $\sqrt{1.48}$ × ${\sigma }_{\mathrm{stat},\mathrm{PSF}}$. Finally, if the positional error ellipse semimajor axis quoted by Source Extractor is larger than our model for ${\sigma }_{\mathrm{stat}}$, we assign ${\sigma }_{\mathrm{stat}}$ to the Source Extractor value.

7.1.4. Total Astrometric Uncertainty

When enumerating quadruplets, we do not simply allow all transients in the merged catalog to be paired with all other transients. We enforce two rules restricting which transients may be linked:

• 1.  
  No two transients may be linked if they are extracted from the same (coadd_id, epoch) coadd image.
• 2.  
  No two transients may be linked if they are separated temporally by <0.4 days. Because of the way that our coadd epochs are determined from WISE visit and Moon-avoidance maneuver boundaries, coadd pixels with mean MJD separated by less than 0.4 days will in general have been constructed from (at least partially) the same set of single-exposure images.

Otherwise all transients may be linked. Importantly, we allow linkages of transients that occur during the same season but are extracted from different coadd images, provided the second rule above is not violated. Note also that we impose no constraints on the relative fluxes of transients that may be linked.

Our basic orbit linking approach will be to enumerate all viable quadruplets of transients, and determine which can be well fit by a Keplerian orbit. We break up this process into two phases, a "prescreen" of a large number of possible candidates, followed by a detailed fit to Keplerian orbital parameters. The prescreen step starts with the complete list of n_transient transients, for each one identifying a region of the sky where that object might be at subsequent epochs. The actual shape of this window is set by our choice of a minimum orbital distance, which we take to be 250 au, and limits on the inclination. Our results of Section 9 place no restrictions on inclination. In this way we create n_transient lists of neighboring transients, with typically 40 neighbors each.

The next phase of our prescreen is to run through the neighbor lists, combining them to form unique quadruplets. We fit the observation times and sky positions of each quadruplet with a simple linear model that includes the magnitude of the parallax and the 2D plane-of-sky orbital drift; this model is appropriate for distant solar system bodies observed on timescales that are short compared with their orbital periods. If this model yields a ${\chi }^{2}\lt 80$, we accept the quadruplet; otherwise we reject it.

The second phase of our linking algorithm is to fit each prescreened quadruplet with Keplerian orbital parameters. For this step we use the Bernstein & Khushalani (2000) code, orbfit, designed for objects at Kuiper Belt distances and beyond. If orbfit successfully resolves a Keplerian orbit (see Section 7.4), we save the quadruplet. Finally we take the list of saved quadruplets and attempt to combine them into quintuplets, sextuplets, etc., if they have transients in common. In this way we start with ∼ $2\times {10}^{5}$ W1 transients, assess over 14 billion quadruplets, and perform Keplerian orbit fits to the ∼12,000 that pass prescreening. The final list contains just under 500 successful quadruplet linkages. We always retain every tuple with low orbfit ${\chi }^{2}$, never making cuts on the best-fit orbital elements or relative fluxes of constituent transients.

We must choose a goodness-of-fit threshold defining the boundary between quadruplets that are well fit versus poorly fit by a Keplerian orbit. We do so by choosing a ${\chi }^{2}$ threshold such that the expected false negative rate is 0.001. Our orbfit quadruplet (quintuplet) fits always have 2 or 3 (4 or 5) degrees of freedom (DOF). To obtain the desired level of false negatives, we therefore choose our ${\chi }^{2}$ threshold to be 13.8, 16.3, 18.5, and 20.5 for 2, 3, 4, and 5 DOF. The impact of our 0.001 failure rate in fitting orbits can be accounted for in combination with the simulations of Section 8, but it leads to a negligible shift in our 90% completeness threshold.

There are two types of completeness relevant to our search. The first completeness refers to the probability that a transient source of a given brightness is successfully detected in a reference-subtracted epochal coadd. This will be denoted f_det. Note that f_det will be a function of coadd integer coverage n_exp. The second completeness is the probability with which an object of a particular brightness is detected as a transient in $\geqslant 4$ epochal coadd difference images, given the available number of coadd epochs and the number of contributing exposures within each coadd epoch. Calculating ${f}_{\geqslant 4}$ is our ultimate goal, because we require that a Planet Nine candidate consist of $\geqslant 4$ transients which are linked along a low-${\chi }^{2}$ orbit. Calculating f_det as a function of integer coverage is one step towards determining ${f}_{\geqslant 4}$.

f_det is expected to approach unity for very bright transients and should decrease to zero for very faint transients, which will be indistinguishable from background noise. The situation is similar for ${f}_{\geqslant 4}$. The rejection of transients due to, e.g., our bright star masking means that ${f}_{\geqslant 4}$ will asymptote to a value slightly less than unity for bright objects.

For these reasons, we parameterize completeness curves, in both the case of f_det and ${f}_{\geqslant 4}$, with the model

$$\begin{eqnarray}&&f(W1)=\displaystyle \frac{A}{2}\times [1+{\rm{erf}}(-({\rm{W}}1-\mu )/\sigma )],\end{eqnarray} \tag{ 10 }$$

where W1 is the W1 magnitude, and there are three free parameters: A, μ, and σ. A is the asymptotic completeness value toward low W1, and is typically very close to unity. When fitting (A, μ, σ) we impose the constraint $A\leqslant 1$. μ is the W1 value at which the completeness is $A/2\approx 0.5$, the center of the transition from f ≈ 1 to f = 0. σ dictates the width of this transition.

In order to most accurately calculate ${f}_{\geqslant 4}$, it is necessary to compile a set of completeness curves ${f}_{\det }(W1;{n}_{\exp })$, for each value of n_exp. We have done so by injecting 300,000 pointlike fakes of varying W1 into our difference images. One hundred real difference images spread throughout our search footprint are used as the basis for 3000 distinct mock difference images, each with 100 fakes injected. We record whether each fake is recovered using a 1 pixel match radius. We then bin the fakes by n_exp, and fit Equation (10) to the fraction of fakes recovered as a function of W1 for each n_exp value.

Not all n_exp values are sufficiently sampled by these fakes, because many n_exp values represent just a tiny fraction of our search area. Most of the fakes have $6\leqslant {n}_{\exp }\leqslant 15$, and these are the n_exp values for which high-quality f_det curves can be created. Table 3 lists the best-fit completeness curve parameters (A, σ, μ) for these n_exp values.

Table 3.  Parameters Tabulated as a Function of n_exp, the Number of Exposures Contributing to a Coadd Pixel

nexp	$p({n}_{\exp })$	μ	σ	A	${\mu }_{\mathrm{pred}}$
3	0.0104	⋯	⋯	⋯	16.156
4	0.0108	⋯	⋯	⋯	16.312
5	0.0117	⋯	⋯	⋯	16.433
6	0.0139	16.512	0.302	0.9970	16.532
7	0.0197	16.624	0.356	0.9914	16.616
8	0.0338	16.689	0.368	0.9940	16.688
9	0.0657	16.766	0.347	0.9955	16.752
10	0.118	16.832	0.342	0.9940	16.809
11	0.175	16.867	0.352	0.9933	16.861
12	0.196	16.908	0.354	0.9945	16.908
13	0.162	16.948	0.355	0.9943	16.952
14	0.0986	16.997	0.360	0.9943	16.992
15	0.0450	17.044	0.368	0.9958	17.030
16	0.0193	⋯	⋯	⋯	17.065
17	9.81 × 10−3	⋯	⋯	⋯	17.098
18	4.70 × 10−3	⋯	⋯	⋯	17.129
19	2.29 × 10−3	⋯	⋯	⋯	17.158
20	1.27 × 10−3	⋯	⋯	⋯	17.186
21	7.12 × 10−4	⋯	⋯	⋯	17.212
22	5.42 × 10−4	⋯	⋯	⋯	17.238

Download table as:  ASCIITypeset image

One might hope that μ(n_exp) scales as expected from statistical noise proportional to $1/\sqrt{{n}_{\exp }}$:

$$\begin{eqnarray}&&{\mu }_{\mathrm{pred}}({n}_{\exp })=\mu ({n}_{\exp ,0})+2.5{{\rm{log}}}_{10}(\sqrt{{n}_{\exp }/{n}_{\exp ,0}}).\end{eqnarray} \tag{ 11 }$$

We take ${n}_{\exp ,0}=12$, given that 12 is both the median and mode of n_exp in our full set of epochal coadds. ${\mu }_{\mathrm{pred}}({n}_{\exp })$ calculated according to Equation (11) is also listed in Table 3, and it clearly agrees very well with the measured $\mu ({n}_{\exp })$ values for $6\leqslant {n}_{\exp }\leqslant 15.$

For our downstream Monte Carlo analysis, we employ ${f}_{\det }({\rm{W}}1;{n}_{\exp })$ completeness curves parameterized by Equation (10). For these parameterizations, we adopt μ = ${\mu }_{\mathrm{pred}}({n}_{\exp })$. Since σ and A appear roughly constant over the range of n_exp for which these have been measured, we adopt their median values (A = 0.9943, σ = 0.355) for all ${f}_{\det }(W1;{n}_{\exp })$ completeness curves during Monte Carlo simulations.

Our Monte Carlo simulation of ${f}_{\geqslant 4}$ requires as input statistics regarding the fraction of sky area in the full search footprint as a function of number of coadd epochs (n_epochs) and number of exposures per epochal coadd pixel (n_exp). Tables 4 and 3 list these values, referred to as $p({n}_{\mathrm{epochs}})$ and $p({n}_{\exp })$, respectively. Note that our analysis discards the small number of coadd pixels with integer coverage ≤2, as it is not possible to perform outlier rejection with fewer than three exposures.

Table 4.  Fraction of Pixels in Search Region as a Function of Number of Coadd Epochs

nepochs	$p({n}_{\mathrm{epochs}})$	nepochs	$p({n}_{\mathrm{epochs}})$
1	1.32 × 10−5	6	0.126
2	1.20 × 10−5	7	0.483
3	1.41 × 10−5	8	0.327
4	1.50 × 10−5	9	0.0635
5	2.50 × 10−4

Download table as:  ASCIITypeset image

We generate a large number of fake objects, each assigned its own magnitude W1. Each object is randomly assigned a number of coadded epochs using a multinomial draw from $p({n}_{\mathrm{epochs}})$. Each epoch of each object is then randomly assigned a corresponding number of exposures n_exp using a multinomial draw from $p({n}_{\exp })$. Given an object's assigned W1 and n_exp within each epoch, each appearance of that object will be detected with a probability given by ${f}_{\det }({\rm{W}}1;{n}_{\exp })$, which is specified by the parameters (μ, σ, A) described in the final paragraph of Section 8.3. Each epoch of each object is then marked as either "detected" or "undetected" by performing a Bernoulli trial with success probability given by ${f}_{\det }({\rm{W}}1;{n}_{\exp })$. ${f}_{\geqslant 4}({\rm{W}}1)$ is then calculated by computing the fraction of objects for which $\geqslant 4$ total detections are accumulated, using a sample of 10,000 objects per value of W1, with W1 = 16.20,16.21, ..., 17.59, 17.60.

Our measured ${f}_{\geqslant 4}$ completeness curve is plotted in Figure 7. Fitting the tabulated ${f}_{\geqslant 4}({\rm{W}}1)$ values with Equation (10) gives (μ, σ, A) = (16.89, 0.18, 0.998). This best-fit model reaches ${f}_{\geqslant 4}=0.90$ ("90% completeness") at W1 = 16.72. Requiring five rather than four detections results in a 0.11 magnitude reduction in sensitivity, with 90% completeness at W1 = 16.61. We can account for the masking of 1.0% of our search area (Section 6.3) by randomly zeroing out 1% of ${f}_{\det }(W1;{n}_{\exp })$ probabilities during the Monte Carlo simulation. We find that this pushes the 90% completeness threshold brighter by 6.4 (9.5) mmag in the case where four (five) linked detections are demanded.

Zoom In Zoom Out Reset image size

We can scale these measured point source completeness thresholds in a manner consistent with the final term of Equation (6) to account for variations in sensitivity as a function of $| \beta |$ and d₉. The resulting 90% completeness threshold is shown in Figure 8. Over the parameter space in which we are sensitive to the maximally W1-luminous Fortney et al. (2016) model, the mean 90% completeness threshold is 16.66 (16.54) when demanding at least four (five) linked detections. The ∼0.05 mag degradation in sensitivity relative to our point source measurement is a result of parallactic smearing, as the average ${\pi }_{9}$ value within our search space is 36.

Zoom In Zoom Out Reset image size

We caution that we might fail to link detections of Planet Nine if it saturates in W1 (${\rm{W}}1\lesssim 8.25$), although this would require a luminosity hundreds of times greater than any of the Fortney et al. (2016) models. We might also fail to detect Planet Nine if it is brighter than W1 = 9.5, due to the AllWISE catalog masking described in Section 6.3.2. Finally, we have not presented a detailed accounting of small-scale variations in our 90% completeness threshold, such as those introduced by bright star masking and the spatially varying number of coadd epochs shown in Figure 3.