Trans-Neptunian Objects Found in the First Four Years of the Dark Energy Survey

The DES observational strategy is fully described in Diehl et al. (2016, 2018) but we summarize the relevant details here. The "wide" portion of the survey aims to completely tile the 5000 deg² footprint with 10 × 90 s exposures in each of the g, r, i, and z bands, and 6 × 45 + 2 × 90 s in the Y band, such that the completed wide survey will comprise ≈80,000 DECam exposures. DES observations are taken in seasons beginning in August and ending in February each year. The wide survey exposures are interleaved with a "deep" survey which images 10 DECam pointings (≈30 deg²) at ≈ weekly intervals in the griz bands, primarily for detection and measurement of high-redshift Type Ia supernovae (Bernstein et al. 2012). The Y4 TNO search reported herein was conducted only on the wide-survey images, but we have also searched much of the deep data and reported TNO detections to the Minor Planet Center (MPC) (Gerdes et al. 2016, 2017; Becker et al. 2018; Khain et al. 2018, 2020; Lin et al. 2019).

Observations are scheduled with the obstac algorithm (Neilsen & Annis 2013), which works to optimize the quality and homogeneity of the final wide-survey products in the face of variable clouds, seeing, and moonlight, while balancing the needs of the deep survey for temporally regular sampling. Importantly, part of this optimization is to avoid imaging the same part of the wide survey more than once in any filter on any single night (the goal is to spread the weather fluctuations across the footprint to homogenize the final survey products). Occasionally obstac will elect to take successive exposures at the same pointing in two different filters, typically g, r, and/or i in dark time and zY in bright conditions. These pairs, with only two minute intervals between them, are not useful for detecting TNO motion. Repeat imaging on any time interval longer than two minutes and shorter than one day is very rare. When we calculate the robustness of our TNO detections against false-positive linkages, we will always consider multiple detections on the same night to be fully correlated (not independent) since most sources of false positives (asteroids, variable stars, artifacts from bright sources) will tend to repeat in successive exposures. Each DES exposure is processed and evaluated the next day. Exposures failing to meet minimum standards for seeing and signal-to-noise ratio (S/N) are discarded, and their pointings placed back on the queue. We do not search rejected images for TNOs.

Nominally the survey goals were to cover half of the footprint with four "tilings" per filter in the first year, cover the other half with four tilings in the second year, and then two full tilings per year for three succeeding years. Variability in weather and availability of fields meant that this plan is only a rough approximation to reality. In particular, very poor weather in year 3 put the survey behind schedule, and an additional half-season's worth of observations were allocated for a sixth year to complete the survey. In the Y4 data under analysis here, a typical point in the footprint has been within the DECam field of view for around seven exposures per band. Within the field of view, 15%–20% of the sky is lost to gaps between the DECam CCDs, regions of defective CCDs, and area lost in the glare of bright sources. Thus a typical TNO will be cleanly imaged ≈6× per filter in the Y4 data.

One other observational detail of note is that DECam has a two-blade shutter which takes ≈1 s to sweep across the focal plane (Flaugher et al. 2015, Section 6.1). Successive exposures sweep the blades in opposite directions. We assign all TNO detections to a time corresponding to the midpoint of the shutter-open interval in the center of the focal plane. The true mean time of the open-shutter interval for any particular TNO's exposure may differ by up to ≈0.5 s in an unknown direction. Since the TNO apparent motion in this time interval is <1 mas, the shutter uncertainty contributes negligibly to errors in orbit determination.

For the analysis presented here, we make use of two distinct types of object catalogs produced during standard DES processing (Morganson et al. 2018): the single-epoch (SE) produced for each individual exposure, and the multi-epoch (a.k.a. "coadd") catalogs produced from an image averaging all exposures.

The SE catalogs used for this search come from the "Final Cut" reduction of the internal "Y4A1" data release. This processing includes the detrending and calibration of the raw CCD data. Artifacts such as saturated pixels close to bright stars, cosmic rays, and streaks are masked, sky background template images are subtracted from the science image, and the point-spread function (PSF) is modeled via the bright stars. Sources are detected in the images using SExtractor (Bertin & Arnouts 1996), with its default 3 × 3 pixel filter first being applied to the image. The detection criteria are that a source needs 6 pixels (${\mathtt{DETECT}}\_{\mathtt{MINAREA}}=6$) above the detection threshold ${\mathtt{DETECT}}\_{\mathtt{THRESH}}={\mathtt{ANALYSIS}}\_{\mathtt{THRESH}}=1.5$. A variety of position, flux, and shape measurements are made on each object and recorded in the output SE catalog.

The coadd catalogs in use here come from an internal data release labeled "Y3A2," i.e., they combine images only from the first three DES seasons. Sources are detected on a sum of all the r, i, and z images using SExtractor settings similar to the SE values above. The coadd catalogs and images used in this work are the same as those in the public DR1 release³¹ (Abbott et al. 2018).

The DES catalogs are exquisitely well calibrated to a common photometric system across the focal plane (Bernstein et al. 2017a) and across all the exposures of the survey (Burke et al. 2017). Comparison of DES magnitudes to Gaia DR2 magnitudes show uniformity across the footprint to ≈6 mmag rms (Abbott et al. 2018, Section 4.2). Trailing of TNO images in the 90 s DES wide exposures has negligible impact on photometry: at the maximual rate of apparent motion for objects ≥30 au of 5'' hr^–1, the object would move 012. The second moment of the trail is then <1% of the second moment of the PSF, and the point-source calibrations should be accurate to <0.01 mag even for this maximal apparent motion.

Every DES exposure has an astrometric map from pixel coordinates to the Gaia DR1 (Lindegren et al. 2016; Brown et al. 2016) celestial reference frame. Here again we benefit from extensive DES calibration efforts. Bernstein et al. (2017b) describes the DECam astrometric model, demonstrating that all distortions induced by the telescope, instrument, and detectors are known to 3–6 mas rms. We apply this model to all exposures in the Y4 DES TNO search. The astrometric uncertainties for bright TNO detections are dominated by stochastic displacements caused by atmospheric turbulence, at 10–15 mas rms on typical exposures. Below we describe the estimation of this turbulence uncertainty for each exposure. Fainter detections' astrometric uncertainties are dominated by shot noise in the centroiding of the image.

The astrometric solution includes terms for differential chromatic refraction (DCR) in the atmosphere and lateral color in the optics, which are calibrated in terms of the g − i color of the source. Similarly the photometric solution includes color terms. The maximal amplitude of the DCR (for airmass X < 2) and lateral color are ≈80 and 40 mas per mag of g − i color, respectively, for g-band observations, and five or more times smaller in other bands (Bernstein et al. 2017b). The g − i colors of all stars are measured directly by the DES, which fixes the reference frame, but the TNO apparent position will depend on its (unknown) color. Our TNO search is executed using positions that assume a nominal color, g − i = 0.61 (a typical stellar color). Only after a TNO is linked can we estimate its color. The final positions, magnitudes, and orbital parameters reported herein are calculated after refinement using the proper chromatic corrections.

The procedure for estimating the atmospheric turbulence contribution to astrometric errors is as follows.

• i.  
  For each high-S/N star i in the survey footprint, we calculate a mean position by averaging the positions predicted by the astrometric model in all DES exposures, as well as any position available in Gaia DR1, to produce a "truth" value ${\bar{{\boldsymbol{x}}}}_{i}$.
• ii.  
  Restricting ourselves to the positions ${{\boldsymbol{x}}}_{i}$ measured on an individual DES exposure, we find the measurement error ${\rm{\Delta }}{{\boldsymbol{x}}}_{i}\equiv {{\boldsymbol{x}}}_{i}-{\bar{{\boldsymbol{x}}}}_{i}=({\rm{\Delta }}{x}_{i},{\rm{\Delta }}{y}_{i})$ from the mean position. The displacement is measured in a local gnomonic projection, with x pointing to local equatorial east and y north.
• iii.  
  A cubic polynomial function of field coordinates is fit to the ${\rm{\Delta }}{{\boldsymbol{x}}}_{i}$ and adopted as the large-scale distortions from atmospheric and turbulent refraction for this exposure. This fit is subtracted from the ${\rm{\Delta }}{{\boldsymbol{x}}}_{i}.$
• iv.  
  We calculate the two-point correlation functions of astrometric error, averaging over all pairs of stars in the exposure versus their separation r:

  $$\begin{eqnarray}&&{\xi }_{x}(r)=\langle {\rm{\Delta }}{x}_{i}{\rm{\Delta }}{x}_{j}{\rangle }_{| {{\boldsymbol{x}}}_{i}-{{\boldsymbol{x}}}_{j}| =r},\end{eqnarray} \tag{ 1a }$$

  $$\begin{eqnarray}&&{\xi }_{y}(r)=\langle {\rm{\Delta }}{y}_{i}{\rm{\Delta }}{y}_{j}{\rangle }_{| {{\boldsymbol{x}}}_{i}-{{\boldsymbol{x}}}_{j}| =r},\end{eqnarray} \tag{ 1b }$$

  $$\begin{eqnarray}&&{\xi }_{\times }(r)=\langle {\rm{\Delta }}{x}_{i}{\rm{\Delta }}{y}_{j}{\rangle }_{| {{\boldsymbol{x}}}_{i}-{{\boldsymbol{x}}}_{j}| =r}.\end{eqnarray} \tag{ 1c }$$

  An example of the behavior of ξ is shown in Figure 2.
• v.  
  We assign a 2D Gaussian positional uncertainty to each detection in the exposure, with a covariance matrix given by

  $$\begin{eqnarray}&&{{\boldsymbol{\Sigma }}}_{\mathrm{atm}}=\left(\begin{array}{l}\langle {\xi }_{x}\rangle \langle {\xi }_{\times }\rangle \\ \langle {\xi }_{\times }\rangle \langle {\xi }_{y}\rangle \end{array}\right)\end{eqnarray} \tag{ 2 }$$

  where the ξ values are averaged over the separation range 24'' < r < 40'' where they typically plateau. Note that this use of the correlation function at r > 0 removes any contribution due to shot noise of individual stellar measurements.

Zoom In Zoom Out Reset image size

Figure 3 shows the histogram of atmospheric turbulence strength in the exposures. The typical atmospheric turbulence is seen to be ≈10 mas. The astrometric errors are typically substantially anisotropic on most individual exposures, with ${{\rm{\Sigma }}}_{{xx}}\ne {{\rm{\Sigma }}}_{{yy}}$ and ${{\rm{\Sigma }}}_{{xy}}\ne 0,$ because astrometric errors are substantially different parallel versus perpendicular to the prevailing wind direction.

Zoom In Zoom Out Reset image size

The final astrometric uncertainty assigned to each transient is the quadrature sum of the atmospheric turbulence ellipse in Equation (2) with the circularly symmetric shot noise error in the centroid reported by SExtractor.

The computational burden of linking TNOs from the transient catalog scales roughly as n³, where n is the density of transient detections per square degree per exposure. It is therefore of great importance to minimize n while retaining all true TNO detections in the transient catalog. Because asteroids are indistinguishable from TNOs at the single-exposure level, the transient catalog must include all detected asteroids, which will greatly outnumber the TNOs.³² Asteroids thus form an irreducible source of false-positive TNO transient detections, and our goal, therefore, is to reduce any other spurious transient detections to a level well below the asteroid density.

We wish to identify all sources in the thousands of DES images that are present in a given sky location for only a single night. We begin by matching the SE in all bands and coadd catalogs (in which each detection has information from all bands) by sky coordinates. All detections in a small region are projected into a tangent plane and then grouped with a kD-tree (Maneewongvatana & Mount 1999) friends-of-friends (FoF) algorithm that links detections within 05 of each other.

For each group of SE detections, we calculate the following quantities:

• 1.  
  whether or not it is matched to a coadd detection, ${\mathtt{COADD}}=1,0;$
• 2.  
  how far apart in time are the first and the last SE detections of the match group, ${\rm{\Delta }}t\equiv {t}_{\mathrm{last}}-{t}_{\mathrm{first}};$
• 3.  
  how much fainter or brighter the SE detection is compared to the coadd detection in the same band (when present), ${\rm{\Delta }}m\equiv {m}_{\mathrm{SE}}-{m}_{\mathrm{coadd}}$.

If the SE detection is a detection of a solar system object, we do not expect any flux to be present at this location on the other exposures contributing to the coadd taken more than a few hours away from the SE detection. If there are K exposures in the coadd, then the averaging process will reduce the apparent flux of the source in the coadd by a factor 1/K, leaving the coadd source with a magnitude fainter than the SE value by $2.5{\mathrm{log}}_{10}K\geqslant 0.75$ mag for K ≥ 2. The coadd source should thus either be absent or significantly fainter than the SE source. We therefore implement the following cut to retain potential solar system measurements while eliminating many variable fixed (stellar) sources:

• 1.  
  ${\mathtt{COADD}}=1$, Δt < 2 days, Δm ≤ −0.4; OR
• 2.  
  ${\mathtt{COADD}}=0$, Δt < 2 days.

Note that Δt < 2 days is by itself an insufficient condition, since it would include many faint sources that are pushed above detection thresholds on only 1–2 exposures because of intrinsic variability or noise fluctuations. The veto by the more sensitive coadd images solves this issue, with the Δm ≤ −0.4 threshold estimated empirically to include nearly all truly moving sources (TNOs), as illustrated in Figure 4.

Zoom In Zoom Out Reset image size

Unlike Goldstein et al. (2015) or Lin et al. (2018), we do not make cuts on object sizes or other morphological features of detections, as these measures are too noisy to be useful discriminants of non-stellar artifacts from TNOs for sources at the threshold of detectability. We are willing to accept higher false-positive rates in order to keep the lowest flux threshold for true TNOs.

The first panel of Figure 5 maps the sky locations of transient detections, revealing many coherent structures that cannot be produced by TNOs, asteroids, or other true celestial moving objects. These structures correspond to image artifacts that were not successfully masked by the SE pipeline, such as meteor/satellite/airplane streaks, bad pixel columns, and artifacts from detections around bright sources. There are also some exposures with incorrect astrometric solutions, creating large numbers of spurious sources which do not match sky coordinates of other exposures' detections. We find empirically that the following steps applied to each exposure serve to reduce the number of spurious transient detections to level comparable to or below the irreducible asteroid density.

• i.  
  The detections are clustered with FoF linking length of 60''.
• ii.  
  Groups with more than 20 detections are removed from the catalog, as these correspond to entire CCDs or even exposures being identified as transient sources.
• iii.  
  Groups that have between 10 and 20 detections could come from either accidental matching of true transients, or from coherent structures like unmasked streaks or arcs around bright stars.
• iv.  
  To remove the latter, we test for a tendency of the points to lie along a line or curve: we evaluate the "bend angle," that is, the angle between the lines connecting each detection on this group to its two nearest neighbors. These angles are clipped to the range 0°–90°. If the average angle of all detections in the cluster is less than 15°, these detections are also removed from the catalog;
• v.  
  We also remove transients within 30 pixels of any edge of a CCD, since source images cut off at CCD edges have mis-measured centroids that cause spurious mis-matches with their corresponding coadd sources.

Figure 6 illustrates the three cases described here. The second panel of Figure 5 shows the identified structures in a region of the survey, and the last panel shows the output, and final, transient catalog. Some spurious structures remain, but since the number of these is now well below the number of asteroids, there is little to be gained by further cuts on the transient catalog.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

There are ≈7 × 10⁹ SE detections in the 60,000 Y4 exposures. Of these, 2 × 10⁷ are designated as transients (see Table 1).

Figure 1 presents a map of the density of bright (r < 23) transients for each DES Y4 r band exposure. At these magnitudes, most exposures are virtually complete (that is, the magnitudes of 50% completeness of each exposure is on average larger than this limit; see Section 2.6 for a discussion on the completeness estimates), and the drop in density from exposures with $\sim 200\,\mathrm{transients}\,{\deg }^{-2}$ close to the ecliptic plane to $\lt 100\,\mathrm{transients}\,{\deg }^{-2}$ far from it suggests that asteroids are comparable in the catalog to other astrophysical transients and artifacts.

This density is consistent with the sky density of 210 asteroids deg⁻² brighter than m_R = 23 reported by SKADS (Gladman et al. 2009) for observations within a few degrees of opposition. DES may encounter fewer objects given that observations are not typically at opposition. We note that this density is consistent with the latitudinal asteroid density, computed from transforming the inclination distribution for all asteroids in the MPC³³ to an ecliptic latitudinal distribution following Brown (2001), plus a background level of 45 transients deg⁻² independent of ecliptic latitude, as presented in the lower left panel of Figure 7.

Zoom In Zoom Out Reset image size

We do not investigate the source of this background. They are not supernovae: from simulations of DES supernovae transient detections (Kessler et al. 2019), using supernovae rates from Jones et al. (2018), we expect only order unity supernovae (Ia and core collapse) per exposure to have r < 23, and yet they appear in only one night's exposure.

The upper right-hand panel of Figure 7 shows the average transient density in all exposures for all transients and, while the lower density far from the ecliptic is still present in the griz bands, the density plateaus at ∼200 per square degree per exposure (r band), defining the minimum false-positive rate for our search. The plateau at $\sim 150\,\mathrm{transients}\,{\deg }^{-2}$ on the faint end corresponds to astrophysical transients (for example, stellar outbursts and supernovae that passed the Δt cut, as well as high-inclination asteroids), faint sources detected only once that passed the Δm cut, as well as image artifacts (e.g., cosmic rays, unmasked streaks). Image inspection on a randomly selected subset of false positives (identified in false TNO linkages; see Section 4) indicates that most of the background comes from unmasked image artifacts. While the single-night transient catalog is not pure, the level of false positives is low enough that the search is feasible.

To characterize the detection efficiency for point sources in a given exposure, we search the coadd catalogs for unresolved sources overlapping the exposure's footprint. Here, unresolved sources are defined as those with $| {\mathtt{SPREAD}}\_{\mathtt{MODEL}}\_{\mathtt{I}}| \lt 0.003$ (see Section 8.1 of Drlica-Wagner et al. 2018). We then note each coadd source's magnitude in the band of the exposure, and record whether the source was detected in the SE processing of that exposure. This yields a list ${m}_{\det ,i}$ of the "true" (coadd) magnitudes of SE-detected point sources and another list ${m}_{\mathrm{non},j}$ of non-detected point sources in the exposure. If we posit a probability of detection p(m) for this exposure, then the total probability of the observed outcome is

$$\begin{eqnarray}&&{p}_{\mathrm{tot}}=\prod _{i}p\left({m}_{\det ,i}\right)\prod _{j}\left[1-p\left({m}_{\mathrm{non},j}\right)\right].\end{eqnarray} \tag{ 3 }$$

We posit that the completeness function for the exposure takes the form of a logit function

$$\begin{eqnarray}&&p(m)=\displaystyle \frac{c}{1+\exp (k(m-{m}_{50}))},\end{eqnarray} \tag{ 4 }$$

and we adjust the parameters m₅₀, c, and k to maximize p_tot. These then define the completeness function we adopt for this exposure. Figure 8 shows the distribution of m₅₀ values for the exposures included in the Y4 wide-survey TNO search, as well as an example of the fit in one exposure.

Zoom In Zoom Out Reset image size

For short arcs of distant objects, we can neglect the gravitational perturbation, and the effect of $\dot{\gamma }$ on observed position is highly suppressed. Rewriting Equations 5(a) and (b) with $g,\dot{\gamma }=0$,

$$\begin{eqnarray}&&\alpha +\dot{\alpha }(t-{t}_{0})=(1-\gamma {z}_{{\rm{E}}}){\theta }_{x}+\gamma {x}_{{\rm{E}}},\end{eqnarray} \tag{ 9a }$$

$$\begin{eqnarray}&&\beta +\dot{\beta }(t-{t}_{0})=(1-\gamma {z}_{{\rm{E}}}){\theta }_{y}+\gamma {y}_{{\rm{E}}}.\end{eqnarray} \tag{ 9b }$$

If we assume a distance d = 1/γ, this allows transformation of the $({\theta }_{x},{\theta }_{y})$ coordinates into an (α, β) system where the motion is linear. Holman et al. (2018) exploit this result for linking of tracklets, but we must link individual detections. From Equation (8), it becomes clear that after a time Δt an object in a bound orbit will be in a circle in the (α, β) plane of radius

$$\begin{eqnarray}&&{r}_{2}(\gamma ,{\rm{\Delta }}t)=\sqrt{2{{GM}}_{\odot }{\gamma }^{3}}{\rm{\Delta }}t\end{eqnarray} \tag{ 10 }$$

centered around the first detection.

We start by selecting all exposures from season k which could potentially contain a TNO with γ within the selected bin and (α, β) within the 35 patch radius. We apply the transformations from Equations 9(a) and (b) for γ_j to all transients detected in each exposure, and each set of resulting positions is used to constructed a kD tree. Trees for exposure pairs that are up to 90 days apart from each other are then searched for pairs of detections with a search radius as defined in Equation (10). This process is repeated at the edges of the distance bin (${\gamma }_{j}\pm \delta {\gamma }_{j}$), resulting in the final list of pairs for the bin j.

The expected number of pairs of unassociated transients grows with the area of the search circles and the time interval between pairs. The mean number of such false-positive pairs between an exposure μ and all later exposures ν is

$$\begin{eqnarray}&&{N}_{2,\mu }(\gamma )=\sum _{\nu \gt \mu }2\pi {{GM}}_{\odot }{\gamma }^{3}{\rm{\Delta }}{t}_{\mu \nu }^{2}{n}_{\mu }{n}_{\nu }{A}_{\mathrm{DECam}}\end{eqnarray} \tag{ 11 }$$

where n_μ is the sky density of transients in exposure μ, Δt_μν is the time interval between these two, and A_DECam is the imaging area of DECam. (This equation ignores the case where the search circle only partially overlaps exposure ν.) The total number of pairs expected in the search is ${N}_{2}={\sum }_{\mu }{N}_{2,\mu }.$ The Y4 search yields ≈10¹² pairs (Table 1).

With a pair of detections we can compute $\dot{\alpha }$ and $\dot{\beta }$ as a function of γ and $\dot{\gamma }$. A pair's nominal expected position at future exposure at time t is determined by Equations 5(a) and 5(b) with α, β, $\dot{\alpha }$, and $\dot{\beta }$ computed for the bin center ${\gamma }_{j}$ and $\dot{\gamma }=0$. We further assume that the future position is approximately linear in $(\gamma -{\gamma }_{j})$ ("parallax" axis) and $\dot{\gamma }$ ("binding" axis) as long as we remain within the distance bin range ${\gamma }_{j}\pm \delta {\gamma }_{j}$ and at $\dot{\gamma }$ small enough to maintain a bound orbit.

To compute the line of variation in (θ_x, θ_y) due to γ deviations (parallax), we calculate the positions for orbits corresponding to α, β, $\dot{\alpha }$ and $\dot{\beta }$ computed for $\gamma ={\gamma }_{j}\pm \delta {\gamma }_{j}$ and $\dot{\gamma }=0$.

To compute the line of variation in (θ_x, θ_y) due to non-zero $\dot{\gamma }$ (binding), we first make the assumption that the target orbit is bound, in which case we can write the system's energy more carefully, and derive (as per BK)

$$\begin{eqnarray}\begin{array}{rcl}{\dot{\gamma }}^{2}\leqslant {\dot{\gamma }}_{\mathrm{bind}}^{2} & \equiv & 2{{GM}}_{\odot }{\gamma }^{3}\\ & & \times {\left(1+{\gamma }^{2}-2\gamma \cos {\beta }_{0}\right)}^{-1/2}-{\dot{\alpha }}^{2}-{\dot{\beta }}^{2},\end{array}\end{eqnarray} \tag{ 12 }$$

where β₀ is the target's solar elongation. The line of variations is then derived from positions for orbits corresponding to α, β, $\dot{\alpha }$ and $\dot{\beta }$ computed for γ = γ_j and $\dot{\gamma }=\pm {\dot{\gamma }}_{\mathrm{bind}}$.

Figure 9 illustrates how the search region for triplet candidates at some time t is the parallelogram constructed from the parallax axis of γ line of variations and the binding axis of $\dot{\gamma }$ line of variation.

Zoom In Zoom Out Reset image size

For each pair, we search for triplets in exposures up to 90 days after the second detection. For each exposure that crosses the region of potential third transients, we use its kD tree to rapidly locate for transients lying in a set of circles that cover the parallelogram of potential orbit extrapolations (plus a small contribution for position measurement errors).

One can show that the leading order behavior in γ of the total number of spurious triplets (from randomly positioned transients) is

$$\begin{eqnarray}&&{N}_{3}(\gamma )\simeq {{nN}}_{2}(\gamma ){A}_{3}(\gamma )\propto {n}^{3}{\gamma }^{7/2},\end{eqnarray} \tag{ 13 }$$

where A₃ is the area of the search parallelogram.

The Y4 search generates 6 × 10¹⁰ triplets (Table 1).

Once a triplet of detections is identified, we can in principle fit an orbit with all six orbital parameters left free. We follow routines similar to those in BK to fit the orbits and determine expected positions and their linearized uncertainties for the circumstances of any DES exposure. We do not make use of the prior in the BK Equations (20) and (21) that favors nearly circular orbits, because it can cause convergence problems for some of the very distant and nearly parabolic TNO orbits that we wish to search. Instead we institute a Gaussian prior on (inverse) distance to force the orbit into our distance bin. The prior contributes to the fit χ² as

$$\begin{eqnarray}&&{\chi }_{\mathrm{prior}}^{2}={s}_{\gamma }\displaystyle \frac{{(\gamma -{\gamma }_{j})}^{2}}{\delta {\gamma }_{j}^{2}},\end{eqnarray} \tag{ 14 }$$

where the strength of this prior is adjusted by setting s_γ, and we choose s_γ ≈ 1. We also include a prior for bound orbits, defined as

$$\begin{eqnarray}&&{\chi }_{\mathrm{bind}}^{2}=b\displaystyle \frac{{\dot{\alpha }}^{2}+{\dot{\beta }}^{2}+{\dot{\gamma }}^{2}}{2{{GM}}_{\odot }{\gamma }^{3}}=\displaystyle \frac{| \mathrm{KE}| }{| \mathrm{PE}| }.\end{eqnarray} \tag{ 15 }$$

Here b, the binding factor, defines the strength of this prior. Note that this is equal to b for a parabolic orbit. We set b = 4 for our fits. The ${\chi }_{\mathrm{bind}}^{2}$ is added to the quantity being minimized in the BK code.

The growing process starts with the orbit fitted to a triplet. The position and error ellipse are calculated at all other exposures that the TNO might cross. The kD tree for each potential exposure is searched for transients that lie within the 4σ error ellipse defined by convolving the error ellipse of the orbit prediction at the time of the search exposure with the measurement uncertainty of each transient position on the exposure. For each such transient found, a new n-let is defined, and is queued for its own orbit fit. The addition of transients is iterated for each n-let until no new transients are found to be consistent with the orbit. We also discard any n-let whose orbit fit is significantly unbound, whose best-fit orbit has unacceptably high χ² value, or which duplicates a set of transients that have already been examined.

Each time we are attempting to grow an n-let, we calculate the FPR for the addition of the (n + 1)th transient. If the 4σ error ellipse on exposure j has area ${A}_{j,\mathrm{search}},$ and the density of transients on this exposure is n_j, then the total probability of a spurious linkage is

$$\begin{eqnarray}&&\mathrm{FPR}=\sum _{j}{A}_{j,\mathrm{search}}{n}_{j},\end{eqnarray} \tag{ 16 }$$

where the sum is over all exposures being considered for the (n + 1)th transient. Our calculation of FPR currently accounts for CCD gaps and other details of geometry only approximately, but this will be sufficient. Note that as the orbit becomes better defined from higher n and longer arc, the FPR for further additions will shrink. If an n-let has a short arc and we are searching for transient n + 1 in a different season's exposure, then the error ellipse may be large, and we will be flooded with false linkages. We therefore do not search an exposure if the contribution of this exposure to the FPR sum would be ${A}_{j,\mathrm{search}}{n}_{j}\gt 10.$ Since the transient catalog density is $\lesssim 200\,{\deg }^{-2},$ this is roughly a requirement that the orbit be localized to an area of $\lesssim 0.05\,{\deg }^{2}$ if we are to proceed. We do not implement this FPR cutoff until the n-let has linked detections on ≥4 distinct nights, since the ${A}_{j,\mathrm{search}}$ is unavoidably large with ≤3 points on the arc.

The growing process terminates when no additional detections are found to match the orbit fit to an n-let. We associate with this terminal n-let the FPR that was calculated for the final successful linkage step. Thus the FPR recorded with the n-let estimates the probability that the last detection linked to the orbit is spurious, i.e., a transient that randomly fell within the orbit's error ellipse. In the future this information will allow us to estimate the number of spurious linkages contaminating our TNO catalog. At present, we only make a very mild cut of discarding individual n-lets with FPR > 1. If a terminal n-let has detections in seven or more unique nights and has FPR < 10⁻³, we call it a secure orbit, and we remove all of its detections from consideration for linkage to further orbits.

At the end of the linking process, we have to decide which sets of linked detections reliably correspond to multiple detections of the same real solar system object, versus spurious linkages of mixtures of detections of distinct sources or artifacts. In this Y4 search, we use the following criteria to cull the set of unique terminal n-lets to high-reliability candidates.

• 1.  
  If the number of unique nights on which these detections were seen is ${\mathtt{NUNIQUE}}\lt 6,$ the n-let is rejected. Of the 132 linkages with ${\mathtt{NUNIQUE}}=6$ meeting the other criteria, only 52 were confirmed as real by the method of Section 4. The number of n-lets with ${\mathtt{NUNIQUE}}=5$ is large and certainly dominated by spurious linkages. We defer until the final DES search an effort to extract a reliable candidate set from these. A total of 1684 distinct linkages satisfy this criterion.
• 2.  
  The time between first and last detections (${\mathtt{ARC}}$) must satisfy ${\mathtt{ARC}}\gt 6$ months, i.e., we demand detections in distinct DES observing seasons. More stringently, we define ${\mathtt{ARCCUT}}$ to be the shortest arc that remains after eliminating any single night of detections, and demand also that ${\mathtt{ARCCUT}}\gt 6$ months. In effect this means that at least two detections must occur outside the season of the triplet that gestates the n-let. Six linkages have ${\mathtt{ARC}}\lt 6$ months, and another 1235 also have ${\mathtt{ARCCUT}}\lt 6$ months.
• 3.  
  The ${\chi }^{2}$ per degree of freedom from the (prior-less) orbit fit must satisfy ${\chi }^{2}/\nu \lt 4$, which rejects another 19 linkages.
• 4.  
  The FPR is less than 1. No surviving linkages fail this criterion.

A total of 424 distinct linkages pass these criteria in the Y4 search (see Table 1). The "sub-threshold confirmation" technique described next will serve to remove spurious linkages from the candidate list defined by the above criteria.

The process for simulating the transient-catalog entries for a TNO of given orbital parameters and magnitude m is as follows.

• i.  
  We find all exposures for which the DECam FOV contains the TNO position given its orbit.
• ii.  
  An observed position for the TNO is derived by adding observational error to the position predicted by the orbit. Currently the observational error is drawn from the distribution of errors for all point sources in the image (including both shot noise and the atmospheric turbulence contribution per Equation (2)); in the future we will properly track the errors versus source magnitude.
• iii.  
  We check whether the "observed" position lies within a functional DECam CCD, i.e., we would have collected an image of this TNO. If this is true, this is considered an observation.
• iv.  
  We compute the value of p(m) from Equation (4) with the fitted parameters for this exposure and compare to a random value r between 0 and 1. If p(m) > r, this observation is considered a detection and is entered into the transient catalog.
• v.  
  For completeness, a random magnitude error is drawn from the magnitude error distribution determined for this exposure and added to the truth magnitude, for each synthetic transient. Our linking algorithm makes no use of the detected magnitude, so this is irrelevant.

This process does not simulate the loss of TNOs due to potential overlap with other sources or image artifacts, which as noted above creates a loss of a few percent of sources. Note also that no light-curve variation is placed on the simulated sources.

We generate a population of synthetic TNOs that is intended to sample the full phase space of possibly detectable orbits and magnitudes, with no intention of mimicking the true TNO population. We generate the fakes by sampling the barycentric phase space $\{{x}_{0},{y}_{0},{z}_{0},{\dot{x}}_{0},{\dot{y}}_{0},{\dot{z}}_{0}\}$ at a reference time t₀ near the survey midpoint. To generate the position vector, we sample the unit sphere by constructing uniformly distributed angles with a Fibonacci lattice (see, e.g., González 2010) in equatorial coordinates, and discard all points well outside the DES footprint. Each fake is assigned a random barycentric distance from a distribution placing half of the fakes uniformly in the range 30–60 au and half logarithmically distributed between 60 and 2500 au. Similarly, the velocities are sampled by placing angles on a spherical Fibonacci lattice, and assigning a velocity $v={f}^{1/3}{v}_{\mathrm{esc}}(d),$ where v_esc(d) is the escape velocity at the barycentric distance, and f is a uniform deviate between 0 and 1.

Each fake is also assigned an r-band magnitude, independent of its distance, sampled from a uniform distribution between m_bright = 20 and m_faint = 24.5, spanning the range in which we expect our Y4 discovery efficiency to go from near unity to zero. Each TNO's magnitude is assumed to be constant (i.e., there is no light curve for the object). The colors are fixed and chosen to be similar (but not identical) to those of (136199) Eris (Brown et al. 2005) as observed by DES: g − r = 0.55, r − i = 0.07, r − z = −0.02, and r − Y = −0.04.

The transient catalog used in the search includes these fakes, and the process is blinded in that the linking algorithms make no distinction between real and fake transients.

Figure 13 shows (at left) the frequency of observations and detections for a population of ≈200,000 fakes. Here we can see that the typical TNO is observed on ≈20 distinct nights of the Y4 data, with a tail to low values for TNOs that move in or out of the survey footprint during the survey. This will occur near the edges and along the thin equatorial stripe of the footprint. Fainter objects are detected in fewer observations, once m_r > 21. In the range $23\lt {m}_{r}\lt 23.5,$ about half of fake TNOs are "recoverable" by our criterion that they have ${\mathtt{NUNIQUE}}\geqslant 6.$ The right-hand panel integrates over all orbital parameters to give the effective survey area as a function of m_r, indicating that the detection and linkage are highly complete for m_r < 23, and 50% complete at m_r ≈ 23.3 (assuming that all colors are Eris-like). Figures 14 shows the recovery efficiency as a function of inclination and barycentric distance, indicating that there is very little dependence on the orbital properties of an object at fixed apparent magnitude. The DES footprint is much broader than a typical TNO's orbital path, except for the narrow equatorial stripe, meaning that the detectability of a TNO is almost entirely a function of apparent magnitude once it has distance >30 au. A minor exception is for TNOs with inclination near 0° or 180°, for which a significant fraction of our coverage is in the narrow strip.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

A search of the MPC database for known TNOs that were within the FOV of at least six of the DES Y4 griz exposures shows that all such objects above our estimated r = 23.3 mag 50% completeness level are indeed among our 316 detections. The brightest of these known TNOs that we do not rediscover is 2013 RY₁₀₈, at r = 23.47, which was in fact discovered from the deep supernova-search images in DES. It is in the FOV for 13 wide-survey exposures on nine distinct nights, but was only detected in three of these exposures. This is fully consistent with our estimated completeness thresholds. The final DES TNO search will have a lower SE detection threshold, and more exposures, that should enable discovery of many similar TNOs.