The DECam Ecliptic Exploration Project (DEEP) VI: first multi-year observations of trans-Neptunian objects

We present the first set of trans-Neptunian objects (TNOs) observed on multiple nights in data taken from the DECam Ecliptic Exploration Project (DEEP). Of these 110 TNOs, 105 do not coincide with previously known TNOs and appear to be new discoveries. Each individual detection for our objects resulted from a digital tracking search at TNO rates of motion, using two to four hour exposure sets, and the detections were subsequently linked across multiple observing seasons. This procedure allows us to find objects with magnitudes $m_{VR} \approx 26$. The object discovery processing also included a comprehensive population of objects injected into the images, with a recovery and linking rate of at least $94\%$. The final orbits were obtained using a specialized orbit fitting procedure that accounts for the positional errors derived from the digital tracking procedure. Our results include robust orbits and magnitudes for classical TNOs with absolute magnitudes $H \sim 10$, as well as a dynamically detached object found at 76 au (semi-major axis $a\approx 77 \, \mathrm{au}$). We find a disagreement between our population of classical TNOs and the CFEPS-L7 three component model for the Kuiper belt.

object found at 76 au (semi-major axis a ≈ 77 au).We find a disagreement between our population of classical TNOs and the CFEPS-L7 three component model for the Kuiper belt.

INTRODUCTION
Beyond Neptune, there exist thousands of known bodies, called Trans-Neptunian Objects (TNOs).Due to current observational capabilities, the majority of the known TNOs have diameters of order 100 km, corresponding to the observational limits of the surveys of the last two decades (see Trujillo & The DEEP Collaboration 2023, for a review).TNOs can be separated into multiple dynamically-distinct subpopulations, including the classical Kuiper Belt, the scattering disk, resonant TNOs, and the detached population (Gladman et al. 2008).Explaining this structure has been an active subject of research, as these small bodies are important probes of the formation and evolution of the early Solar System (Nesvorný 2018;Gladman & Volk 2021), providing key observational constraints to models of an early giant planet instability in the Solar System (Tsiganis et al. 2005).The migration of the giant planets can lead to objects being captured in mean motion resonances (Malhotra 1995) or scattered into orbits with significant eccentricities and inclinations (Duncan & Levison 1997), and their present day orbital distribution are indicative of dynamical processes in the early Solar System (Morbidelli & Nesvorný 2020).
Since the detection of the first KBO in 1992 (Jewitt & Luu 1993), astronomical surveys have consistently added to the number of known TNOs (e.g.Petit et al. 2011;Sheppard et al. 2016;Bannister et al. 2018;Bernardinelli et al. 2022).As of June 2023, the JPL small body database 1 lists 4463 known TNOs.The majority of these objects have been discovered in surveys that detect objects in individual images, link those to detections in other images at different epochs, and use the collection of observations to determine or constrain the orbit of the discovered TNO.The depth of this type of survey is limited by the depth of the individual exposures: the motion of the objects restricts the duration of individual exposures, as once an object has moved more than the angular extent of the image quality or "seeing," trailing limits the resulting signal-to-noise ratio (S/N ).A number of deeper surveys, typically limited to narrow fields (e.g.Gladman et al. 1998;Allen et al. 2001;Bernstein et al. 2004;Fraser et al. 2008;Fuentes et al. 2009) or relatively short time baselines (e.g.Smotherman et al. 2021;Whidden et al. 2019), overcome this limitation using a technique known as "shift and stack" or "digital tracking".This technique works by combining the signal from a series of short, untrailed exposures along hypothesized sky plane trajectories.Because the orbits are not known a priori, a range of rates and directions of motion are tested.This allows the detection of sources significantly fainter than the single image detection limit.Here, we present the first multi-year (2019-2021) linkages of TNOs detected with digital tracking in the 30 deg 2 of the B1 fields of the DECam Ecliptic Exploration Project (DEEP, Trilling & The DEEP Collaboration 2023).These time spans lead to high-quality orbit determination that enable dynamical classifications (Gladman et al. 2008).The combination of the dynamical properties with the well-understood observational biases (Bernardinelli & The DEEP Collaboration 2023) make statistical studies of our objects possible.
In Section 2 we briefly discuss the data processing of these fields as well as the single-night observations.In Section 4, we discuss the methods used to link these single-night observations to multi-year arcs, then fit orbits to these arcs.In Section 5, we discuss the implication of our recovered TNOs and their orbital properties when compared to existing TNO distributions.In Section 6, we summarize our results and discuss future work.

DATA
As its acronym indicates, DEEP uses the Dark Energy Camera (DECam) on the 4m Blanco telescope at the Cerro Tololo Inter-American Observatory (CTIO) (Flaugher et al. 2015).The data are taken in the VR band, which has a similar band center to, for example, the Pan-STARRS r p1 band, but broader wavelength coverage.This makes it well suited to detecting faint sources, although it is less preferred for precision photometry.
The search area is separated into four quadrants, "A0", "A1", "B0" and "B1".Each quadrant is composed of a number of fields that spread out in a triangular pattern that follows the typical dispersion of on-sky motions of the trans-Neptunian over time (Trujillo & The DEEP Collaboration 2023).All of our fields are aligned with the invariable plane (the solid black line in Figure 1) of the Solar System, to maximize object discovery.Each of these quadrants was covered with a series of 3 square degrees pointings that were visited over multiple nights and years.A field visit in a night consists of a "long stare", a series of between 52 and 130 consecutive or nearly-consecutive 120 second exposures, Table 1.Fields, nights, RA, Dec, and number of observations (N obs ) for all of the nights in the DEEP B1 quadrant that are searched with KBMOD.RA and Dec correspond to the central RA and Dec of DECam for the first exposure in the long stare.Each exposure in each long stare is 120 seconds long.
Figure 1 shows the pointing of each field in the first night that field was observed in the 2021 observing campaign.Over various nights, the pointing center for any given gield can vary by up to 0.817 degrees, with a median value of 0.39 degrees, from other observations of that field.These intra-field offsets were designed in part to mitigate the effect of chip gaps as well as the gaps between fields.Within a long stare, the median standard deviation is σ RA = 0.629 arcsec and σ Dec = 0.070 arcsec, and the images are aligned to a common WCS, as required by the shift and stack procedure.

Data pre-processing
We use the Legacy Survey of Space and Time (LSST) Science Pipelines to process the DEEP B1 data (Jurić et al. 2017).For these data, we start with DECam raw (uncalibrated) image files, flat field images, and bias images.Each long stare has associated flats and biases from the same night, with the exception of the 2020-10-15 visit of field B1b and the 2020-10-19 visit of field B1d, where we use flats and biases from the visits in 2020-10-16 and 2020-10-18, respectively.We then use the LSST Science Pipelines to perform instrument signal removal (ISR), photometric calibration, astrometric calibration, and image differencing.Photometric calibration is performed using Pan-STARRs DR1 data (Flewelling et al. 2020) as a reference catalog, mapping the V R band to the reference r p1 band.However, we elected to be explicit in our reported magnitudes that we are using this wide-band filter, so in all references below we have the V R subscript when referring to magnitudes.Astrometric calibration is performed using Gaia DR2 data (Lindegren, L. et al. 2018) as a reference catalog.We will follow the convention in the LSST Science Pipelines and refer to these photometrically and astrometrically calibrated exposures as "calexps".At this point, cosmic rays are masked in the individual images, so that they do not introduce correlated noise in the shift and stack procedure.
After calibration, we inject a population of synthetic TNOs (or "fakes") into the calexps, as described in detail in Bernardinelli & The DEEP Collaboration (2023).The goal of this population is to test all physically plausible TNO orbits, rather than attempting to reproduce any real existing TNO population structure.To accomplish this, we create a joint distribution from a low-e, low-i Kuiper belt-like population that spans 30 < a < 80 au; a moderate-eccentricity population with 0 < e < 0.4 and inclinations uniformly spanning 0 deg to 90 deg; and an isotropic distribution with randomly-sampled Cartesian coordinates between 25 and 1000 au.We propagate these objects to the epoch of each image, adding objects that intersected with the plane of a CCD to the synthetic catalog.This creates a catalog of 5737 unique synthetic objects.These objects each overlap with a CCD image 1 to 9 times.We assign values for m V R from two uniform distributions, 20 < m V R < 24 and 24 < m V R < 28, with 20% and 80% of the total objects in the brighter and fainter distributions, respectively.As our goal is to measure our discovery efficiency as a function of magnitude, these uniform distributions ensure that we have enough sources at a given magnitude to properly reconstruct our discoverability.Note that the discovery efficiency is equivalent to the ratio of the number of detected objects to the number of injected objects, so it is independent of the shape of the underlying injected population -see Bernardinelli & The DEEP Collaboration (2023) for a thorough discussion.We propagate magnitudes with respect to the changing heliocentric and geocentric distances, such that the reference magnitude for the object coincides with the magnitude at the reference epoch for the orbital elements.Half of the synthetic population included simulated lightcurves with periods between 2 and 100 hours and amplitudes between 0 and 0.5 mag.
After injecting synthetic TNOs, we generate template images for the difference imaging procedure by coadding the long-stare images from each night.Notably for this dataset, we use the default LSST Science Pipelines coaddition statistic of MEAN, where the coadded image corresponds to the mean flux in each pixel that contributes to the stack.Smotherman et al. (2021) describes in more detail the the image differencing technique, as well as the astrometric performance of the LSST Science Pipelines applied to DECam.Because we generate coadded templates from the images of each long stare, we are self-subtracting some flux from the TNOs, both real and synthetic, in the individual images.This removal will depend on the speed of the TNO and the time baseline of the specific long stare.This has the effect of reducing our S/N (and thus reducing our effective depth) and complicating the magnitude measurement for each candidate object.For detected candidate objects, we model and account for this flux loss, as described in Section 3. Future DEEP data releases will investigate the use of a non-default coaddition statistics to mitigate our flux self-subtraction loss.We then subtract the template from each calexp to create difference images, in which we conduct our object search.

DETECTION AND CHARACTERIZATION USING SHIFT AND STACK
We use the Kernel Based Moving Object Detection (KBMOD, Whidden et al. 2019;Smotherman et al. 2021) pipeline to perform digital tracking on each long stare in the data set.KBMOD is a GPU-accelerated digital tracking pipeline that uses likelihood images, as defined below, to estimate the likelihood that there is a source following any given trajectory, represented by a set of velocities in which the images are stacked, using all pixels in the first image as an origin point.Our grid choice for the KBMOD trajectories as applied to DEEP is as follows: velocity ranges from 130 pixels/day to 400 pixels/day (1.4 arcsec/hr to 4.5 arcsec/hr, for comparison, at 40 au, a TNO's typical range of motion is 3 arcsec/hr) in 50 uniform steps.Angles are ±45 • of the direction of decreasing ecliptic longitude, divided uniformly into 30 steps, so that the maximum separation between neighboring trajectories would be less than about 2 PSF FWHM over a 4 hour time baseline.This grid choice exhausts the entirety of the Kuiper belt, and the majority of bound orbits between 30 and 80 au (see Bernardinelli & The DEEP Collaboration 2023, for a detailed presentation).
We define where σ j is the variance at pixel j, n(x j ) is the number of counts in pixel j at position x j in a given difference image, and T (x j − y) is the point spread function (PSF) centered at the true position y of a given point source.This means that Ψ, the likelihood image, represents the cross-correlation of the PSF and the data, a quantity equivalent to the convolution of the transpose of the PSF with the data (for a symmetric PSF, this reduces to the convolution of the PSF and the data).Ψ is simply the effective area of the PSF when weighted by the inverse variance of the counts in each pixel.These quantities lead to an optimal coaddition for point sources (Whidden et al. 2019).In this framework, sources in the image represented by ν = Ψ/ √ Ψ are χ2 distributed, and so, sources above a certain threshold ν thres are detections of this significance.
When coadding along images i in a trajectory, we have that: where ν is the S/N of a source in the coadd.Following the notation of Smotherman et al. (2021), we denote by LH the summed likelihood (equivalent to the S/N in the pre-convolved images) of a single trajectory.In other words, the summed likelihood is the S/N as given by equation 5 for a initial pixel (x, y) and trajectory.We refer the reader to Whidden et al. (2019) for more detailed description of the KBMOD procedure.
This formalism allows the core GPU algorithm in KBMOD to evaluate more than 10 10 candidate trajectories in a few minutes using consumer-grade GPUs (Whidden et al. 2019).For DEEP, a KBMOD search of a stack of O(100) images for a single CCD has a median runtime of 576s using an NVIDIA A40 GPU.This includes loading the data from disk, running the core GPU algorithm, filtering false positives and generating image stamps.Our total runtime for this search was about 282 hours.
In order to characterize our false positive rate, we also run "reverse" searches in the direction of increasing ecliptic longitude, which is opposite from the direction of a TNO due to the reflex motion of the Earth.These reverse searches will never yield a real TNO detection, but will lead to false positive detections due to the (trajectory independent) correlated noise (e.g.difference imaging artifacts and poorly subtracted sources) in the image series, and thus represent an accurate characterization of the false positive rate as a function of LH.We set the minimum sum likelihood threshold for a candidate trajectory to be considered to LH > 7.This limit was we defined by studying the number of candidates yielded in the reverse searches (where we expect that the majority of detections are false positives) vs LH: below LH = 7, the number of detections coming from the reverse searches increases by an order of magnitude.
Due to the number of false positives at very slow speeds (130 to 150 pixels/day, where the effect of the correlated noise sources is more pronounced) or in dense CCDs where the images were poorly differenced (when the number of trajectories with LH > 7 at velocities larger than 150 pixel/day 2 ), we increase the likelihood threshold from LH > 7 to LH > 10.This empirically determined limit significantly reduces the number of false positive trajectories that result from the search.
In order to further reduce the number of false positive trajectories in the search, we apply the same convolutional neural network (CNN) from Smotherman et al. (2021) on the resulting coadded stamps from the trajectories.This CNN discards candidates with a CNN-assigned probability of being from a "good" source of less than 50%.We note that in practice this threshold is merely a user-defined cutoff, rather than any kind of robust statistical likelihood -this choice ensures a good trade-off between true positives/false negatives and false positives.For the candidate trajectories that remain after the aforementioned cuts, we assign a randomized candidate id and generate images for human vetting.This is done with the combined results from both the forward and the reverse searches, so that we can estimate the false negative and false positive rate of the vetting procedure, as shown below.We show two examples of the images used for human vetting in Figure 2.This procedure helps further identify artifacts such as poorly substracted sources that otherwise produced an acceptable, PSF-like stack, as well as identify other cases where the CNN does not identify a false positive (e.g.unmasked diffraction spikes).This procedure rejected about half of the candidates coming from the CNN filtering.After human vetting, we further improve the quality of the determined trajectory and to characterize our errors in the starting and ending position, flux, and magnitude.We perform this analysis on all candidates that pass human vetting.This method uses the PSF determined for each image, as opposed to a constant-width gaussian PSF as used by KBMOD.The process for applying this "joint-fit" characterization is as follows.First, we load 31x31 pixel warped difference imaging and calexp cutouts, PSFs, magnitude zero points, stamp centers in R.A./dec and pixel space, and associated world coordinate systems (WCSs).Next, we rescale the stamps to a common magnitude zero point of m V R,0 = 31.After loading the data, we run a high-precision positional fit.We use the same formalism as KBMOD for flux and likelihood estimation, and minimize LH to find the topocentric best-fit trajectory, that is, the statistically optimal four dimensional solution of the starting and ending positions, which are each bounded to be no more than 10 pixels from the center of the corresponding stamp.By computing the 4 × 4 Hessian matrix of this LH minimization, we also determine the (astrometric) covariance matrix (defined as the inverse of the Hessian) of the starting and ending positions.These will be incorporated in our orbit fitting procedure in Section 4.4.
Finally, we measure the magnitude of each candidate object.When fitting magnitudes, we adjust the PSF of each image to model the effect of the "negative well" in the difference images.This "negative well" is caused by flux that was self-subtracted in the coaddition procedure.We model this effect by subtracting the relevant PSF scaled by the inverse of number of images at the location the candidate object for each image in the stack.Using this approach, we first measure the flux in each individual image in a single candidate trajectory.We use sigma-clippping to iteratively reject 3-σ flux outliers, and, for the cutouts that are not clipped, we then measure the variance-weighted flux and flux uncertainty using the same formalism as KBMOD.Next, we add an additional flux term in quadrature to the flux uncertainties such that the sigma-clipped standard deviation of the distribution of the residuals of the implanted and measured fluxes for the simulated objects corresponds to a unit Gaussian distribution.This helps to ensure we do not underestimate our uncertainties in flux.Future processings of the DEEP data will mitigate this problem by employing more refined flux measurement techniques in the presence of background sources (see Bernardinelli et al. 2023, for a recent example).We refer to candidates that are potential moving objects (as opposed to implanted fakes or false positives) as "real".First, to further ensure no fakes are in this sample, we run fake association on both starting and ending RA and Dec. From our list of candidate objects that pass human vetting, we remove any candidate object that is within 5 arcsec of the starting position or the ending position of a fake in the implanted fakes catalog.Only six additional candidates are associated with fakes when using either starting or ending position (rather than solely starting position).Regardless, we exclude these six additional candidates from the real distribution in the following analysis.
First, we again run the joint-fit, this time on the reverse search candidates that also passed human vetting.We then use these joint-fit results to assign a signal-to-noise ratio to each detection.We use this second run of the joint-fit to identify an appropriate S/N cutoff.Here, S/N is functionally identical to LH, except it uses the Φ and Ψ values from the joint-fit and is therefore more closely aligned to the actual S/N of the source.Namely, S/N = i Ψ i / i √ Φ i , where Ψ i and Φ i are generated using the negative well modelling described in Section 3.
In order to find the ideal S/N cutoff, above which we can classify a non-fake forward search candidate as real, we use our recovered fake distribution and our reverse search distribution.To characterize our false positive rate, we ran reverse searches on all stacks of CCD exposures in the long stares.This resulted in a total of 6973 reverse search candidates included in human vetting.Of these, we erroneously labeled 350 as good.Similarly, to characterize our vetting efficiency, we compare the distribution of fakes before and after human vetting.These distributions are shown in Figure 3 as a function of S/N .
Our goal is to find a S/N cutoff that keeps the greatest number of (forward search) fakes while rejecting the greatest number of reverse search candidates.We quantify this cutoff by iteratively increasing the S/N cutoff until the fractional decrease in the number of fakes kept is consistently greater than the fractional decrease in the number of reverse search candidates kept.Note that, as the samples of each group (real and false positives from the forward and reverse searches) are imbalanced, this cutoff does not apply equally to the two samples.Here, we define "consistently" to mean that the aforementioned comparison happens three times in a row.Above this cutoff, increasing the S/N threshold discards more fakes than it discards reverse search candidates.For the purpose of finding the S/N cutoff, we consider the range of S/N from 2 < S/N < 10, with an S/N step size of 0.05.As shown in Figure 3, we find this cutoff point to be S/N = 4.10.We will use this cutoff as one of several cuts in linking and orbit fitting.The number of all fakes that were included in vetting (blue) and the number of fakes that were labeled "good" in human vetting.The ideal case is that the orange histogram matches the blue histogram.S/N is bounded to 0 < S/N < 100 for visibility.Top right: The number of all reverse search candidates that were included in vetting (blue) and the number of reverse search candidates that were labeled "good" in human vetting (orange).The ideal case is that the orange histogram goes to zero.We bound to 0 < S/N < 100 for visibility.Bottom: A parametric plot of the number of reverse search candidates kept after human vetting (lower is better) versus the number of fakes kept after human vetting (higher is better) as a function of S/N .The blue dash-dot vertical line corresponds to the selected S/N cutoff value of 4.10.7) value of each distribution.
In Figure 4, we show the astrometric residuals and the photometric residuals for the joint-fit results with S/N > 4.1.We study the total astrometric residuals where α corresponds to RA at the start of the long stare, δ corresponds to Dec at the start of the long stare, d α is the residual in RA and d δ the residual in Dec and, for simplicity, this measurement ignores the directional dependence of the errors, that is, the direction of the (d α cos δ, d δ ) vector.In addition to the mean and standard deviation, we also compute σ G , an estimator for the underlying standard deviation of a Gaussian distribution in the presence of outliers.This is a more robust estimator of the "characteristic" scatter one might expect when selecting a random member of a data set (Ivezić et al. 2014;Smotherman et al. 2021).Given two quantile values from a data set, q i and q j , we have that σ G ≡ C i,j (q j − q i ) , where C depends on the choice of quantiles.While the choice of quantiles is rather arbitrary, the 25% and 75% quantiles (so C 25,75 ≈ 0.7413) offers a good balance between a large enough number of samples from the distribution without risking significant contamination from the outliers.The mean of the astrometric residuals is µ = 140.5 mas, the standard deviation is σ = 190.2 mas and σ G = 93.9mas.Given that the distribution of astrometric residuals is bounded to be positive, these are overestimating the astrometric quality of our data.We revisit this concept in Section 4.4.For the photometric residuals, the mean is µ = −0.0048mag, the standard deviation is σ = 0.1027 mag, with σ G = 0.0551 mag.The standard deviation implies that our typical source has S/N ≈ 10 (S/N ≈ 18 if using σ G ), a reasonable result considering that we are limiting our detections to sources with S/N ≳ 7 and that our fakes are uniformly distributed in magnitude space.

Linking and Orbit Fitting Procedures
Now that we have validated single-night recoveries by comparing to our population of fakes detected in a single-night long stare, we can link our single-night observations to fit detailed orbits.Furthermore, we can analyze the population of linked fakes to validate the real objects and real linked orbits that we detect.
We parameterize each single-night detection as a tracklet of two observations defined by the starting and ending on-sky positions, and, for simplicity, we take the largest diagonal value of the covariance matrix for each set of positions.Tracklets for which the joint-fit failed to return a covariance matrix (about 10% of the total) due to numerical instabilities in calculation of the Hessian matrix of the likelihood are given a default uncertainty in starting and ending RA and Dec of 270 mas (approximately 1 pixel).Note that, in Section 4.4, we revisit these uncertainties, so the overestimate here is not an issue.
We apply the linking methodology of Bernardinelli et al. (2020Bernardinelli et al. ( , 2022) ) 3 .We refer the reader to these papers for detailed presentations of the methods.The application of this technique to the temporally sparse Dark Energy Survey (DES) wide field data has proved to be > 99% efficient in an exhaustive population of synthetic objects subject to the DES geometric selection functions and the magnitude efficiency of individual exposures (Bernardinelli et al. 2022).The initial linking uses only the information of a single position per night and proceeds in bins of inverse geocentric distance γ ≡ 1/d.We limit our search to 20 < d < 2500 au, exhausting the range where we are able to recover TNOs in the digital tracking procedure, and enabling us to determine potential mislinkages at distances where we do not expect to see detections.
The process proceeds by clustering individual detections (represented by either end of a tracklet) in pairs in future observing nights, where the clustering radius defined by the limiting motion of a bound orbit at the chosen γ for the time span.In other words, given an individual detection and γ, detections from other exposures that lie inside a circular area with this limiting radius are clustered into pairs.Triplets are found by decomposing the motion of the pair into the parallactic and binding axes defined in Bernardinelli et al. (2022), representing deviations from the nominal search γ and by the limiting line of sight velocity given the first two detections.We proceed by fitting the orbits of all triplets in the distance bin, with a Gaussian prior on γ.Unlike the DES application of this methodology, we find pairs and triplets of detections across the multiple observing seasons of the survey.
The major difference between the two data sets that enables this generalization is that the DEEP candidate list has a low transient density and high purity, that is, most detections are real objects, while the DES transient catalog is dominated by asteroids and image artifacts (Bernardinelli et al. 2020).Each typical DEEP pointing has O(100) detections, with a typical false positive rate less than one false positive per real source, while the DES search had O(10 3 ) detections per exposure, with O(1) real TNOs per image.
After further detections are found in the other nights of data, we fit the candidate's orbit using the orbit fitting routines of Bernstein & Khushalani (2000), part of the orbitspp4 package.Each candidate orbit is propagated to all other exposure epochs, and we search for further detections within the projected on-sky 5σ error ellipse.Once an additional single night candidate is found, we iterate this process with the updated orbit, until all possible combinations have been exhausted.This leads us to an initial χ 2 with ν = 2n−6 degrees of freedom (n being the number of detections in the fit).
Once we have fit candidate orbits for our single-night detections, we apply cuts to all possible candidates in order to create our set of trustworthy multi-night linked orbits.We require that the orbital arcs span > 0.8 yr, and that the shortest arc remaining once removing either the first or the last night of the data to span > 0.5 yr.We also require that the object be detected in 4 or more nights.These choices were optimized so that the majority of our fakes are accepted by the linking process, while still minimizing the potential for false positives.
We require each real candidate orbit to be composed solely of real single-night detections and each fake candidate orbit to be composed solely of fake (implanted synthetic) single-night detections.We label candidate orbits with a mix of real and fake component tracklets as "bad", because no correct linkage can be composed of both real and fake component tracklets.
Next, we iteratively remove possible duplicate orbits.We iterate through the list of all real orbits, fake orbits, and bad orbits separately.For each orbit in the list, we find any other orbits that share at least one tracklet with the orbit.From this group, we save the orbit with the lowest χ 2 /ν and discard any other orbits that share tracklets with this orbit.If there are any remaining orbits, we again check if there are any shared tracklets and repeat the process if so.Then, we remove any candidate orbits in the bad population if that candidate orbit shares any single-night component tracklets with an orbit in the real or fake orbit populations.This process requires that each single-night long stare tracklet corresponds to only one orbit in any of the three (real, fake, or bad) populations.
At this point, we investigate which of our fakes are properly linked and which are mislinked.A mislinked fake orbit occurs when single-night tracklets from two different implanted fakes link to a single orbit, or when a fake connects to a real tracklet (that potentially belongs to a real TNO).Out of 250 candidate orbits, 20 are mislinked fakes and 230 are properly-linked fakes.We add these 20 mislinked fakes to the "bad" population.With these populations in mind, we can identify one final cut before investigating the real candidate orbits.
Given that we have recovered 230 of the total 244 possible candidates, we have a linking efficiency of 94%.We investigate these 14 synthetic objects that were not recovered in the linking process.Of these, ten had significant (> 2 ′′ ) astrometric residuals and, thus, are outliers of our residual distribution.Another two had poorly constrained positional uncertainties, with at least one sky coordinate having σ ≈ 1 ′′ , and so were rejected by the linking due to uncertain orbit fits.The final two objects did not have any obvious reason to be rejected by the linking algorithm, and so we believe those were true missed linkages.Accounting for these factors, the "true" linking efficiency is closer to the 99% reported by Bernardinelli et al. (2022), but any statistical characterization of our recovered objects must use the smaller 94% figure (plus any cut derived from the poor quality orbit fits of Section 4.4), as these problems also affect real TNOs.

Tracklet fitting
As we demonstrated in Section 4.2, our shift and stack "tracklets" have strong correlations not only between each pixel coordinate for any given time, but also between the positions at the beginning and the end of the long stares.If we do not account for these correlations, for example, by using a circular error matrix where the radius is the semi-major axis of the error ellipse or only include the correlations between each component in a given exposure, this would make our orbit fits not properly reflect our astrometric capabilities.This would potentially lead to uncertain or incorrect dynamical classifications, and it also makes any "fit quality" tests (i.e.selecting on χ 2 /ν) complicated, as we do not know if we can trust our uncertainties.
To address these problems, we expand upon the methodology of Bernstein & Khushalani (2000) to include the full correlation matrix obtained by the joint-fit.Our goal is to fit the Bernstein & Khushalani (2000) parameter set α ≡ {α ≡ x/z, β ≡ y/z, γ ≡ 1/z, α ≡ ẋ/z, β ≡ ẏ/z, γ ≡ ż/z}, where the z coordinate corresponds to the line-of-sight distance d between target and observer.
We transform the starting and ending sky positions to tangent plane coordinates {θ x , θ y , θ ′ x , θ ′ y }.The difference between each tracklet for long stare µ with starting observing time t µ and ending time t ′ µ , and the model coordinates { θx (t|α), θy (t|α)} is Thus, our goal is to find α that minimizes where C µ is the covariance matrix of the tracklet.We add the covariance determined by the joint fit procedure to the identity matrix scaled by the excess astrometric uncertainty determined in Section 4.2, σ 2 ast = (93.9mas/ √ 2) 2 .Note that the √ 2 is needed as this value was determined for both RA and Dec simultaneously, while here we need an error per coordinate.
In order to maintain the numerical stability of the fit, we first fit the orbit with γ = 0 and circular errors corresponding to the semi-major axis of each error ellipse.After this fit converges, we re-fit the full orbit with the χ 2 in Equation 9 using the same routines as Bernstein & Khushalani (2000) and modifications presented in Bernardinelli et al. (2020), namely the weak bound orbit prior in their Equation (15).We do not include their prior on γ, as with 4+ temporally sparse tracklets the solution is well-behaved and can reach its global minimum easily.The covariance matrices for α are calculated using the analytical derivatives of the model, as in Bernstein & Khushalani (2000).The reported state vectors are ICRS aligned with the Solar System barycenter as the origin, with a common epoch for all objects.In addition to the Sun, the four giant planets are treated as active particles for the orbit fit, and their positions are derived from the DE-440 ephemerides (Park et al. 2021).The derived orbital elements are barycentric and use the total mass of the eight planetary systems in addition to the Solar mass for the conversions from state vectors.Now, as we've determined a new, more reliable χ 2 /ν for our objects, we can apply one additional cut to select our high-confidence orbits.Figure 5 shows the distributions for the our three subsets of linkages (real, fake and bad linkages).We choose a threshold that accepts 95% of the correctly linked fakes, χ 2 /ν < 4.5, leaving us with 218/244 = 89.3%succesful linkages.Of our 38 bad linkages, only 3 have χ 2 /ν ≤ 10 (4.1, 5.9 and 6.6, from smallest to largest), so our contamination from real mislinkages, if it exists, is minimal.Only two of our real objects have 4.0 ≤ χ 2 /ν ≤ 4.5 and, beyond this threshold, the distribution increases slightly (at 5 objects with χ 2 /ν ≈ 5) before dropping again.
We also note that our χ 2 /ν distributions (for both real and fake objects) peak very close to χ 2 /ν = 1, indicating that the joint-fit procedure, with this additional 93 mas scatter per coordinate pair, is fully capturing our astrometric precision.The right panel of Figure 5 shows the fractional error in semi-major axis a for both the real and fake linkages.A typical orbit of a real object has σ a /a ≈ 10 −4 , an orbit quality that is comparable to more "traditional" surveys like OSSOS (Bannister et al. 2018) and DES (Bernardinelli et al. 2022) that discover objects in individual images.We note that the differences in both distributions are direct effects of the differences between the magnitude and orbital distributions of the real and fake populations: the uniform magnitude distribution of the fakes leads to a higher fraction of bright, better measured objects when compared to the approximately power law distribution of the real objects (on average fainter, thus with worse astrometric measurements).However, as the real objects (as will be seen below) are primarily classical KBOs with relatively low eccentricities, their semi-major axes are better constrained than the more eccentric population of fake objects, which reproduces the well known correlation between a and e in trans-Neptunian orbit recovery (see, e.g.Bernstein & Khushalani 2000). .Left: Histograms of χ 2 /ν determined with the full covariance matrices for the subsets of linkages corresponding to real detections (blue), fake detections (red), and the bad linkages, that is, those that we know are mislinked (orange).The vertical line indicates the χ 2 /ν = 4.5 threshold.Right: Histogram of the fractional error in semi-major axis (σa/a) for the real (blue) and fake (red) orbits with χ 2 /ν < 4.5.We note that the different structure in the distributions is not an inherent problem, but stems from the very different populations of orbits that belong to these two samples.

Object sample
After applying the aforementioned cuts, we have a catalog of 110 detected real objects with multi-year linkages.We use skybot (Berthier et al. 2006) to identify if there are known objects in the DEEP B1 data.We associate all single-night candidates input into linking with KBOs reported with skybot.If both the starting RA and the starting Dec of the detection are within 20 arcsec of the RA and Dec predicted by skybot in multiple DEEP pointings, then we consider a candidate associated with a known KBO.This leaves us with 4 known KBOs (2003QL 91 , 2002 PV   .The median topocentric discovery distance for all correctly-linked fake orbits and all 110 real orbits that pass all orbital cuts.The median is taken from the distribution of topocentric distances for each long stare (and thus each tracklet) in the linked orbit.
2002 PX 170 , and 1999 RX 215 ), all of which were recovered, further increasing our confidence in our search.The orbital elements, magnitudes, dynamical classifications and MPC identifiers for all of our 110 objects are presented in Table 2.
In Figure 6, we show the median topocentric discovery distance for the 230 correctly-linked fake orbits and the 110 real orbits which pass all orbital cuts.This distance is defined as the median of the topocentric distances reported by the orbit fit at each long stare for which the object was observed.As expected, our results for the fake objects are consistent with the distances for the fake objects detected in single long stares, shown in Figure 1 of Bernardinelli & The DEEP Collaboration (2023).The distribution of discovery distances for the fake orbits encompasses the distribution of the discovery distances of the 110 real orbits, as one would expect for correctly-linked real orbits.
In order to determine dynamical classifications for our objects, we follow the classification scheme of Gladman et al. (2008) and Khain et al. (2020) using the implementation of Bernardinelli et al. (2022).We create 20 clones of each object, and then use Rebound (Rein & Liu 2012) with the WHFast symplectic integrator (Wisdom & Holman 1991;Rein & Tamayo 2015) to integrate these for 10 Myr from the starting time t 0 .Clones that librate around the resonance argument for more than 90% of the simulation time are considered resonant.Clones where the semi-major axis varies by more than 3.75% from a(t 0 ) are considered scattering.Clones where a(t 0 ) < a N are considered inner Centaurs, where a N is the semi-major axis of Neptune.Clones with e(t 0 ) > 0.24 are considered detached.Remaining clones are considered classicals.Objects are classified based on the classification of the majority of their clones (> 50%).Objects where more than 80% of the clones show the same mean-motion resonance (MMR) for more than 90% of the time are classified as securely resonant, and those with 50 − 80% or more of the clones showing resonant behavior are classified as insecurely resonant.Objects where no single classification reaches a majority are considered insecurely classified; these are presented in Table 2.All 110 TNOs and their classifications are shown in Figure 7.
We identify six resonant TNOs, in the 3:2, 7:4, 2:1 resonances, with more objects in the 3:2.This is in line with previous surveys (e.g.Bernardinelli et al. 2022), but we will defer any statistical claims to future releases, when fields that are expected to contain more resonant objects are analyzed.The B1 fields are close to Neptune's position in the sky, so these fields are generally far from the 3:2, 7:4, and 2:1 resonant populations, which explains the proportionatelysmaller number of resonant objects compared to (for example) Bannister et al. (2018) or Bernardinelli et al. (2022).
We identify 96 classical TNOs, the majority of our sample, unsurprisingly, given our survey geometry.66 of these objects are dynamically "cold" (i ≤ 5 deg), another feature expected from the survey geometry, as our fields are in the invariable plane.The four recovered known objects are in this population.One of our "classical" objects has a ≈ 56 au, close to the 5:2 MMR, adding to the population of high-q non-resonant objects near these MMRs (see Sheppard et al. 2016;Lawler et al. 2019;Bernardinelli et al. 2022).
We also identify six detached objects.One of these detached objects (a = 77.53au, e = 0.52, i = 20.69• ) was detected with a median topocentric discovery distance of 76.23 au, being our most distant object discovery.Finally, we also have two scattering objects.With the detection of these scattering objects, we demonstrate that we are able to detect and characterize the four main dynamical populations of TNOs (classical, scattering, detached, and resonant).
Finally, we investigate the observed H magnitude distributions.We compare the residuals of the simulated values for median H magnitude and the median of the observed values for H magnitude.Note that, as the fakes were generated without phase coefficients, we likewise do not include these here.The median of the residuals in H is −0.001, indicating that there is no observable systematic offset.The standard deviation is σ = 0.086 (σ G = 0.031), also indicating that we are correctly recovering our absolute magnitudes.
As our real objects have between 4 and 6 shift and stack detections with similarly low phase angles (as expected from the survey design), and are in their majority cold Classicals, a population that has higher variability than the other TNO subpopulations (Strauss et al. 2023;Bernardinelli et al. 2023), we elect to report the mean absolute magnitude H over all detections; any choice of phase coefficient and typical light curve amplitude would lead to systematic offsets in H for an undetermined subset of the population (for which these parameters are not representative).We show the distribution of these values in Figure 8.Our smallest object (that is, the one with the largest H) has H = 10.2, being one of the smallest objects with q > 30 au with multi-year arcs to date.

Structure of the classical Kuiper belt
We investigate the distribution of semi-major axis of our classical TNOs and their relationship to models of the Kuiper belt.We compare our data to the L7 model of the Canada-France Ecliptic Plane Survey (CFEPS, Petit et al. 2011), where the classical Kuiper belt is divided into a dynamically excited, "hot" population, and two dynamically cold ones: the "kernel", concentrated in the 43.8 < a < 44.8 au range, and the "stirred" population covering the 40 < a < 48 au range.While this model is known to not reproduce the inclination distribution of these objects (Petit et al. 2017;Bernardinelli et al. 2022), this model reproduces well their semi-major axis distribution (Bannister et al. 2016;Bernardinelli et al. 2022).Since we are probing a different size range than the one in the published CFEPS-L7 model, we have re-sampled the model and changed the size distribution of the cold and hot classical populations to those of Fraser et al. (2014), to provide a sufficient sample of objects in our size range.
We have used our survey simulator (Bernardinelli & The DEEP Collaboration 2023) to produce a population of objects based on this model when subject to our survey biases, and compare their a distribution with our data.The results are shown in Figure 9.We do not see a sharp increase in the number of objects at 43.8 < a < 44.8 au predicted by the narrow Kuiper belt kernel in the CFEPS-L7 model.The model distribution increases sharply in the 42.6 ≲ a ≲ 44.6 range, and systematically over-predicts objects with a ≳ 46 au when compared to the observed data.We compute a Kolmogorov-Smirnov (KS) test (Press et al. 2007) between these two distributions, and obtain a p-value ≤ 10 −5 that the data and the model are drawn from the same distribution, confirming the visual impression that these distributions do not agree.
An implication here is that the kernel either must be wider than the 1 au currently in the CFEPS-L7 model, since our object density is ∼constant (within 1σ) in the 42.6 au to 44.6 au range, or that the cold Classical population must have a different structure than that implied by the stirred+kernel model.A more complicated alternative would be for such differences to be as a function of size: the CFEPS-L7 model was derived with H < 8 objects, while here the majority of our objects are fainter than such H.However, by restricting ourselves to H < 8 in both the model and data, the K-S test leads to a p-value = 0.0047, which still indicates a discrepancy.Understanding these differences will have implications for dynamical models that seek to reproduce the Kuiper belt.In particular, the formation of the kernel requires, for example, that Neptune "jumps" during its migration (Nesvorný 2015).We defer more detailed tests of these dynamical models to future releases of the DEEP data, when our object counts are expected to increase by a factor of 10 (Trilling & The DEEP Collaboration 2023).
We also note that Napier et al. (2023) used the DEEP B1 data to infer the absolute magnitude distribution for the cold Classical population.Further analysis of DEEP data will, then, enable more detailed tests of the joint radial and absolute magnitude distribution of the Kuiper belt.

SUMMARY
We have demonstrated that we are able to systematically link and fit high-quality multi-year orbits to sources identified with digital tracking.We have recovered 110 TNOs (4 already known), including objects as faint as m V R ≈ 26.2 and H V R ≈ 10.2.This means that we are able to probe the dynamical structure of the trans-Neptunian region to systematically smaller sizes than any other TNO survey with more than 100 objects to date.The population is, as expected from the survey design and geometry, dominated by classical Kuiper belt objects in the 40 < a < 48 au range, but we have also recovered a detached object (a ≈ 78 au) with H V R ≈ 6.1, and a few TNOs in the other dynamical classes.We used our discovered Kuiper objects and our survey simulator to compare the distribution of our classical objects with the CFEPS-L7 model, and conclude that our data is inconsistent with that model.Future investigations of the Kuiper belt should lead to improved model constraints of this region.
As this is the first multi-year arc release of the DEEP survey, a number of improvements are expected, including a more complete astrometric characterization and an additional year of observations for the B1 field.We also have more fields in which this analysis will be conducted.We can expect that the number of discovered TNOs by the DEEP survey will significantly increase, and this increase in statistical power will make these data set even more valuable for testing models of the trans-Neptunian region.

Figure 1 .
Figure 1.On-sky positions of the B1 field.The pointings included here are the first observation of each field taken in 2021.Each field is centered on the RA and Dec of the first exposure in the long stare.The solid line shows the invariable plane, and the dashed line shows the ecliptic.

Figure 2 .
Figure2.Two examples of objects recovered with KBMOD, the first an injected TNO with magnitude mV R = 24.9 and a barycentric distance of 60 au, the second a real cold Classical TNO with mV R = 26.0.These images are used for human vetting.The title includes information from the KBMOD search, including LH, flux, starting x and y pixel location, and speed.The two postage stamps in the upper left are sum and median coadds, as well as an approximate S/N estimated from these postage stamps.The graph in the upper right shows an approximate flux lightcurve.The postage stamps included below are ther first 20 our of 100 cutouts from individual images in the long stare.
. S/N clipping of candidate single-night detections

Figure 3 .
Figure 3.A visualization of the signal-to-noise ratio cutoff used to determine if a candidate should be considered real.Top left:The number of all fakes that were included in vetting (blue) and the number of fakes that were labeled "good" in human vetting.The ideal case is that the orange histogram matches the blue histogram.S/N is bounded to 0 < S/N < 100 for visibility.Top right: The number of all reverse search candidates that were included in vetting (blue) and the number of reverse search candidates that were labeled "good" in human vetting (orange).The ideal case is that the orange histogram goes to zero.We bound to 0 < S/N < 100 for visibility.Bottom: A parametric plot of the number of reverse search candidates kept after human vetting (lower is better) versus the number of fakes kept after human vetting (higher is better) as a function of S/N .The blue dash-dot vertical line corresponds to the selected S/N cutoff value of 4.10.

Figure 4 .
Figure 4. Astrometric (left) and photometric (right) residuals for single-night detected fakes with S/N > 4.1.These results come from the output of the joint-fit process described in Section 3. The title of each subplot shows the mean µ, standard deviation σ, and σG (see Equation7) value of each distribution.
Figure5.Left: Histograms of χ 2 /ν determined with the full covariance matrices for the subsets of linkages corresponding to real detections (blue), fake detections (red), and the bad linkages, that is, those that we know are mislinked (orange).The vertical line indicates the χ 2 /ν = 4.5 threshold.Right: Histogram of the fractional error in semi-major axis (σa/a) for the real (blue) and fake (red) orbits with χ 2 /ν < 4.5.We note that the different structure in the distributions is not an inherent problem, but stems from the very different populations of orbits that belong to these two samples.
Figure6.The median topocentric discovery distance for all correctly-linked fake orbits and all 110 real orbits that pass all orbital cuts.The median is taken from the distribution of topocentric distances for each long stare (and thus each tracklet) in the linked orbit.

Figure 7 .
Figure7.Eccentricity versus semi-major axis (top) and inclination versus semi-major axis (bottom) for the 110 objects from the DEEP B1 fields.Color denotes dynamical classifications.Vertical dotted lines correspond to the approximate location of select Neptunian p:q mean motion resonance.Curved dotted lines (top panel only) correspond to lines of constant pericenter.

Figure 8 .
Figure8.The inverse variance weighted mean H magnitudes of the 110 detected real objects.This is computed using an H-G model with an assumed value of G = 0.15.

Figure 9 .
Figure 9. Cumulative (top) and differential (bottom) histograms of semi-major axis distribution of the predicted CFEPS detections (orange) given our survey biases and our recovered classical objects (blue), limited to the 40 < a < 48 au range.The differential histogram was normalized to the number of real objects, and we also include Poisson error bars in the counts of the real objects.

Table 2 .
170 , Orbital properties of the TNOs from DEEP B1 Note-The table is provided in its entirety in machine-readable format (FITS), here we are only describing the columns included in the data release.All elements correspond to Julian epoch 2019-08-30, and are barycentric and ecliptic-aligned.