The Atacama Cosmology Telescope: A Search for Planet 9

The most common way to discover solar system objects is to look for objects that have moved between two different exposures of the same patch of sky. This method is fast, but is limited by the depth of each image, since the object needs to be independently detected in both. This depth can be improved with longer exposures, but this is limited by the angular velocity of the object itself. Integrating longer than the time it takes the object to move by the size of the beam will just smear it out without any further gains in signal-to-noise ratio (S/N). This is the regime ACT is in for Planet 9.

For ACT sky coverage and sensitivity, it would take 3–4 yr of observations just to have a chance of detecting a Planet 9–like object that was not moving in the sky. By Kepler's laws, a planet with semimajor axis a, eccentricity e, and current solar distance r will have a Sun-centered angular speed of

$$\begin{eqnarray}&&v=1\buildrel{\,\prime}\over{.} 932\,{\mathrm{yr}}^{-1}\cdot \sqrt{\displaystyle \frac{a}{500\ \mathrm{au}}}{\left(\displaystyle \frac{r}{500\ \mathrm{au}}\right)}^{-2}\sqrt{1-{e}^{2}}.\end{eqnarray} \tag{ 1 }$$

At the same time, the Earth's orbit sweeps out a yearly parallax ellipse with a semimajor axis of

$$\begin{eqnarray}&&{\theta }_{\pi }=\displaystyle \frac{\mathrm{au}}{r}=6\buildrel{\,\prime}\over{.} 875\cdot \displaystyle \frac{500\,\mathrm{au}}{r}\end{eqnarray} \tag{ 2 }$$

corresponding to a maximum angular speed of

$$\begin{eqnarray}&&{v}_{\pi }=\displaystyle \frac{2\pi {\theta }_{\pi }}{\mathrm{yr}}=43\buildrel{\,\prime}\over{.} 20\,{\mathrm{yr}}^{-1}\cdot \displaystyle \frac{500\,\mathrm{au}}{r}.\end{eqnarray} \tag{ 3 }$$

For comparison, ACT has an angular resolution of 2.05'/1.40'/0.98' FWHM at f090/f150/f220, respectively. To avoid excessive smearing we need (v_π + v)Δt ≈ v_π Δt ≪ FWHM. For the smallest beam (f220) and a closest possible distance of ${r}_{\min }=300\,\mathrm{au}$, this gives us Δt ≪ 5 days. With 5 days of integration time and the current ACT survey strategy the expected Planet 9 S/N would be ∼1, more than 5 times too low for a detection, or more than 25 times too low in terms of observing time!

The smearing could be eliminated if one knew the orbit of the object one was looking for, since that would allow one to shift each exposure to track the object as it moves across the sky. In practice, while the Planet 9 hypothesis makes some predictions about its orbit, they are far too vague to allow for simple tracking like this. However, with enough computational resources it is possible to loop through every reasonable orbit, make a shifted stack of individual short exposures using that orbit, and then look for objects in the resulting image. This is the shift-and-stack algorithm, and has been used to successfully detect objects below the single-exposure sensitivity limit (Gladman et al. 1998; Bernstein et al. 2004; Holman et al. 2018).

Planet 9's orbit is characterized by its six orbital elements: semimajor axis a, eccentricity e, inclination i, longitude of ascending node Ω, argument of periapsis ω, and true anomaly ν. However, because of its large distance and corresponding slow motion, it is sufficient for us to consider its motion to be drifting linearly on the sky, modulated by parallax. This gives us the following five free parameters:

• 1.  
  The heliocentric R.A. (α) and decl. (δ) of the planet at a reference time t₀.
• 2.  
  The horizontal and vertical components of the heliocentric angular velocity v = [v_x , v_y ]. We define these in the local tangent plane, such that ${\alpha }_{\mathrm{obs}}=\alpha +{v}_{x}(t-{t}_{0})/\cos {\delta }_{0}$ and δ_obs = δ + v_y (t − t₀).
• 3.  
  The planet's current distance from the Sun, r, which we treat as constant in time.

With these, the shift-and-stack algorithm takes the following general form:

• 1.  
  Split the data into chunks with duration Δt, and make a sky map of each.
• 2.  
  For each reasonable value of r, v_x , and v_y , use these with the time t of each map to shift them according to their constant heliocentric angular velocity and parallactic motion, and stack them to produce a combined map.
• 3.  
  Use a filter matched to the noise and signal properties to look for point sources in each combined map.

We will go through the details of this process in the following sections.

4.3.1. The Sky Maps

ACT observes the sky by sweeping backwards and forwards in azimuth while the sky drifts past. As it does so, the temperature registered by the detectors is read out hundreds of times per second, forming a vector of time-ordered data d. We model d as

$$\begin{eqnarray}&&d={Pm}+n,\end{eqnarray} \tag{ 4 }$$

where m is the (beam-convolved, pixelated ⁴¹ ) sky in μK CMB temperature units, P is a response matrix that encodes the telescope's pointing as a function of time, and n is instrumental and atmospheric noise, which we model as Gaussian covariance N. The maximum-likelihood estimate for m given d is

$$\begin{eqnarray}\begin{array}{rcl}\hat{m} & = & {\left({P}^{T}{N}^{-1}P\right)}^{-1}{P}^{T}{N}^{-1}d\\ M & = & {\left({P}^{T}{N}^{-1}P\right)}^{-1}\end{array}\end{eqnarray} \tag{ 5 }$$

Here, M is the noise covariance matrix of the estimator $\hat{m}$.

4.3.2. The Matched Filter

To look for point sources in $\hat{m}$ we start by assuming that all sources are far enough apart that they can be considered in isolation. Our data model for a map containing a single point source in some pixel p is then

$$\begin{eqnarray}&&\hat{m}={{RQs}}_{p}+u\end{eqnarray} \tag{ 6 }$$

where s_p is the point-source flux density in pixel p in mJy at a reference frequency ν₀ = 150 GHz, and Q_i = δ_ip is a vector that is unity at the source location in pixel p and zero elsewhere. It takes us from just a single flux density value to a map with that value in a single pixel.

R = Bg(ν₀, ν)f(ν)A_p is a response matrix that takes us from that map to beam-convolved μK at the observed frequency. Here, B is the instrument beam normalized to have a pixel-space integral of one, g(ν₀, ν) is the conversion factor from flux density at the reference frequency ν₀ to the observed frequency ν, f(ν) is the conversion from flux density in mJy to beam-convolved peak height in μK, and A_p is pixel area in steradians. Since we expect Planet 9 to be a blackbody with temperature T ≈ 40 K, we have g(ν₀, ν) = b(ν, T)/b(ν₀, T), where

$$\begin{eqnarray}&&b(\nu ,T)=\displaystyle \frac{2h{\nu }^{3}}{{c}^{2}}\displaystyle \frac{1}{\exp \left(\tfrac{h\nu }{{k}_{B}T}\right)-1}\end{eqnarray} \tag{ 7 }$$

is the Planck law for surface brightness b; and ⁴²

$$\begin{eqnarray}\begin{array}{rcl}f(\nu ) & = & {\left(\displaystyle \frac{2{x}^{4}{k}_{B}^{3}{T}_{\mathrm{CMB}}^{2}}{{h}^{2}{c}^{2}}\displaystyle \frac{1}{4\sinh {\left(x/2\right)}^{2}}{10}^{23}{A}_{b}\right)}^{-1}\\ x & = & \displaystyle \frac{h\nu }{{k}_{B}{T}_{\mathrm{CMB}}}.\end{array}\end{eqnarray} \tag{ 8 }$$

Finally, u is the noise in $\hat{m}$ and has a covariance matrix U. For the purposes of point-source detection, u consists of everything in $\hat{m}$ that is not the point source, which includes both the instrumental and atmospheric noise described by the M covariance matrix from before, but also the CMB, Cosmic Infrared Background (CIB), Galactic dust, etc.

Given this model for $\hat{m}$, the maximum-likelihood estimate for the point-source flux density at the reference frequency is

$$\begin{eqnarray}\begin{array}{rcl}{\hat{s}}_{p} & = & {\rho }_{p}/{\kappa }_{p},{\rho }_{p}={Q}^{T}{R}^{T}{U}^{-1}\hat{m}\\ {\sigma }_{{\hat{s}}_{p}} & = & 1/\sqrt{{\kappa }_{p}},{\kappa }_{p}={Q}^{T}{R}^{T}{U}^{-1}{RQ}\end{array}\end{eqnarray} \tag{ 9 }$$

Here, ${\sigma }_{{\hat{s}}_{p}}$ is the standard deviation of ${\hat{s}}_{p}$, κ_p is the corresponding inverse variance, and ρ_p is the inverse variance weighted flux density.

So far we have only estimated the point-source flux density in some pixel p. But since we do not a priori know where on the sky the planet could be, we need to estimate the flux density in every pixel, resulting in the flux density sky map $\hat{s}$ and corresponding uncertainty ${\sigma }_{\hat{s}}$ given by: ⁴³

$$\begin{eqnarray}\begin{array}{rcl}\hat{s} & = & \rho /\kappa ,\rho ={R}^{T}{U}^{-1}\hat{m}\\ {\sigma }_{\hat{s}} & = & 1/\sqrt{\kappa },\kappa =\mathrm{diag}({R}^{T}{U}^{-1}R)\end{array}\end{eqnarray} \tag{ 10 }$$

where the division is done pixel by pixel. The corresponding S/N is

$$\begin{eqnarray}&&{\rm{S}}/{\rm{N}}={\boldsymbol{s}}/{\sigma }_{{\boldsymbol{s}}}=\rho /\sqrt{\kappa }=\frac{{R}^{T}{U}^{-1}\hat{m}}{\sqrt{\mathrm{diag}({R}^{T}{U}^{-1}R)}}\end{eqnarray} \tag{ 11 }$$

which we recognize as the matched filter for $\hat{m}$. This S/N map is what one would usually use for object detection, e.g., by identifying peaks with S/N > 5. As we shall see in Section 4.7 the shift-and-stack parameter search complicates this, but the general idea stays the same.

4.3.3. Stacking

If we have multiple estimates { s _i } built from independent chunks of data, such as the few-day chunks we will use in the shift-and-stack algorithm, these combine straightforwardly: ⁴⁴

$$\begin{eqnarray}\begin{array}{rcl}{\rho }_{\mathrm{tot}} & = & \sum _{i}{\rho }_{i},{\kappa }_{\mathrm{tot}}=\sum _{i}{\kappa }_{i}\\ {{\boldsymbol{s}}}_{\mathrm{tot}} & = & {\rho }_{\mathrm{tot}}/{\kappa }_{\mathrm{tot}},{\sigma }_{{{\boldsymbol{s}}}_{\mathrm{tot}}}=1/\sqrt{{\kappa }_{\mathrm{tot}}}.\end{array}\end{eqnarray} \tag{ 12 }$$

Sadly, the presence of the same CMB, CIB, etc. in each chunk of data breaks the assumption of independence that this expression builds on. It would be possible to build a more complicated expression that takes this into account, but given the computationally expensive parameter search we perform we need the stacking operation to be as fast and simple as possible. Thankfully we can eliminate these correlated components by simply subtracting the time-averaged mean of the sky from each chunk of data.

4.3.4. Mean Sky Subtraction

We can avoid the complications of the CMB, CIB, etc. acting as correlated noise common to the data chunks by subtracting a high S/N estimate of the mean sky from each chunk of data before mapping it. This eliminates any static part of the sky such as the CMB, CIB, Galactic emission, etc. (including any we do not know about), and leaves only time-dependent signals such as the planet we are looking for, as well as variable point sources (which can be masked) and transients (which are rare enough that we can ignore them). The cost is a small increase in the noise if the mean sky model is not noise free, and a partial subtraction of the signal itself that must be estimated and corrected for. For this search we use the ACT+Planck combined maps described in Naess et al. (2020), but extended to include the 2019 season of data.

Aside from letting us stack using Equation (12), mean sky subtraction has the effect of removing all but the instrumental and atmospheric noise from the individual sky maps, and hence the matched filter noise covariance matrix U reduces to M. Inserting this into Equation (10) we get:

$$\begin{eqnarray}&&\rho ={R}^{T}{\overbrace{{P}^{T}{N}^{-1}d}}^{\mathrm{rhs}}\end{eqnarray} \tag{ 13 }$$

The part labeled "rhs" is a map that is much cheaper to compute than $\hat{m}$ because it avoids the expensive inversion ${({P}^{T}{N}^{-1}P)}^{-1}=M$ that must usually be done using iterative methods like conjugate gradients. ⁴⁵ That leaves us with κ, which we approximate as

$$\begin{eqnarray}&&{\kappa }_{i}={R}_{{ji}}{M}_{{jk}}^{-1}{R}_{{ki}}\approx \alpha {R}_{{ji}}^{2}{w}_{j}\end{eqnarray} \tag{ 14 }$$

where the map w is an approximation pixel-diagonal of M⁻¹ built assuming white (uncorrelated) noise and α is a factor that compensates for the mean error we make by replacing M⁻¹ with $\vec{w}$. We determine α by evaluating a few pixels of the exact κ.

4.3.5. Ad-hoc Filter

4.3.6. Point-source Handling

4.3.7. Dust Masking

The distance to and velocity of Planet 9 are relatively poorly determined, but we can infer rough limits on the acceptable fit from Figure 15 of B19, as shown in Table 2. We see that Planet 9's current distance is limited to 300 au ≲ r ≲ 1300 au. Equation (1) for the heliocentric angular velocity of the planet can be re-expressed as

$$\begin{eqnarray}\begin{array}{rcl}v & = & {v}_{\mathrm{ref}}{\left(\displaystyle \frac{r}{500\,\mathrm{au}}\right)}^{-2}\\ {v}_{\mathrm{ref}} & = & 1\buildrel{\,\prime}\over{.} 93\,{\mathrm{yr}}^{-1}\cdot \sqrt{a(1-{e}^{2})/(500\ \mathrm{au})}\end{array}\end{eqnarray} \tag{ 15 }$$

and from Table 2 we see that v_ref is in the range 16 to 23 yr⁻¹ for all of the acceptable fits. ⁴⁶ This means that only a hollow cone in our r, v_x , v_y parameter space needs to be explored.

While in theory there is a continuum of possible parameter values inside this cone, in practice the limited angular resolution of the telescope means that very similar parameters are indistinguishable. From Equation (2) we see that getting the distance wrong by δ r results in a parallax ellipse that is bigger by

$$\begin{eqnarray}&&\delta {\theta }_{r}\approx -0\buildrel{\,\prime}\over{.} 14\cdot \displaystyle \frac{\delta r}{10\,\mathrm{au}}{\left(\displaystyle \frac{500\,\mathrm{au}}{r}\right)}^{2}.\end{eqnarray} \tag{ 16 }$$

Thus shift-stacking with the wrong distance leaves a residual ellipse with a radius of ∣δ θ_r ∣. If we step through distances in steps of Δr, then δ r will take on values in the range $[-\tfrac{{\rm{\Delta }}r}{2},\tfrac{{\rm{\Delta }}r}{2}]$. Using Equations (16) and (A4) from Appendix A, we see that on average, this increases the beam FWHM in quadrature by:

$$\begin{eqnarray}&&{\rm{\Delta }}{\mathrm{FWHM}}_{{\rm{r}}}=0\buildrel{\,\prime}\over{.} 093\cdot \displaystyle \frac{{\rm{\Delta }}r}{10\,\mathrm{au}}{\left(\displaystyle \frac{500\,\mathrm{au}}{r}\right)}^{2}.\end{eqnarray} \tag{ 17 }$$

Similarly, getting the speed wrong by δ v will over a time T = t − t₀ accumulate to a position error of

$$\begin{eqnarray}&&\delta {\theta }_{v}=0\buildrel{\,\prime}\over{.} 3\cdot \displaystyle \frac{\delta v}{0\buildrel{\,\prime}\over{.} 1\,{\mathrm{yr}}^{-1}}\displaystyle \frac{T}{3\,\mathrm{yr}}.\end{eqnarray} \tag{ 18 }$$

For a velocity step of Δv we get, using Equation (A6),

$$\begin{eqnarray}&&{\rm{\Delta }}{\mathrm{FWHM}}_{v}=0\buildrel{\,\prime}\over{.} 29\cdot \displaystyle \frac{{\rm{\Delta }}v}{0\buildrel{\,\prime}\over{.} 1\,{\mathrm{yr}}^{-1}}\displaystyle \frac{T}{3\,\mathrm{yr}}\end{eqnarray} \tag{ 19 }$$

where we have included a factor $\sqrt{2}$ in the numerical factor to take into account the smearing in both the x and y directions. The factor T depends on when in the ACT observing campaign each observation was taken, but will at most be three years if we choose t₀ to be the midpoint of ACT observations. The integration time Δt also results in smearing,

$$\begin{eqnarray}&&{\rm{\Delta }}{\mathrm{FWHM}}_{t}=0\buildrel{\,\prime}\over{.} 080\cdot \displaystyle \frac{500\mathrm{au}}{r}\displaystyle \frac{{\rm{\Delta }}t}{\mathrm{day}}\end{eqnarray} \tag{ 20 }$$

as does the pixel window

$$\begin{eqnarray}&&{\rm{\Delta }}{\mathrm{FWHM}}_{\mathrm{pix}}=0\buildrel{\,\prime}\over{.} 68\cdot \displaystyle \frac{\mathrm{res}}{1^{\prime} }\end{eqnarray} \tag{ 21 }$$

where res is the pixel side length. ^{47, 48} Together these effects make up our smearing budget, and each must be chosen small enough that their combined effect does not overly degrade the S/N. We choose

• 1.  
  ${\rm{\Delta }}v=0\buildrel{\,\prime}\over{.} 1\,\mathrm{yr}\ \Rightarrow {\rm{\Delta }}{\mathrm{FWHM}}_{v}=0\buildrel{\,\prime}\over{.} 29$ for T = 3 yr.
• 2.  
  ${\rm{\Delta }}t=3\ \mathrm{days}\ \Rightarrow {\rm{\Delta }}{\mathrm{FWHM}}_{t}=0\buildrel{\,\prime}\over{.} 40$ for r = 300 au, which is the closest distance we will consider.
• 3.  
  ${\rm{\Delta }}r=33\,\mathrm{au}\cdot \left(\tfrac{500\,\mathrm{au}}{r}\right)\Rightarrow {\rm{\Delta }}{\mathrm{FWHM}}_{r}=0\buildrel{\,\prime}\over{.} 31$. This results in the discrete set of distances: 300, 321, 346, 375, 409, 450, 500, 563, 643, 750, 900, 1125, 1500, and 2000 au. The last two distance bins are more distant than Planet 9 is likely to be, but are included because of their low computational cost.
• 4.  
  $\mathrm{res}=1^{\prime} \Rightarrow {\rm{\Delta }}{\mathrm{FWHM}}_{\mathrm{pix}}=0\buildrel{\,\prime}\over{.} 68$. That is, we use a pixel size of 1'. ⁴⁹

These combine in quadrature to ${\rm{\Delta }}\mathrm{FWHM}=0\buildrel{\,\prime}\over{.} 90$, which when combined with our beams represents a 9/19/34% increase in beam size and loss in S/N in the f090/f150/f220 bands respectively. The largest contribution to this is the 1' pixel size. With a $0\buildrel{\,\prime}\over{.} 5$ pixel size these numbers would instead have been 5/11/21% at a cost of 4×as high CPU and memory budgets. We might consider using smaller pixels when we revisit this in the future.

The full, quantized search space is visualized in Figure 2. In total the r, v_x , v_y parameter space has 25,837 cells.

Zoom In Zoom Out Reset image size

After splitting the 2013–2019 ACT data set into 3 day chunks and building matched filter maps for each, we were left with 3834 pairs of ρ and κ maps taking up a total of 1.9 TB of disk space. Since these maps are in units of equivalent flux density at the reference frequency ν₀ = 150 GHz, maps from different arrays and bandpasses that were observed at the same time, and hence all have the same shifts, can be directly combined before the main shift-and-stack search. This resulted in a more manageable 787 pairs taking up 220 GB.

The analysis was performed in 10° × 10° tiles with an additional 1° padding on all sides using data "belonging" to neighboring tiles to avoid discontinuities at tile edges. For each tile we loop over our parameter space and keep track of the highest-S/N value of v_x and v_y in each pixel for each value of r. This is illustrated in the pseudo-code below:

Here, r takes on the values 300, 321, 346, 375, 409, 450, 500, 563, 643, 750, 900, 1125, 1500, and 2000 au. For each value we visit all velocities v_x = iΔv, v_y = jΔv, where i and j are integers and

$$\begin{eqnarray}&&\begin{array}{l}1\buildrel{\,\prime}\over{.} 6\,{\mathrm{yr}}^{-1}\cdot {\left(\displaystyle \frac{r}{500\,\mathrm{au}}\right)}^{-2}-{\rm{\Delta }}v\leqslant \sqrt{{v}_{x}^{2}+{v}_{y}^{2}}\\ \quad \leqslant 2\buildrel{\,\prime}\over{.} 3\,{\mathrm{yr}}^{-1}\cdot {\left(\displaystyle \frac{r}{500\,\mathrm{au}}\right)}^{-2}.\end{array}\end{eqnarray} \tag{ 22 }$$

The function shift applies the coordinate transformation from observed coordinates at time t = t₀ + T to heliocentric coordinates at time t₀, taking into account both parallax for the distance r and the planet's angular velocity v_x , v_y . We use bilinear interpolation to allow for fractional pixel shifts. This function is the most time-critical part of the search, so it was implemented in optimized C using AVX intrinsics and OpenMP parallelization. Since the distance and direction each pixel is displaced changes slowly as a function of position in the map, we use the same displacement for blocks of 8 × 8 pixels, saving a large number of trigonometric operations at no loss of S/N. Overall our implementation is 480 times faster than a straightforward numpy/scipy implementation.

The function update updates result to maintain a running record of the highest S/N observed in each pixel, and what value of ρ_tot, κ_tot, v_x , and v_y that occurred for. We maintain one such result for each value of r because both bias from mean sky subtraction and the appropriate S/N threshold for a detection (which depends on the effective number of trials) depend on r.

Mean sky subtraction mainly removes the static parts of the sky, but it also subtracts some of the signal from moving objects. These appear as a smeared-out tracks in the mean sky map, and since part of an object's track necessarily overlaps with its position in each individual exposure, mean sky subtraction will always lead to a loss of signal power. The size of the bias is both distance dependent (because more distant objects move less and hence overlap more with the mean sky) and position dependent (because areas with less coverage will see less of the object's motion).

To characterize this bias we performed a set of signal injection and recovery simulations. We created a catalog of fake planets with the same flux density forming an equispaced grid in heliocentric R.A. and decl. at t = t₀, with values of r stepping through the 14 values we consider in the parameter search for every 14 grid positions in R.A., and v taking on the corresponding 14 values 180, 157, 135, 115 , 097, 080, 065 , 051, 039, 029, 020, 013, 007, and 004 yr⁻¹. An example of such a grid of simulated objects is shown in panel A of Figure 3.

Given this catalog, we compute the Earth-centered position of each object for each 3 day map, and use the instrument beam, bandcenter, and map-maker filter to simulate a signal-only image ${m}_{\mathrm{sim}}$ of these sources in the same units as our data maps: CMB temperature increment in μK. This is shown in panel B of the figure.

These were used to build new matched filter maps ${\rho }_{\mathrm{sim}}^{\mathrm{raw}}={R}^{T}w\circ {m}_{\mathrm{sim}}$, where ◦ is the element-wise product and w is the white noise inverse variance map from Section 4.3.4. ⁵⁰ Using $w\circ {m}_{\mathrm{sim}}$ instead of Equation (13) is an approximation, but based on a small number of full time-domain simulations it appears to be accurate to <5%.

To capture the effect of mean sky subtraction we define the mean flux density map for the simulations

$$\begin{eqnarray}&&{F}_{\mathrm{mean}}={\rho }_{\mathrm{sim},\mathrm{tot}}/{\kappa }_{\mathrm{tot}},{\rho }_{\mathrm{sim},\mathrm{tot}}=\sum _{i}{\rho }_{\mathrm{sim},i}^{\mathrm{raw}}\end{eqnarray} \tag{ 23 }$$

where i loops over all of the individual maps and the division is element-wise. Panels C and D of Figure 3 shows examples of individual and mean matched filter maps, but differ from the ones described here by not being noise free (see Section 4.9). Using these maps we define the mean sky subtracted simulations

$$\begin{eqnarray}&&{\rho }_{\mathrm{sim}}={\rho }_{\mathrm{sim}}^{\mathrm{raw}}-\kappa {F}_{\mathrm{mean}},\end{eqnarray} \tag{ 24 }$$

This was done individually for each bandpass, both to avoid mixing maps with different beams and to reflect what was done to the actual data.

Finally, we ran the shift-and-stack procedure from Section 4.5 on the simulated data set, and read off the recovered flux density for each simulated source. We find practically no dependence on the direction of the velocity, and therefore average the data points for different velocities for our final bias model, resulting in bias maps b (r) with a resolution of 7° × 4°. These are shown for the closest and furthest Planet 9 distance considered in Figure 4. The bias changes smoothly with position and is well resolved even with these large pixels. The standard deviation of the data points going into each pixel is about 0.5%, which we take as the uncertainty on our bias maps. We use this to define ${\rho }_{\mathrm{tot}}^{\mathrm{debiased}}={\rho }_{\mathrm{tot}}\circ {\boldsymbol{b}}$ and ${\kappa }_{\mathrm{tot}}^{\mathrm{debiased}}={\kappa }_{\mathrm{tot}}\circ {{\boldsymbol{b}}}^{2}$, from which bias-free flux densities can be recovered via Equation (10).

Our search method results in a map for each r where each pixel has the maximum S/N across all of the velocity parameters for that r. To construct a list of detection candidates and detection limit maps we need to know the background distribution of these S/N values. This is made difficult by the varying depth, and varying temporal and spatial distribution of the data used in the search. The effective number of trials is a strong function of r, and the individual trials are correlated, with the correlation depending on how densely the ACT observations covers each spot of the map. The S/N distribution should therefore vary both as a function of r and position.

The simple approach of multiplying the number of beams in the map (∼30 million) with the total trial number (25,837) to get a total number of trials (∼10¹²) and a corresponding Gaussian quantile (7σ) does not work. Aside from overestimating the effective number of trials, it would also lead to the search grossly preferring candidates with low r by not penalizing the much larger parameter space for low r compared to high r.

Instead we will take the approach of transforming S/N into an overall detection statistic z that follows a simple, uniform Gaussian distribution, at least for its high-z tail. This procedure is described in Appendix B, where we find that

$$\begin{eqnarray}&&z=\xi ({\rm{S}}/{\rm{N}})=({\rm{S}}/{\rm{N}}-{\mu }_{z})/{\sigma }_{z}\end{eqnarray} \tag{ 25 }$$

where μ_z and σ_z are functions of distance r and position in the map.

With the normalized detection statistic z in hand, we build a set of preliminary candidate detections by selecting peaks with z > 3.5. Given the large sky area covered, this low threshold will result in a large number of candidates, the vast majority of which would of course simply be noise fluctuations (especially considering that we expect at most one real object), but that allows us to get a good handle on the background distribution that any real objects would stand out from.

To better understand the background, we took advantage of the fact that the planet signature would be positive in our maps and repeated the whole search with the sign of all the data flipped. No signal is expected in the sign-flipped search, but it shares the same noise properties and many of the systematics (e.g., variable point sources and edge artifacts), so it gives a good estimate of the background detection rate.

We classified each candidate as Planet 9–like or general based on whether they satisfied the expected bounds on Planet 9's orbital inclination, 10° < i < 30° (B19). ⁵¹ The inclination is not one of the free parameters of our fit, but we can approximate this selection by transforming the candidate coordinates and velocities into ecliptic coordinates, and requiring

$$\begin{eqnarray}&&10^\circ \lt \hat{i}\lt 30^\circ \,,\,\hat{i}=\sqrt{{\beta }^{2}+{v}_{\beta }^{2}/{v}_{\lambda }^{2}}\end{eqnarray} \tag{ 26 }$$

and where λ, β are the ecliptic longitude and latitude, respectively, and v_λ , v_β are the velocity components in those directions. ⁵² The formula for $\hat{i}$ assumes that orbits have $\beta (\lambda )=i\sin (\lambda -{\lambda }_{0})$, which is a decent approximation as long as i is small.

Finally, the top 100 candidates from each list were visually inspected using the diagnostic maps below, and candidates with obvious problems were cut. ⁵³

• 1.  
  Maps of the detection statistic z, to ensure that there is a sensible local maximum at the candidate location. The edges of the mask sometimes cause extended artifacts that are easily identified this way.
• 2.  
  The shift-stacked maps using the best-fit orbital parameters for each candidate, again to ensure a compact, isolated maximum.
• 3.  
  Raw sky maps for the f090 and f150 bands at the candidate location, to ensure that there is not an unmasked point source or dust structure there.
• 4.  
  Individual 3 day matched filter maps, which were useful for disqualifying false positives caused by variable point sources and transients.

With the exception of the 3 day maps, these are shown as the image columns in Tables 3 and 4.

Table 3. Top 10 Planet 9–Like Candidates, Sorted by the Detection Statistic z (see Section 4.7 or Appendix B for Definition)

#	z map	Stack	f090	f150	R.A.	Decl.	z	F	Δ F	r	v x	v y
					(deg)	(deg)		(mJy)	(mJy)	(au)	(' yr−1)	(' yr−1)
1					−167.54	1.04	5.17	8.3	1.8	375	2.2	−2.9
2					−50.84	−9.16	5.05	11.5	2.4	375	0.2	3.0
3					−70.32	0.34	5.00	14.8	3.1	321	0.6	4.5
4					−150.37	−4.86	5.00	23.2	5.5	643	−0.1	−1.1
5					−179.17	−0.23	4.98	8.5	1.9	1125	0.0	0.4
6					179.01	4.34	4.92	6.3	1.4	500	1.3	−1.8
7					−173.55	15.20	4.92	4.1	0.9	346	−0.7	4.1
8					5.25	−0.70	4.87	5.6	1.3	643	−0.1	1.3
9					52.66	−2.60	4.87	10.1	2.3	563	−0.2	−1.7
10					−42.35	−45.80	4.87	8.4	1.7	500	0.3	1.8

Note. The columns are: #: the rank in terms of peak z value. z map: a thumbnail of the z map centered on the candidate. Stack: the shift-and-stack (i.e., motion-corrected) map for the best-fit parameters. f090/f150: filtered versions of the mean sky model in the f090/f150 band. Because these do not include any motion correction, no Planet 9 signal is expected here, but they are useful for seeing how "clean" each candidate's neighborhood is, e.g., if there are any bright point sources, dust clumps, or map edges at or near the candidate's location. All thumbnails are $45^{\prime} \times 45^{\prime}$ centered on the candidates. R.A., decl.: candidate's J2000 heliocentric equatorial coordinates on modified Julian day (MJD) 57688. z: the candidate's detection statistic z. F , ΔF : flux in the f150 band in mJy, assuming a 40 K blackbody, and its uncertainty. r : distance from the Sun, in au. v _x , v _y : intrinsic motion in arcminutes per year.

Download table as:  ASCIITypeset image

Table 4. Like Table 3, but for the General Candidates

#	z map	Stack	f090	f150	R.A.	Decl.	z	F	Δ F	r	v x	v y
					(deg)	(deg)		(mJy)	(mJy)	(au)	(' yr−1)	(' yr−1)
1					−162.40	12.65	5.65	4.4	0.8	300	0.7	5.9
2					94.55	−29.48	5.64	9.7	1.8	500	1.6	1.5
3					116.58	−46.50	5.60	25.1	4.3	300	4.3	4.3
4					36.91	−12.81	5.51	13.7	3.4	1500	0.0	−0.1
5					59.20	1.52	5.48	11.3	2.4	643	0.6	0.6
6					69.03	−21.10	5.40	9.0	1.7	300	3.9	1.4
7					179.90	13.94	5.28	4.8	0.9	346	−3.7	−2.3
8					−8.19	−17.78	5.28	12.1	2.8	1125	−0.3	0.0
9					−69.90	−19.31	5.15	14.9	4.4	2000	0.0	0.1
10					−102.24	13.04	5.14	6.5	1.4	500	−1.4	−1.6

Download table as:  ASCIITypeset image

It is useful to be able to translate the survey depth into detection limit maps. To do this, we need the false negative rate as a function of the detection statistic z. We found this by repeating the signal injection, search, and detection procedure from Section 4.6 with two important differences:

• 1.  
  Simulated sources were added to the data instead of replacing it, resulting in noisy simulations.
• 2.  
  For each simulated source we chose a target z ∈ {4.50, 4.75, 5.00, 5.25, 5.50, 5.75, 6.00, 6.50, 7.00, 8.00, 10.00, 15.00}, translated this into an S/N ratio using ξ⁻¹(z); (see Section 4.7 and Appendix B), and combined it with the local survey depth to define a simulated flux density$s={\xi }^{-1}(z)/\sqrt{{\kappa }_{\mathrm{tot}}^{\mathrm{debiased}}}$.

We then ran the standard mean sky subtraction and candidate search on the maps, and computed the fraction of the injected sources that were ultimately recovered as a function of z. The result is plotted as the curve "detection chance" in Figure 5. Overall we find that a source bright enough to correspond to z = 5.3 has a 95% chance of being detected. Hence, the 95% flux density detection limit map is given by

$$\begin{eqnarray}&&{s}_{\mathrm{lim}}^{95 \% }={\xi }^{-1}(5.3)/\sqrt{{\kappa }_{\mathrm{tot}}^{\mathrm{debiased}}}.\end{eqnarray} \tag{ 27 }$$

Aside from its position dependence this limit is also distance dependent, since both ξ and ${\kappa }_{\mathrm{tot}}^{\mathrm{debiased}}$ depend on r.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Given a model for Planet 9's luminosity we can translate the flux density limit to a distance limit. Since the flux density falls with the square of the distance, the distance limit ${r}_{\mathrm{lim}}^{95 \% }$ can be found as the solution to the equation

$$\begin{eqnarray}&&{({r}_{\mathrm{ref}}/r)}^{2}{s}_{\mathrm{ref}}={s}_{\mathrm{lim}}^{95 \% }({r}_{\mathrm{lim}})\end{eqnarray} \tag{ 28 }$$

with Table 1 showing examples of the reference flux s_ref for r_ref = 500 au for different Planet 9 scenarios. ⁵⁴

Consider a Gaussian beam with standard deviation σ, such that its profile is

$$\begin{eqnarray}&&b(r)={e}^{-\tfrac{1}{2}\tfrac{{r}^{2}}{{\sigma }^{2}}}={e}^{-\tfrac{1}{2}\tfrac{{x}^{2}+{y}^{2}}{{\sigma }^{2}}}\end{eqnarray} \tag{ A1 }$$

where $r=\sqrt{{x}^{2}+{y}^{2}}$. Partially uncorrected parallax smears this beam along an ellipse with some semimajor axis μ. The simplest and worst case of this is smearing along a circle with radius μ, so that is what we will consider here. This results in the smeared beam

$$\begin{eqnarray}&&\begin{array}{l}{b}_{2}(r,\mu )=\displaystyle \frac{1}{2\pi }{\displaystyle \int }_{0}^{2\pi }d\theta {e}^{-\tfrac{1}{2}\tfrac{{\left(x-\mu \cos \theta \right)}^{2}+{\left(y-\mu \sin \theta \right)}^{2}}{{\sigma }^{2}}}\\ \quad \propto b(r){\displaystyle \int }_{0}^{2\pi }d\theta {e}^{\tfrac{x\mu \cos \theta +y\mu \sin \theta }{{\sigma }^{2}}}\\ \quad \approx b(r){\displaystyle \int }_{0}^{2\pi }d\theta \left(1+\displaystyle \frac{x\mu \cos \theta +y\mu \sin \theta }{{\sigma }^{2}}\right.\\ \quad +\left.\displaystyle \frac{1}{2}{\left(\displaystyle \frac{x\mu \cos \theta +y\mu \sin \theta }{{\sigma }^{2}}\right)}^{2}\right)\\ \quad \propto b(r)\left(1+\displaystyle \frac{1}{4}\displaystyle \frac{{r}^{2}{\mu }^{2}}{{\sigma }^{4}}\right)\end{array}\end{eqnarray} \tag{ A2 }$$

where we have assumed μ ≪ σ and have ignored any factors that just scale the overall amplitude of the function. If all values $\mu \in [-\tfrac{{\rm{\Delta }}}{2},\tfrac{{\rm{\Delta }}}{2}]$ occur with equal weight, then the average beam across all of these will be:

$$\begin{eqnarray}&&\begin{array}{l}{b}_{\mathrm{eff}}(r)=\displaystyle \frac{1}{{\rm{\Delta }}}{\displaystyle \int }_{-\tfrac{{\rm{\Delta }}}{2}}^{\tfrac{{\rm{\Delta }}}{2}}d\mu {b}_{2}(r,\mu )=\displaystyle \frac{1}{{\rm{\Delta }}}b(r){\displaystyle \int }_{-\tfrac{{\rm{\Delta }}}{2}}^{\tfrac{{\rm{\Delta }}}{2}}d\mu \\ \quad \times \left(1+\displaystyle \frac{1}{4}\displaystyle \frac{{r}^{2}{\mu }^{2}}{{\sigma }^{4}}\right)=b(r)\left(1+\displaystyle \frac{1}{24}\displaystyle \frac{{r}^{2}{{\rm{\Delta }}}^{2}}{{\sigma }^{4}}\right)\\ \quad \approx {e}^{-\tfrac{1}{2}\tfrac{{r}^{2}}{{\sigma }^{2}}}{e}^{\tfrac{1}{24}\tfrac{{r}^{2}{{\rm{\Delta }}}^{2}}{{\sigma }^{4}}}={e}^{-\tfrac{1}{2}\tfrac{{r}^{2}}{{\sigma }_{\mathrm{eff}}^{2}}}\end{array}\end{eqnarray} \tag{ A3 }$$

where

$$\begin{eqnarray}&&\begin{array}{l}{\sigma }_{\mathrm{eff}}^{2}=\displaystyle \frac{{\sigma }^{2}}{1-\tfrac{1}{12}\tfrac{{{\rm{\Delta }}}^{2}}{{\sigma }^{2}}}\approx {\sigma }^{2}+\displaystyle \frac{1}{12}{{\rm{\Delta }}}^{2}\\ \quad \iff \,{\mathrm{FWHM}}_{\mathrm{eff}}^{2}={\mathrm{FWHM}}^{2}+\displaystyle \frac{2\mathrm{log}(2)}{3}{{\rm{\Delta }}}^{2}.\end{array}\end{eqnarray} \tag{ A4 }$$

So circular smearing adds in quadrature to the beam size.

We here smear the beam linearly in the x direction with $\mu \in [-\tfrac{{\rm{\Delta }}}{2},\tfrac{{\rm{\Delta }}}{2}]$.

$$\begin{eqnarray}&&\begin{array}{l}{b}_{\mathrm{eff}}=\displaystyle \frac{1}{{\rm{\Delta }}}{\displaystyle \int }_{-\tfrac{{\rm{\Delta }}}{2}}^{\displaystyle \frac{{\rm{\Delta }}}{2}}d{\mu }^{-\tfrac{1}{2}\tfrac{{\left(x-\mu \right)}^{2}+{y}^{2}}{{\sigma }^{2}}}\approx b(r){\displaystyle \int }_{-\tfrac{{\rm{\Delta }}}{2}}^{\displaystyle \frac{{\rm{\Delta }}}{2}}d\mu \\ \quad \left(1+\displaystyle \frac{x\mu }{{\sigma }^{2}}+\displaystyle \frac{1}{2}\displaystyle \frac{{x}^{2}{\mu }^{2}}{{\sigma }^{4}}\right)\displaystyle \frac{{e}^{-\tfrac{1}{2}\tfrac{{\mu }^{2}}{{\sigma }^{2}}}}{{\rm{\Delta }}}\\ \quad \propto b(r)\left(1+\displaystyle \frac{1}{12}\displaystyle \frac{{x}^{2}{{\rm{\Delta }}}^{2}}{{\sigma }^{4}}\right)\approx {e}^{-\tfrac{1}{2}\left(\tfrac{{x}^{2}}{{\sigma }_{x,\mathrm{eff}}}+\displaystyle \frac{{y}^{2}}{{\sigma }^{2}}\right)}\end{array}\end{eqnarray} \tag{ A5 }$$

with ${\sigma }_{x,\mathrm{eff}}^{2}={\sigma }^{2}+\tfrac{1}{6}{{\rm{\Delta }}}^{2}$. Since only the x direction was smeared, the beam is now slightly elliptical. For the purposes of S/N, what matters is not the shape of the beam, but its area, which has gone from 2π σ² to 2π σ σ_x,eff. We can use this to define an effective overall beam standard deviation:

$$\begin{eqnarray}&&\begin{array}{l}{\sigma }_{\mathrm{eff}}^{2}=\sigma {\sigma }_{x,\mathrm{eff}}={\sigma }^{2}\sqrt{1+\displaystyle \frac{1}{6}\displaystyle \frac{{{\rm{\Delta }}}^{2}}{{\sigma }^{2}}}\approx {\sigma }^{2}+\displaystyle \frac{1}{12}{{\rm{\Delta }}}^{2}\\ \quad \iff \,{\mathrm{FWHM}}_{\mathrm{eff}}^{2}={\mathrm{FWHM}}^{2}+\displaystyle \frac{2\mathrm{log}(2)}{3}{{\rm{\Delta }}}^{2}\end{array}\end{eqnarray} \tag{ A6 }$$

which happens to be the same result as what we got for circular smearing.

For each r we measure the mean μ_spat(r) and standard deviation σ_spat(r) of the S/N map as a function of position. ⁵⁷ Using these, we define the spatially normalized detection statistic z₁ as

$$\begin{eqnarray}&&{z}_{1}(r)=[{\rm{S}}/{\rm{N}}(r)-{\mu }_{\mathrm{spat}}(r)]/{\sigma }_{\mathrm{spat}}(r).\end{eqnarray} \tag{ B7 }$$

This normalization process is shown in Figure 11.

Zoom In Zoom Out Reset image size

We then build the empirical survival function N(z₁ > x) for all peaks in the tile with z₁ > 1 and map this to a corresponding Gaussian quantile

$$\begin{eqnarray}&&{z}^{* }=-\sqrt{2}\ {\mathrm{erf}}^{-1}(2N/n-1).\end{eqnarray} \tag{ B8 }$$

The value of n controls how far into the tail of a Gaussian survival function we map our empirical survival function. Its exact value is not important as long as $n\gt \max (N)$. It should be kept constant for all values of r to ensure that equally rare values of z₁ map to the same z^* for all distances. We chose $n={A}_{\mathrm{tile}}/{A}_{\mathrm{peak}}\approx {(12^\circ /2^{\prime} )}^{2}\,=\,1.3\cdot {10}^{6}$, with A_tile = (12°)² being the area of the tile, and ${A}_{\mathrm{peak}}={(2^{\prime} )}^{2}$ being the approximate feature size in the z₁ map. In principle Equation (B8) could be used to directly normalize z₁, but in practice there are too few samples in the tail. However, z^* turns out to be very well approximated as a linear function of z₁. ⁵⁸ We use this to define the fully normalized detection statistic

$$\begin{eqnarray}&&z(r)=[{z}_{1}(r)-{\mu }_{\mathrm{dist}}(r)]/{\sigma }_{\mathrm{dist}}(r)\end{eqnarray} \tag{ B9 }$$

where μ_dist and σ_dist are the offset and slope of the function z^*(z₁). This process is illustrated in Figure 12. Inserting the expression for z₁ into Equation (B9), we get the full normalization

$$\begin{eqnarray}&&z(r)=({\rm{S}}/{\rm{N}}(r)-{\mu }_{z}(r))/{\sigma }_{z}(r)=\xi ({\rm{S}}/{\rm{N}})\end{eqnarray} \tag{ B10 }$$

where we have defined μ_z = μ_spat + σ_spat μ_dist and σ_z =σ_spat σ_dist; and we have implicitly defined the function ξ.

Zoom In Zoom Out Reset image size