ABSTRACT
Digital co-addition of astronomical images is a common technique for increasing signal-to-noise ratio (S/N) and image depth. A modification of this simple technique has been applied to the detection of minor bodies in the solar system: first stationary objects are removed through the subtraction of a high-S/N template image, then the sky motion of the solar system bodies of interest is predicted and compensated for by shifting pixels in software prior to the co-addition step. This "shift-and-stack" approach has been applied with great success in directed surveys for minor solar system bodies. In these surveys, the shifts have been parameterized in a variety of ways. However, these parameterizations have not been optimized and in most cases cannot be effectively applied to data sets with long observation arcs due to objects' real trajectories diverging from linear tracks on the sky. This paper presents two novel probabilistic approaches for determining a near-optimum set of shift vectors to apply to any image set given a desired region of orbital space to search. The first method is designed for short observational arcs, and the second for observational arcs long enough to require nonlinear shift vectors. Using these techniques and other optimizations, we derive optimized grids for previous surveys that have used "shift-and-stack" approaches to illustrate the improvements that can be made with our method, and at the same time derive new limits on the range of orbital parameters these surveys searched. We conclude with a simulation of a future application for this approach with LSST, and show that combining multiple nights of data from such next-generation facilities is within the realm of computational feasibility.
Export citation and abstract BibTeX RIS
1. INTRODUCTION
Detecting faint objects in the solar system is a technically challenging problem. The approaches taken in the past can be broken into two basic categories: interframe source detection and linking, and the "shift-and-stack" approach. Interframe detection, where transient sources are detected in individual exposures and motion is linked between multiple exposures, is widely employed and favored in planned proposed surveys like Pan-STARRS (Denneau et al. 2007) and Large Synoptic Survey Telescope (LSST; Axelrod et al. 2009) due to its simplicity and robustness. This method, however, does not exploit imaging data to the fullest extent possible, because detections are subject to the completeness limit of individual frames.
Since solar system objects are moving across the sky, their images trail and the noise contribution from the sky background increases. Trailing losses create an upper limit to individual frame exposure times for efficient detection of solar system objects. It is advantageous to combine multiple exposures, each short enough that trailing effects are negligible. Co-addition is routinely done to improve the S/N of stationary sources and to remove cosmic ray events, but in order to apply it to the detection of moving objects, several other steps must be performed. First, contamination from stationary sources can be removed by subtracting a high-S/N template image from each image in the set. Next, the sky motion of the objects of interest must be predicted and this motion compensated for by shifting each image's pixels in software. When the images are then combined, only the flux from sources that moved at the predicted rates will be constructively added.
This method first appeared in literature in Tyson et al. (1992), though details of the application were limited. An early detailed description of this technique was presented by Cochran et al. (1995). A set of images from the Hubble Space Telescope (HST) Wide Field Planetary Camera 2 (WFPC2) were combined to allow statistical detection of very faint trans-Neptunian objects (TNOs). The sky motion of the sources was parameterized as angular rates (
) and angles on the sky (ϕ). Gladman et al. (1997) used the same technique in order to combine a series of images from the Palomar 5 m telescope in the hope of constraining the size distribution of Kuiper Belt Objects (KBOs) by discovering objects with smaller radii than had been detected in most previous wide-area surveys that employed interframe detection. This survey was followed by further searches from Keck (Luu & Jewitt 1998; Chaing & Brown 1999), CTIO (Allen et al. 2001 and 2002; Fraser et al. 2008), Palomar (Gladman et al. 1998), VLT and CFHT (Gladman et al. 2001), which employed the same technique to varying degrees of success. Recently, Fraser & Kavelaars (2009) and Fuentes et al. (2009) both used similar techniques to search Subaru data for KBOs as faint as R ∼ 27.
All of these surveys had observational arcs less than 2 days in length. Over periods this short, the motion of outer solar system sources can be well approximated (with respect to ground-based seeing) by the linear trajectories that the method these surveys implemented implicitly assumes. For longer arcs or higher-resolution data, however, the approximation of a linear trajectory is no longer adequate. Bernstein et al. (2004) employed a more advanced approach, which they called "digital tracking," to extremely high-resolution HST Advanced Camera for Surveys (ACS) data. Even in arcs on the order of a day in length, nonlinear components of motion were nonnegligible relative to the ACS psf. In order to stack frames such that distant moving sources add constructively in this regime, it becomes necessary to generate nonlinear shift vectors via full ephemerides from orbits of interest. They stacked arcs up to 24 hr in length using shift vectors parameterized as a grid over three parameters: two linear components of motion and heliocentric distance d. By sampling densely over this grid, they were able to recover orbits with 25 AU ≤ d ≤ ∞ and i < 45°, and detect objects as faint as R ∼ 28.5. However, this dense sampling, somewhat unconstrained discovery volume, and extremely small psf required ∼7 × 105 shift vectors, which translated into ∼1014 pixels to be searched. Bernstein et al. found that stacking more than 24 hr of their data became computationally prohibitive.
In the age of ground-based gigapixel CCD mosaic imagers, how can data obtained over long arcs be efficiently exploited to search for faint sources? In the interest of exploring massive data sets with long observation arcs to as deep a limit as possible, we have created a simple probabilistic method for generating an optimized set of shift vectors for any imaging survey and any interesting volume of orbital element space. By generating shift vectors in a probabilistic way, exploring an explicitly demarcated region of orbital parameter space becomes straightforward. The targeted region of orbital parameter space is searched with uniform sensitivity, simplifying the accurate characterization of any detected sample of objects. This method also avoids searching nonphysical or uninteresting parameter space, thereby maximizing scientific return for minimum computational cost.
In § 2.2, we define the maximum tracking error, the basic parameter that defines the density of a grid. Section 2 describes previous analytical methods for search-grid parameterization, our probabilistic approach, and a grid-tree structure for optimizing multiple-night arcs. Section 3 compares the efficiency of our method to those used in previous surveys, and in § 4 we show an example of application of this method to the kind of data that may be acquired by the LSST.
2. GENERATING SEARCH GRIDS
2.1. Estimates of On-Sky Motion
Previous on-ecliptic surveys (e.g., Fraser & Kavelaars 2009) have used simple analytical estimates for the apparent on-sky rates of motion for distant solar system objects in order to generate shift rates. For an on-ecliptic observation of a field β degrees from opposition, the range of on-sky angular rates of motion 
and angle from the ecliptic ϕ for a distant object at heliocentric distance d, and geocentric distance Δ with inclination i and eccentricity e, can be approximated by the vector addition of the reflex motion of the object due to Earth's motion (inversely proportional to Δ) and the intrinsic on-sky orbital motion of the object by Kepler's law (proportional to d-3/2)
where 
, representing the fraction of the mean angular velocity of an orbit at pericenter (+e) or apocenter (-e), compared to a circular orbit at heliocentric distance d.
This motion can also be parameterized as two components of angular rates of motion, one parallel and one perpendicular to the ecliptic. Rearranging equation (1) into parallel and perpendicular components, we find
The 
and ϕ parameterization is convenient for visualization purposes, but in § 2.2 we show that the 
and 
parameterization is more appropriate for generating a search grid.
The angular rate of motion is lowest (for a given d and Δ ≃ d) for i = 0° orbits at pericenter, since vd is at its highest and the parallax and orbital motion vectors are antialigned.1 The highest angular rate of motion at opposition is reached in two cases, depending on the maximum inclination and eccentricity considered (imax and emax):
where
and 
, 
.
Given this approximation, the highest ϕ any orbit can achieve is a strong function of d. If we consider an orbit with e → 1 at pericenter with i = 90°, and approximate Δ ∼ d, then by equation (2) it can be shown that
2.2. Grid Geometry and Maximum Tracking Error 
max
In using digital tracking, the motions of real sources are approximated by a grid or some distribution of shift- or motion-vectors. In order to quantize the effectiveness of this strategy, we define a maximum tracking error to be the largest possible error tolerated between any real object's motion and the nearest predicted motion vector. To estimate its effect, we compare the final signal-to-noise ratio S/N' for a faint circular source with flux f and area A0 before and after a linear tracking error is added in an image with exposure time t. The linear tracking error adds an additional rectangular region to the area of the source, with width 2R (where R is the aperture radius applied to the initial circular source) and length . This rectangle has area 2R, and the total area of the blurred source is now
The final signal-to-noise S/N' is given by
So, for a search grid with maximum tracking error max, the worst possible S/N degradation is a factor of 
. Given a Gaussian psf with a full-width at half-maximum of Γ, the radius at which the S/N of a sky-dominated source is maximized is R ≃ 0.68Γ. If we define this as our initial aperture radius, then the maximum allowable tracking error given a desired limit on F becomes
Previous surveys have used widely varying values for max with respect to the typical seeing of that survey, ranging from ∼0.1Γ–3Γ. Surveys that searched data by eye were limited in the number of shift vectors they could apply, and leveraged the eye's robust noise discrimination to compensate for the loss in S/N due to using large max (e.g., Fraser & Kavelaars 2009). Surveys using automated detection pipelines have decreased max to limit their S/N loss, and leveraged massive computational resources to compensate for the increased number of required shift vectors (e.g., Fuentes et al. 2009).
In order to evaluate the maximum tracking error of any search grid of shift vectors applied to a given data set, we determine the distribution of shifts from the initial image to the final image, where the shift for orbit k in image i is
The maximum separation any real orbit's motion from one of these offsets determines the maximum tracking error of the search grid. After time t has elapsed, the on-sky angular separation Ω of two objects starting from the same initial position with identical angular rates of motion 
and on-sky angles of motion separated by small angle dϕ is given by
The maximum tracking error between a search grid of shift vectors and a real distribution of orbits depends on the geometry of the search grid. Consider the grid spacing of a geometrically regular search grid in the plane defined by the total changes in RA and DEC made by moving sources during the observations, which we will refer to as the (dα, dδ) plane. If the points of the search grid are arrayed at the vertices of identical adjoining geometric cells on the (dα, dδ) plane, the point most distant from any point on the search grid is the point at the center of each cell. The distance from any vertex of a cell to the center of that cell defines the maximum tracking error of that search grid.
For points arrayed on the vertices of a grid of adjoining rectangles with sidelengths dA and dB, the maximum tracking error is 
. In the case where dA = dB, this becomes 
. This geometry well approximates the search grids used by most previous surveys. The left panel of Figure 1 illustrates a grid with this geometry.
Fig. 1.— Comparison of grid geometries. Left panel: Packing of a square lattice, with scales illustrated. Right panel: packing of a hexagonal lattice, with scales illustrated.
If a fixed (
, dϕ) grid spacing is adopted such that a maximum tracking error max is satisfied after time t for objects with angular rate of motion 
, objects with 
will see increasing tracking errors. To illustrate, consider a roughly rectangular grid where 
and 
. Inserting these values into our definition of max for a rectangular grid, we see that 
. A fixed (
, dϕ) grid parameterization oversamples motions with low 
, and undersamples motions with high 
.
The over- and undersampling problem can be remedied by parameterizing the sky motion as two components of angular rates of motion, one parallel to, and one perpendicular to the ecliptic (
and 
). By adopting a fixed spacing in (
, 
) along with angles from the ecliptic (ϕmin, ϕmax) between which motions are confined, a uniform max can be maintained for all angular rates of motion. This is the approach that was adopted by Fuentes et al. (2009).
The packing efficiency of a search grid can be improved over the regular rectangular case. If the area of (dα, dδ) to be searched by a grid is much larger than the unit of area searched by an individual grid point such that we can largely ignore boundary conditions, the proof by Gauss (1831) regarding the most efficient regular packing of circles on a plane holds. In this case, defining the points of the search grid as the vertices of a grid of adjoining equilateral triangles will produce the highest packing efficiency. Considering a grid of equilateral triangles with sidelength dA, the maximum tracking error is 
(where the effective dB spacing becomes 3max/2). The right panel of Figure 1 illustrates a grid with this geometry. Both our short- and long-arc solutions described in the next section produce a grid with approximately this geometry, and typically require ≥20% fewer grid points than a similarly defined rectangular grid with the same max.
Regardless of the grid geometry, each cell is basically covering a small area of the final (dα, dδ) plane. The final area of the (dα, dδ) plane covered by real objects is approximately ∝ t2 (equations [3] and [4], neglecting acceleration), while the area covered by a single grid point is 
. The total number of grid points N required to search a given population scales as 
.
2.3. Probabilistic Methods for Generating Shift Vectors
The methods for generating shift vectors described here are cumbersome to analytically generalize over any region of the sky, arbitrary regions of orbital parameter space, and arbitrarily long observational baselines as they still suffer from the problem of nonlinear components of motion that become significant as t grows. In this regime it becomes necessary to add additional dimensions to the search grid (
and 
, for example), and selecting an optimized set of grid-points becomes a challenge. Instead of pursuing an analytical solution, in the following section we outline a simple probabilistic approach for defining an optimum set of shift vectors for any region of the sky and over any length of observational baseline at any solar elongation.
Instead of approximating the minimum and maximum angular rates of motion and angle on the sky for an ad hoc subset of the orbits of interest, our method takes the following approach: First, generate a very large sample of synthetic orbits that fill (in a prescribed way) the orbital parameter space of interest. Second, use ephemeris software to determine the (α, δ) shift between frames for every synthetic orbit in every imaged epoch. Finally, search for the optimum (smallest) subset of orbits whose motions accurately represent the motions of the entire set of synthetic orbits within a maximum tracking error tolerance. The set of (α, δ) shifts (translated into [x, y) pixel shifts] from this representative set of orbits is our optimum set of shift vectors. A more detailed description of the method follows.
2.3.1. Initial Sample Generation
To create our initial sample, we use a randomly generated set of synthetic orbits that fill a region of orbital parameter space of interest such that the orbital ephemerides intersect the imaged field during the times of observation. The approach to filling orbital parameter space is subject to some consideration, as any orbital-space biases may translate into a skewed on-sky motion distribution, resulting in a biased selection of "optimum" orbits. Over small ranges of d, a uniform sampling is adequate, but for larger ranges we sample uniformly with respect to d-c, where the index c = 2 is appropriate to generate a uniform distribution in 
for observations near opposition. Once we have selected a given d, we select a and e from a uniform distribution between the minimum and maximum values set for those parameters, in the following order: first select e, then a such that it falls within its limits and also a≥d(1 - e).
In cases where a population has a particular orbital parameter that is well-constrained, it is preferable to generate this parameter first (e.g., a for a mean-motion resonance [mmr]), and then force the other parameters onto as uniform a grid as possible. In the case of populations in mean-motion resonances, a is pinned at the resonant value ammr, then d is selected from a uniform distribution between ammr(1 - emax) and ammr(1 + emax). Finally, e is selected such that ammr(1 - e) ≤ d ≤ ammr(1 + e). In either case, a, e, and d define two possible values for M, which are given equal probability.
Sampling inclination from the distribution of an isotropic sphere (probability of i ∝ sin(i)) may seem appropriate to avoid biasing the on-sky motion distribution to that of low-inclination orbits. However, because the component of motion perpendicular to the ecliptic is roughly ∝ sin(i), and because the small area of any imaged field with respect to the entire sky places stronger limits on the phase-space volume available to high-inclination objects compared to low-inclination objects, we contend that is preferable to sample i from a uniform distribution. The minimum inclination orbit possible to be observed is 
, where lmin is the minimum ecliptic latitude of the field. In the case where retrograde orbits are of interest, duplicate the set of inclinations with ir = 180 - i.
The remaining orbital elements to be generated are the node and the argument of perihelion. These can be rotated to force the object to fall on the imaged field during the dates of observation, and should fill the parameter space set by these constraints smoothly. M and the argument of perihelion are coupled for objects in mean-motion resonances, but for the purposes presented here it is convenient to simplify the treatment by treating the two as independent.
Depending on the size of the orbital region of interest and the length of the observational arc, the number of synthetic orbits required to produce a statistically adequate sample varies. The basic requirement is that the density of points in the final (dα, dδ) plane is such that the max spacing between any two nearest neighbors in the initial sample is ≪ max, ensuring that there are no empty bins in (dα, dδ) given circular bins with radius max. In trials, we found that this required the initial sample size to range from several thousand to hundreds of thousands of synthetic orbits.
Once a set of synthetic orbits has been created, the shift vectors for each orbit can be determined. We define the shift vector for orbit k in image i to be a point in the (dαi, dδi) as given by vk,i in equation (8). The list of shift vectors [vk,0...n] is stored as a list linked to the orbit k for which they were generated. We then compare each list of shift vectors to determine the optimum set, which is the smallest subset of vectors having motions that accurately represent the whole set within some spatial tolerance in the image plane. We set this tolerance by specifying that all synthetic orbits in the initial sample must be matched to the motions of at least one orbit in the optimum sample to within our maximum tracking error tolerance max in every imaged epoch. In other words, every synthetic orbit k must be matched to at least one optimum orbit h in every image i by |vk,i - vh,i| ≤ max.
2.3.2. Short-arc Solution
For observational baselines over which any nonlinear component of motion is negligible compared to e, the largest offset between any orbit and the nearest predicted motion vector will always occur in the final image n. To generate shift vectors in this regime, we use the following approach:
- 1.Propagate all synthetic orbits in the initial sample forward to the final image (n) plane and record the final shift vector vk,n for each orbit k. Find the median of the distribution of shift vectors and set the corresponding point in (dαn, dδn) to be the center of a geometrically-optimum triangular grid with sidelength
that extends over a region of (αn, dδn) significantly larger than that occupied by any synthetic orbits in the initial sample. - 2.From this initial large grid, select only those grid points that have one or more synthetic orbits whose final shift vector lies within of it in (dαn, dδn). Record a list of the orbits that are matched to each grid point.
- 3.After the initial grid has been cut down to just those grid points that are matched to at least one synthetic orbit, find the smallest subset that are matched to unique synthetic orbits: First, select the grid point that is matched to the largest number of synthetic orbits, set it and its list of orbits aside, then remove the references to each synthetic orbit recovered by that grid point from all other grid points. Repeat this process for the remaining grid points, until all the original synthetic orbits are accounted for in the set-aside list. All set-aside grid points are now identified as the optimum set for recovering the initial sample of synthetic orbits, given this grid geometry and orientation.
- 4.Iterate Steps 1–3 several times, each time rotating the initial grid around its central point by some small angle dμ (up to a total rotation of 60°). Record the number of optimum grid points required for each orientation, and select the orientation that requires the fewest. As a small refinement, if any of the optimum grid points are found to lie with centers outside the distribution of synthetic orbits in the (dαn, dδn) plane, we shift its center to the nearest (dαn, dδn) point from a synthetic orbit. After this refinement, the set of optimum grid points from this orientation is defined as our final optimum set.
These grid points are points in the final (dα, dδ) plane. We can back out the implied set of optimum (
, 
) rates by simply dividing the respective components of each grid point by the total observational baseline tb.
The rotation in Step 4 is included to address boundary conditions. If the region of (dα, dδ) searched by a grid is not ≫max, then there will be a preferred orientation of the grid. However, if max is small relative to the total region, this step is unnecessary.
If, instead of trimming the search grid to as small a number of points as possible, it is preferable to create a "safe" grid that oversearches the boundaries in rate-space by some comfortable margin, then modify Step 2 to select grid points that lie within a larger distance N × max of any final shift vector. Step 3 must be skipped in this case. This maintains the grid spacing and geometry, but adds a "buffer" region around each synthetic orbit.
2.3.3. Long-arc solution
In the case that the observational baseline tb is long enough such that nonlinear components of motion are nonnegligible with respect to , a different treatment is required. Instead of presuming that the largest offset between the motions of a given orbit and the nearest grid point will occur in the final frame n, we must instead ensure that for every synthetic orbit there is one nonlinear shift vector which matches its motions within in every frame. While nonlinear components of motion are nonnegligible, the dominant factor driving the number of shift vectors required is still the linear components of motion.
- 1.After creating an initial sample of synthetic orbits, generate a search grid via the short-arc (linear) solution for the observations. When generating the search grid, set the effective maximum tracking error to be slightly smaller than the desired final maximum tracking error. This is to take into account that the final grid will be selected from a discrete distribution; namely, the motions of the synthetic orbits in the initial sample.
- 2.For each grid point p on the optimum linear search grid, take the list lp of orbits uniquely matched to that grid point and consider their shift vectors in every imaged epoch. For each orbit k in list lp, find all the other orbits within lp that have motions that are matched to the motion of k within the desired maximum tracking errormax in every frame.
- 3.Select the single synthetic orbit k1 with motions matched to the largest number of unique synthetic orbits in lp and set it aside as an optimum orbit.
- 4.If some of the orbits linked to the grid point p are not recovered by the first optimum orbit k1, remove all orbits recovered by the orbit k1 from consideration, then repeat Step 3 to select a second optimum orbit k2.
- 5.Repeat Step 4 until all of the orbits linked to the grid point p are accounted for by one or more optimum orbits.
- 6.Repeat Steps 2–5 for each grid point in the initial solution.
The shift vectors of the final optimum set of orbits are then identified as the optimum set of nonlinear shift vectors to apply to the given data set.
2.3.4. Grid Tree: Optimization for Multinight Arcs
Other significant optimizations can be made to reduce the total number of pixel additions required to search a given data set. One particular optimization that our probabilistic method lends itself to is the use of a multilevel grid tree described by Allen (2002). This approach recognizes that for multinight stacks, there are essentially two distinct timescales: the unit time periods over which images are obtained (length of a single night), and the total observational baseline. Allen (2002) points out that a significant reduction in total number of required shift vectors can be obtained if, instead of combining all images obtained over several nights with one search grid, images are combined on a per-night basis (with a single-night grid), then these stacks are combined with a second search grid to remove the motion that would have occurred over the entire period. If all else remains equal, this tree-of-grids structure results in fewer overall pixel additions, as redundant combinations of images taken on a single night are only performed once. Since this method relies on combining images that are the combination of other images, it can only be applied for combination algorithms that can nest: for example, weighted averaging and co-adding are usable, but not medians.
It is important to note that in combining a tree of grids, the resulting maximum tracking error max goes as the sum of each level's tracking errors 
. As such, each level's tracking error must be limited to i = max/N, where N is the total number of levels. In the case described here, the tree has two levels, and so the effective tracking error in each level must be limited to i = max/2. Thus, while this two-level case requires fewer pixel additions, it results in significantly more pixels to search than in the one-level case.
While we have not implemented this method, we can estimate its results. Since the number of grid points N is approximately proportional to t2/2, and the two-level case requires i = max/2, we can estimate that this method will require 4 times as many shift vectors per unit time interval as a one-level grid. If we consider n nights of data, where the length of a night is equal to the length of a day such that the total observational baseline is tb = (2n - 1) in units of the length of a single night. Thus, if the one-level grid requires N1 shift vectors for the full observational time line, we can estimate that the two-level grid for the full length of time will require N2,full ≃ 4N1. The two-level method will require ∼4 times more pixels be searched in the final stacks than the one-layer method. However, to determine the improvement in pixel additions, we need to determine the size of each nightly grid:
For surveys where the number of observations taken per night Nobs≫N2,full/N2,nightly = (2n - 1)2 such that the number of pixel additions is dominated by the nightly stacks, we can estimate the total number of pixel additions required by:
Whereas for the one-level grid, the number of pixel additions required is N1,add = N1 Npixels Nobs n. The resulting improvement in the number of required pixel additions is roughly
So for a two-night arc, we estimate a factor of ∼2.25 improvement in the number of pixel additions, while for a three-night arc this increases to ∼6.25.
3. COMPARISON OF PROBABILISTIC SOLUTION TO PREVIOUS SURVEY GRIDS
In § 3.1, we will apply our probabilistic shift-vector generation method to previous surveys' observations, and compare the results to the methods used by the authors of the surveys. In order to determine the limits of parameter space searched by each survey, we generate our initial sample of synthetic orbits over very large ranges of heliocentric distance (20–500 AU), inclination (0°–180°), and eccentricity (0–0.999). The sky motion of each synthetic orbit generated is then tested to ensure that it is within the (
, ϕ), or (
, 
) ranges searched by the authors; as such, we only perform this characterization for surveys for which we can accurately reproduce the original search grid from the literature.
Three surveys are selected for comparison: Fraser et al. 2008, Fraser & Kavelaars 2009, and Fuentes et al. 2009. These surveys are selected because they cover relatively large areas (0.25–3 deg2) to relatively faint limits (R ∼ 26–27), and because they contain very complete information regarding their respective search grids and targeted orbital parameter space, which are outlined in Table 1. The comparisons between the original and optimized grids and the heliocentric distance range our analysis indicates each survey was sensitive to are listed in Table 2.
 |
 |
3.1. Fraser et al. 2008
Fraser et al. 2008 searched ∼3 deg2 (taken over baselines ranging from 4 to 8 hr) for TNOs using fixed grids of rates and angles, searching the final stacks by eye. The images were acquired from several facilities, and detailed grid information is only supplied for the MEGAPrime observations. These observations spanned ∼4 hr, and the search grid applied (described in Table 1) required 25 rates and angles be searched. The adopted search-grid spacing was fixed in 
and dϕ, resulting in a maximum tracking error as a function of 
as described in § 2.1. We estimate that the resulting maximum tracking error was ∼1.6'' ≃ 2.2Γ, resulting in a maximum S/N degradation factor of F ≃ 0.57. As the authors reported no sensitivity loss as a function of rate of motion, we adopt this as our uniform maximum tracking error. The authors state that their selection of rates and angles were designed to detect objects on circular orbits with heliocentric distances from ∼25–100 AU with inclinations as high as 70°.
Due to the short arc and large tracking error of this survey, our optimum grid (with 19 shift vectors) is not radically more efficient (∼25%) than the original grid of 25 shift vectors. The original and optimum grids, along with our initial sample, are illustrated in the left two top panels of Figure 2.
Fig. 2.— Results from survey simulations. Top panels: Fraser et al. 2008. Middle panels: Fraser & Kavelaars 2009. Bottom panels: Fuentes et al. 2009. Gray distributions represent density in (dα, dδ) of synthetic orbits generated as described in text. Left panels include overlay of each survey's search grid. Middle panels include overlay of our derived "optimum" search grid. Right panels show inner and outer limits in heliocentric distance d vs. inclination i derived for each survey. Cross-hatched region represents limit variation due to eccentricity. Black triangle represents the distance at discovery and inclination of 2008 KV42.
Based on our simulation of this survey, we find that the minimum heliocentric distance the original grid is sensitive to is a strong function of inclination, with dmin ≃ 22 AU for i = 0° orbits, climbing to dmin ≃ 30 AU for i = 70°. This grid is in fact sensitive to inclinations as high as 180° outside of 36 AU. The outer edge is similarly modified by inclination, but in all cases greater than the 100 AU goal: at the lowest, it is sensitive to objects at distances as high as 164 AU for i ∼ 0°, and at its highest it is sensitive to objects as distant as 184 AU for i ∼ 180°. These limits are illustrated in the top right panels of Figure 2.
3.2. Fraser & Kavelaars 2009
Fraser & Kavelaars 2009 searched ∼0.25 deg2 (taken over a ∼4 hr baseline) for TNOs with a grid of 95 rates and angles, aiming to be sensitive to objects on circular orbits with heliocentric distances from ∼25–200 AU. As in Fraser et al. 2008, the adopted search-grid spacing was fixed in 
and dϕ, and the final stacks were searched by eye. The authors state that they choose the (
, dϕ) spacing to limit the maximum tracking error to ∼2Γ; we verify that at maximum 
after 4 hr, the maximum tracking error of this search grid is approximately 1.25'' ∼ 1.8Γ, resulting in a maximum S/N degradation factor of F ≃ 0.61. As the authors reported no sensitivity loss as a function of rate of motion, we adopt this as our uniform maximum tracking error.
Based on our simulation of this survey, we find that the original search grid is overambitious, as the high-ϕ regions of the grid are not populated by any real objects. No bound orbits observed on the ecliptic at opposition with heliocentric distance outside of d > 28 AU can have angles of motion as high as 15° from the ecliptic (equation [5]). This, coupled with the oversampling at low 
, relative to the maximum tracking error of 1.25'', resulted in excess shift vectors compared to our optimum solution. The optimum grid requires 28 shift vectors, which represents an improvement of over a factor of 3 compared the 95 shift vectors searched by the authors. The original and optimum grids, along with our initial sample, are illustrated in the left two middle panels of Figure 2.
Similar to Fraser et al. 2008, the minimum heliocentric distance the original grid is sensitive to is a function of inclination, with dmin ≃ 24.5 AU for i = 0° orbits, climbing to dmin ≃ 31 AU for i = 70°. This grid was also sensitive to inclinations as high as 180° outside of 36.5 AU. The outer edge is significantly more distant than claimed: at the lowest, it is sensitive to objects at distances as high as 350 AU for i ∼ 0°, and at its highest, it is sensitive to objects as distant as 385 AU for i ∼ 180°. These limits are illustrated in the middle right panels of Figure 2.
3.3. Fuentes et al. 2009
Fuentes et al. 2009 searched 0.25 deg2 (taken over a ∼8 hr baseline) for TNOs with a grid of 732 rates parallel and perpendicular to the ecliptic, aiming to be sensitive to objects with heliocentric distances from ∼20–200 AU. The large number of resulting stacks necessitated the use of an automated detection pipeline. Their grid spacing was 
and motions were limited to ± 15° from the ecliptic. Over an ∼8 hr observational baseline, this grid spacing translates into a maximum tracking error of ≃ 0.6'' ∼ 0.8Γ, resulting in a maximum S/N degradation factor of F ≃ 0.76.
Like Fraser & Kavelaars 2009, this survey also oversearches high-ϕ motions. Because of the small tracking error and nonoptimum grid geometry, any small oversearched regions translate into a significant number of additional search vectors. The optimum grid requires 494 shift vectors, which represents an improvement of ∼33% over the original search grid. The original and optimum grids, along with our initial sample, are illustrated in the left two middle panels of of Figure 2.
Similar to Fraser et al. 2008, the minimum heliocentric distance the original grid is sensitive to is a function of inclination, with dmin ≃ 21.5 AU for i = 0° orbits, climbing to dmin ≃ 26.5 AU for i = 70°. This grid was also sensitive to inclinations as high as 180° outside of 33.5 AU. The outer edge is significantly more distant than claimed: at the lowest, it is sensitive to objects at distances as high as 220 AU for i ∼ 0°, and at its highest it is sensitive to objects as distant as 245 AU for i ∼ 180°. These limits are illustrated in the bottom right panels of Figure 2.
4. EXAMPLE OF APPLICATION: LSST DEEP FIELDS
A component of the LSST survey strategy will be to point to a single field and take hundreds of ∼30 s exposures, then later combine them in software to search for faint moving objects (Chesley et al. 2009). On a single winter night, a single 9.6deg2 field at opposition can be observed for roughly 8 hr, resulting in 850 images with a combined depth of r' ∼ 28. We will determine the feasibility of stacking multiple nights of data to search for specific TNO populations at even fainter magnitudes.
We have run simulations of similar LSST observations to demonstrate the utility of our method for generating shift vectors. One of this method's chief advantages is the ability to strictly limit the orbital parameters of interest, being certain to sensitize the resulting stacks only to those orbits of interest without creating stacks at extraneous rates of motion. As such, the simulation described here is designed to detect a population that has some orbital parameters that are well-defined; namely, objects in the Neptune 2:1 resonance, or "twotinos." The orbital ranges used to simulate this population are listed in Table 3.
 |
We compare observations made at opposition to those made at 45° away from opposition, which have increasing contributions from nonlinear components of motion. If an opposition field can be observed for 8 hr on a single night, fields at this elongation can be observed for ∼7 hr per night. We simulate observing both fields for one, two, and three nights, resulting in the total observational baselines listed in Table 4. The seeing is assumed to be the projected 75th percentile in r', roughly 0.89'' (Tyson et al. 2009). We adopt max = 0.8Γ, similar to that used in previous automated-detection-based pencil-beam surveys, resulting in max ≃ 0.7''.
 |
Table 4 contains the results of our simulations. We have estimated the required number of shift vectors, pixel additions, and pixels to search (without the additional optimization from the tree of grids discussed in § 2.3.4). Figure 3 illustrates the optimum nonlinear grids generated for both elongations. Also illustrated is the magnitude of the error between the final positions of real sources and the linear extrapolation of their position from their initial motion (
, ϕ0).
It is useful to compare the total required number of pixel additions (Nadd) and pixels to search (Nsearch) to the most computationally-intensive digital-tracking survey to date. The HST survey performed by Bernstein et al. (2004) required Nadd ∼ 1016 additions and Nsearch ∼ 7 × 1013 pixels with 2004 computer technology. Because the LSST psf is significantly larger than that of HST, fewer shift vectors are required for the same observational baseline. In our LSST Twotino simulations, even the three-night opposition solution requires as many or fewer Nadd (6.4 × 1016 additions) and Nsearch (2.5 × 1013 pixels). When we include the potential improvement by a factor of ∼6.25 in the total number of pixel additions required with the addition of a second-level grid to the three-night arc (equation (12)), it is clear that even with technology that will be over a decade out of date by the time LSST comes online, multinight stacking of LSST data will be feasible for some TNO populations.
Generalizing to the entire TNO population does not vastly increase the computational demand. We repeated the 8 hr opposition simulation, allowing a to vary uniformly over the same range as d, and the number of required shift vectors increased by approximately 15%. This increase is largely driven by the addition of low-eccentricity objects at low heliocentric distances.
5. CONCLUSIONS
We have shown that by generating search grids for pencil-beam surveys in a probabilistic way and implementing an optimized grid geometry, significant reductions (25%–75%) can be made in the total number of operations required, and using a tree-of-grids structure can reduce the total pixel additions required even further. In addition to improvements in computational requirements, this method also ensures uniform sensitivity to the entire targeted volume of orbital parameter space.
This method also provides a simple way to include nonlinear components of motion necessary to track objects over long observational baselines. This combination of reduced computational requirements and simplicity of extending the observational arc makes the prospect of stacking multiple nights of data appear completely feasible. The additional depth and improved precision of measured orbital properties derived from multinight stacks would be hugely beneficial if used in upcoming large-area surveys like LSST and Pan-STARRS. Simulations of the LSST deep observations show that, for certain TNO populations, it would require less computational effort to obtain the same depth as Bernstein et al. (2004) over 500 times the area by stacking multiple nights of LSST data with the methods described here.
Additionally, this method allows for simple characterization of the orbital sensitivity of existing pencil-beam surveys in literature. From our characterization of the surveys by Fraser et al. (2008), Fuentes et al. (2009), and Fraser & Kavelaars (2009), we have shown that the orbital sensitivity of these surveys varied from what was advertised. Their extra sensitivity to the motions of distant sources makes these surveys sensitive to populations like distant members of the Scattered Disk and Sedna-like objects, and their nondetection provide useful upper limits on these populations. We will explore these upper limits in future work.
Since the discovery of extremely high-inclination, low-pericenter objects like 2008 KV42 (i ≃ 110°), it is important to understand these surveys' sensitivity to such populations. We note that Fuentes et al. (2009) is the only survey characterized here that was sensitive to objects with the same inclination and distance at discovery as 2008 KV42 (see Fig. 2, bottom-right panel), but we contend that it is unlikely that they would have recognized any detection as belonging to this retrograde population. Since follow-up has been performed only rarely for these deep surveys, orbital inclination and heliocentric distance usually remain degenerate, with the prograde solution lying at lower distances (though this degeneracy is rarely acknowledged). It is possible that detections labeled as high-inclination prograde objects may in fact be objects with retrograde orbits like 2008 KV42 at greater distance. Only multinight arcs or follow-up at later epochs can break this degeneracy and clearly identify objects belonging to rare TNO subclasses.
Alex Parker is funded by the NSF-GRFP award DGE-0836694.
Footnotes
- 1
This assumes that the object's distance and eccentricity are not such that
–in other words, that the object is not overtaking Earth at opposition. If this were the case, is lowest at i = 180°.
