The Kuiper Belt Luminosity Function from m(R)=21 to 26

We have performed an ecliptic imaging survey of the Kuiper belt with our deepest and widest field achieving a limiting flux of m(g') = 26.4, with a sky coverage of 3.0 square-degrees. This is the largest coverage of any other Kuiper belt survey to this depth. We detect 72 objects, two of which have been previously observed. We have improved the Bayesian maximum likelihood fitting technique presented in Gladman et al. (1998) to account for calibration and sky density variations and have used this to determine the luminosity function of the Kuiper belt. Combining our detections with previous surveys, we find the luminosity function is well represented by a single power-law with slope alpha = 0.65 +/- 0.05 and an on ecliptic sky density of 1 object per square-degree brighter than m(R)=23.42 +/- 0.13. Assuming constant albedos, this slope suggests a differential size-distribution slope of 4.25 +/- 0.25, which is steeper than the Dohnanyi slope of 3.5 expected if the belt is in a state of collisional equilibrium. We find no evidence for a roll-over or knee in the luminosity function and reject such models brightward of m(R) ~ 24.6.


Introduction
The study of extrasolar debris and dust disks has revealed that, for at least some of these disks to exist as we see them, there must be a source which is responsible for replenishment of small-grain dust in the disk. Otherwise, due to radiation pressure, the small grain dust would be blown out of their stellar systems on time scales shorter than the age of the star. A disk of planetesimals embedded in the dust disk which is undergoing collisional evolution is a likely source of this dust. Disruptive collisions could produce the necessary influx of dust to extend the lifetime of the disk beyond the dust blow-out time. See Meyer et al. (2006) for a current review of debris disks.
The Kuiper belt is analogous to these extra-solar planetesimal disks, and provides an excellent laboratory to study and understand the properties of these planetesimals and the processes that affect them, including collisional processes, tensile strengths, compositions, and the mechanisms by which they formed. Knowledge of the size distribution in the belt can constrain much of this information. The size distribution of small objects provides information on the bulk material properties of the objects, (O'Brien and Greenberg, 2003;Kenyon and Bromley, 2004) while the size distribution of large objects can provide information on the conditions under which these bodies formed (Kenyon, 2002).
We performed a 3.0 square-degree survey of the Kuiper belt to determine the size-distribution of large (D 50 km) Kuiper belt objects (KBOs) via a measure of the belt's luminosity function. This is the deepest photometric Kuiper belt survey of this size ever completed.
In section 2 we describe our observations. In section 3 we describe our search technique and data analysis. In section 4, we derive a relation between the size-distribution and the luminosity function. In section 5, we describe the statistical analysis used. We present our results in section 6. In section 7 we discuss the implications of our findings, and in section 8 we present our conclusions.

Observations
The observations used in our survey were taken with the 3.6 m Canada-France-Hawaii Telescope (CFHT) and the Cerro-Tololo-Interamerican Observatory (CTIO) 4m Blanco telescope. Observations at CFHT were acquired with both the CFH12k (0.33 square degree fov.) and MEGAPrime (0.88 square degree fov) mosaic cameras while observations at CTIO were made with the Mosaic2 camera (0.38 square degree fov), providing 0.67, 0.84, and 1.54 square degrees of searchable area respectively for a total of 3.0 square degrees of searchable area. Details of the observations are provided in Tables 3 and 4. All of these observations were originally acquired for use in searches for satellites of Uranus and Neptune Holman et al., 2004) and prior to this work, none of these fields had been searched for KBOs. The observations were made when Uranus and Neptune (and all KBOs in each field) were at or near opposition, and covered the projected Hill-spheres of each planet. The CFH12k and Mosaic2 observations excluded the area closer than 3 ′ to the planets, to avoid scattered light in the images. These observations all occurred near the ecliptic and are well suited to a deep search for Kuiper belt objects (KBOs).
Approximately 10-20 bright, non-saturated stars (∼ 20 mag.) per chip were used as reference stars for the image reductions. The variation in the average reference star magnitudes follow approximately that expected from the varying airmass of the observations (see Fig. 1) indicating photometric conditions during the observations. The remaining scatter is ∼ 0.02 − 0.04 mag, significantly less than the shot noise for the brightest objects detected.

The Data Processing and Image Search
To determine the behavior of KBO size distribution requires the discovery of a large number of faint moving sources. If long exposures are used, then the sources will move far enough (more than the size of the seeing disk) that a trailed image results and no additional depth will be achieved. In our deep searches we have adopted a strategy of taking exposures short enough that trailing losses will be negligible. We then shift the individual exposures to account for the expected sky motions of the objects of interest. These shifted images are then combined together to achieve the depth needed. To account for the various possible sky motions of KBOs, we have shifted the images at a variety of rates and angles and then visually examined each of the combined images (stacks), searching for point like sources (see Fig. 2).

Data Preprocessing
MEGAPrime: CFHT provides wide field imaging data from the MEGAPrime camera in a 'preprocessed' format, ready for science exploitation. The frames provided have been processed using the CFHT ELIXIR/FLIPS (Magnier and Cuillandre, 2004) processing system. As part of this processing, unique bias, flat-field, fringe, and scattered light images are produced for each 'camera-run' (typically matching a 14 night dark-run) and applied to all frames acquired during that camera run. Dark runs are broken into multiple camera-runs if a significant change in the camera performance or image characteristics is detected. All the MEGAPrime images used in this project were acquired within a single camera-run and have been 'detrended' with a common set of calibrator images. Standard calibration from CFHT results in a flux conserved image with constant sky level across the image that has a typical variation of ±3%. CFH12K: These data were originally acquired as part of a search for irregular satellites of Uranus  and are different from the data used for the wide field satellite search presented in Petit et al. (2006). Due to the strong time constraints of the imaging, the images were acquired in 'classical' observing mode and not automatically detrended by the CFHT queued service observing team. In November 2004 the Canadian Astronomical Data Center (CADC) acquired the ELIXIR/FLIPS software and calibrators for the CFH12K data and subsequently 'detrended' the entire set of CFH12K images for which global calibrator frames (bias, flats, fringe, scattered light) were available. Part of the re-processing effort included the re-processing of the CFH12K frames used in this project. These reprocessed images were used in this search and provide a sky flatness of typically ±4%. CTIO: The CTIO images were originally acquired as part of a project to search for satellites of Neptune . Bias frames were acquired on each observing night and a combination of the overscan strip and an average of a dozen bias frames was used to remove the instrumental ADC bias from each frame. During the period of observations, a number (∼ 15) of independent fields were observed in order to measure the astrometric positions of some previously-known Kuiper belt objects. These sequences of fields were 'median combined' using the IRAF images.median task with high pixel clipping. These combined images provided an excellent flat-field frame that was divided into the search images using the IRAF (Tody, 1993) mscred.ccdproc task. This process resulted in images which have flat sky across the entire mosaic with a typical scatter of ±3% in the sky flux level.
The curvature of the focal plane for all instruments used in this project was small enough compared to the spatial shifts applied to the images during image processing, such that spatially flattening the images was unnecessary.
For the MEGAPrime and CFH12k data, zeropoint calibrations were done with the standard Elixir routines and were provided by the CFHT Elixir QSO (Magnier and Cuillandre, 2004) and were reported to be Z f ield,g ′ = 26.46±0.02 mags. and Z f ield,R = 26.22 ± 0.02 mags. Hence these magnitudes are presented in the Sloan g' (Smith et al., 2002) and Kron-Cousins R. For the N10032W3, N11033, NEP0813NW3, and NEP0815NE3 fields, calibration frames were unavailable, and a nominal zeropoint for the CTIO mosaic of Z V R = 26.0 was used presenting the magnitudes of the objects discovered in these frames in the VR filter.

Artificial Object Planting
To determine the search efficiency as a function of magnitude, artificial objects representative of KBOs were added (implanted) in the images.
The sky rate of motion (θ) of an object on a circular orbit observed at opposition on the ecliptic, at heliocentric distance r can be approximately given byθ (1) We implanted 100-150 moving sources into each CCD image. All sources were added blindly to the data before searching began; the artificial source lists were revealed to the operator only after searching was completed. The rates and angles of motion of the artificial sources (approximately 1.3 − 4.1 arcsec. hr −1 and ±10 o from the ecliptic) were typical of objects on circular orbits at heliocentric distances from 25 − 100 AU with inclinations between 0 − 70 o . Each implanted source was given a randomly selected flux, equivalent to that of a source between 23 −27 mag. The distribution of artificial sources is shown in Fig. 3. Additionally, five artificial sources with flux levels equivalent to 21st magnitude sources were implanted, with zero inclination. These objects have sufficient flux to allow us to flag errors in the image combining algorithms (failed image subtraction, wrong mask limits, etc.) in advance of searching the data, and provided reference moving sources for computing aperture corrections in the final image combinations.
To account for image-to-image flux variations due to changes in airmass and possible sky transparency, the flux of the planted artificial objects was varied with respect to a reference image (usually the middle image in the list) to match the average flux variation of the reference stars (see Fig. 1).
The point spread function (PSF) for each frame was approximated on a perchip basis using the stellar profile of a single bright isolated star. The PSF variation across individual images was small enough that the use of a PSF that varied with position was unnecessary.
Artificial sources whose sky motion would result in a drift of the centroid of the PSF by more than one pixel during an exposure were represented by a series of fainter sources. The number of faint objects was equal to the number of pixels the source would drift during the exposure. The total flux of all the fainter sources was equal to the flux of the original source. In this way, we fully included trailing effects into the implanted sources, and thus accounted for the effects of trailing in our final search efficiencies.

Image Subtraction
During our previous deep searches (Gladman et al., 1998(Gladman et al., , 2001Holman et al., 2004;Kavelaars et al., 2004) we found that the two largest inhibitors to this search method are trails in the combined images (caused by bright stars) and the enormous human effort required to visually examine the broad range of rate/angle combinations needed to ensure detection of KBOs. To combat these issues, we utilize an image subtraction routine to remove most stationary background sources and implement a new image display method which greatly eases the strain on the user during the visual image search. We describe the image subtraction routine, and the image display method here.
A template image was subtracted from each image to remove stationary sources from individual images prior to being stacked, thus improving our ability to detect moving sources. To create the image subtraction template, an artificial skepticism weighted average of all images (Stetson, 1989), with the artificial moving sources already added, was created. This method creates a per-pixel weighted average, with weights calculated iteratively, and quickly converges on an average which places very little weight on spurious pixels such as cosmic ray hits.
For our image subtractions, we chose to use psfmatch3 developed by one of us (CJP) for the Supernova Legacy Survey (Pritchet, 2005;Astier et al., 2006). This subtraction routine compensates for variations between the image quality of the template image and that of the individual images. Several programs that perform this task already exist (see Alard (2000) and references therein). Psfmatch3 has several advantages. The subtraction kernel can have arbitrary form, and does not require representation by a set of basis functions to perform the subtraction. In addition, it can automatically remove both spatial and background variations between images. The result of the subtraction routine, is smooth subtracted images with zero average backgrounds. A detailed description of the method is presented in Appendix A.
While the image subtraction removed any stationary non-saturated sources, cosmic ray spikes and other spurious hot pixels remained. To compensate for this, an image mask was applied, such that if a pixel had a value outside a chosen range, then a 3 pixel by 3 pixel box about that pixel was set to have zero value. The lower limit to the range was chosen to be -5 times the background noise (the average background was set to zero in the image subtraction). The upper limit was chosen to be ∼ 25000 counts, such that most saturated regions and cosmic ray spikes were masked out of the data. This procedure masked the centres of the brightest KBO sources. This however, did not hinder the detections of these sources as they were still glaringly obvious even in individual frames.

Image Stacking
The subtracted, masked images containing artificial sources were shifted before they were stacked, such that each subsequent image was spatially shifted to compensate for the predicted motion of KBOs in the frame. Sources whose sky motion was well matched by the shifts applied to the images appear nearly round in the stack (see Fig. 2) while any residual flux from stationary sources was trailed.
If a source's rate of motion was not well matched by the spatial shifts applied to the subtracted images, that source appeared extended in the stack. This characteristic trailing provided a very robust means of discriminating between real sources and noise. Noise sources were produced during the shifting and stacking when positive flux regions from different images were caused to overlap by the choice of spatial shifts. These false sources, however, did not show the image shape variation characteristic of a real source; false sources did exhibit trailing, but of a different length and width compared to a real (or artificial) source. As such, false sources were not selected as candidates by the trained operator.
The template image used for the Psfmatch3 subtraction process contains both the real and artificial moving sources. These appear as faint trails in the template image. Subtracting this template from the input images results in a low flux area behind each moving source. This feature only occurs around the positions of real (and artificial) moving sources and provides an additional and robust means of source-noise discrimination (see Fig. 2) and was required to be present in order to mark a source as a candidate.
The quality of the image subtraction, and hence the final searchable images suffered from a few effects: bad columns in the CCD, bleeding from saturated stars, and regions of poor subtraction around bright galaxies. These problems were exacerbated by image shifting; the unsearchable area of these regions were expanded, reducing the overall searchable area. Gradients of the background caused by bad image subtraction around bright sources, produced images that were difficult to display, reducing the search efficiency. Our efficiency of detecting artificial, planted sources accounts for all reductions in area as the artificial sources were planted at random locations occupying the full spatial extent of each CCD image.

Visual Search
To search the combined stack of shifted images, the stacks were divided up into ∼ 200 image subsections with a 20 pixel overlap between neighboring subsections. One-by-one, a grid of rates and angles like that shown in Fig.  2 was displayed for each image subsection. Each grid was searched by eye, and sources were recognized by their characteristic trailing and subtraction wells. We found that a five-by-five grid of shift rates and angles maximized the detection efficiency while minimizing the time spent searching. The low variations in the sky background, resulting from the template subtraction, and an image display tool developed specifically for this searching allowed us to rapidly search many data sets.
The pixel coordinates of potential candidates were recorded along with the rate and angle combination that produced the most circular image. A candidate was selected as a planted object if its marked image location was within 10 pixels of an artificial source location. The list of detected artificial candidates was used to determine the detection efficiency as a function of artificial source brightness. The detection efficiencies for each chip of a given field were averaged together to provide a global, per field, detection efficiency which we then modeled (see Fig. 4 and Table 5) using the equation Our deepest and widest field has a η(m) = 50% threshold at m(g ′ ) = 26.4 and a sky coverage of 0.84 square degrees. The sky coverage of our combined data is 3.0 square-degrees.
Remaining candidates not marked as artificial source detections were re-examined at a finer grid of 8 rates and 8 angles to further discriminate between real sources and noise. Candidates were rejected as false positives during this process if the variation of their image shape and trailing length with rate and angle did not match that of brighter sources with similar apparent rates of motion, whose detection was more robust. While this approach allowed us to isolate and remove false detections, the artificial sources were treated differently from the real sources, as we did not examine candidates marked as artificial sources in the finer rate and angle grid. The brightnesses of all false positives were found to be faintward of the 50% threshold of their discovery field, when their fluxes were measured in the same way as all real detections.
Therefore, below the 50% threshold, our search efficiency was no longer representative of the true search efficiency, and subsequently, sources fainter than the 50% threshold of each field are ignored in our analysis. [ NOTE: Here we define the 50% threshold as the point where η (m) = 50%.] A series of follow-up images of the MEGAPrime field as well as the NEP0813NW3 CTIO field were obtained one or two nights after the initial discovery images.
Using the short first night discovery arcs, we projected the motion of each detection forward to the time of the follow-up observations in order to predict the sky location of the source on the second night's images. For 3 of the 25 sources in the MEGAPrime field, the predicted location placed the sources in the gaps between the CCDs of the MEGAPrime MOSAIC. Aside from these 3 sources, all of our initial detections brighter than the 50% threshold were confirmed. Only the faintest initial detection in the MEGAPrime field was not confirmed on the second night, likely due to poorer average seeing conditions during the second nights observations. All 6 detections in the NEP0813NW3 field were successfully confirmed on the second night.
None of our other search fields had follow-up observations. We were confident, however, that all detections in the other fields brightward of the 50% threshold were real because all detections brightward of the 50% threshold (excluding those that fell on chip gaps) in the MEGAPrime and NEP0813NW3 fields were confirmed, and all fields were processed and searched with identical procedures.
For the range of rates and angles that we planted objects into the data, we found no significant variation in detection efficiency versus rate or angle of motion for any of our fields (see Fig. 3). This is important to note because any significant variation in efficiency with rate of motion needs to be included when deriving the size distribution from the observed luminosity function. We did not search for sources with rates of motion consistent with objects at distances further than 100 AU as the search method we employed is not sensitive to sources beyond this distance; such distant objects would not show enough apparent motion over the time-span of our observations to detect the faintest sources, and a degradation of the detection efficiency with distance would occur. Therefore, our determination of the luminosity function only applies to KBOs closer than 100 AU.

Source flux measurements and detection confidence
To measure the flux from each detected source, we stacked the non-subtracted images containing the artificial objects. These images were shifted at the rate and angle that produced the roundest image for the source whose brightness was to be measured. Three sets of stacks were produced, by averaging the images, while using a pixel by pixel cut that threw away pixel values outside 3-sigma of the mean value. Three averages were made from images from the first half, middle half, and last half of the observing period for each field. These average images provided three separate, mostly uncorrelated, measurements of each source's flux. Aperture photometry was performed using the IRAF daophot.phot task, set to use an aperture of 4 pixels radius (close to the size of the FWHM of point-sources in our images). Using the five 21st magnitude planted 'reference' objects, and the mkapfile IRAF routine, aperture corrections were determined for each of the three image combinations (first, middle last) and were used to correct individual flux measurements for each source. The IRAF mkapfile task reported our aperture corrections to an accuracy of σ aper 0.03. The magnitudes, as measured on each exposure set, were averaged for the final reported magnitudes (see Table 6). We did not measure the flux in a particular stack for sources that were within a few seeing disks of bright or bloomed stars or galaxies (∼ 1 in 10 possible source measurements), in that particular stack.
For two objects in the MEGAPrime data, the measured magnitudes varied significantly more than the uncertainty of each individual measurement, indicating a possible light curve. For these two objects, the magnitudes measured from the first image stack were used in our determination of the luminosity function, as was done in Bernstein et al. (2004). The second night's data provided an additional 3 magnitude measurements of any followed-up source. For all but the variable sources, we used the average magnitude measured from all image stacks including those from the second night.

Characterization
As required by the maximum likelihood routine, we parametrize the magnitude uncertainty and detection efficiency of our survey. The scatter between each artificial source's measured and inserted magnitudes was used to determine the precision of our flux measurements.
For background limited sources, the uncertainty of a source's measured magnitude can be represented by where Z is the telescope zeropoint. γ is a function of the background b, the camera gain g, and the on-source integration time, and is fit to the scatter of the artificial source's measured and inserted magnitudes. The fit values of γ for each field can be seen in Table 5. The magnitude measurement scatter for each source was then drawn from Eq. 3 using these best fit values.
The theoretical value of γ is γ * = 2.5 ln 10 where N is the number of exposures. If we use the typical values from our MEGAPrime observations, we expect a value of γ * = 0.12. The measured value of γ = 0.27 is clearly larger. Newberry (1991) however, has shown that the noise calculation from Eq. 4 is incomplete, and does not fully account for the noise introduced into observations during the data reductions. They found that the noise estimate given in Eq. 4 can be too small by a factor of 2. Similarly, source fluxes were measured off of images that were shifted and stacked. If the shifts were slightly different than the source's true motion, the measured magnitude would differ slightly from the true magnitude. Thus we find that the uncertainties measured from the artificial sources are close to the expected values. The net uncertainties used in our luminosity function determinations are the 1-sigma shot-noise, ∆m, aperture correction (σ aper 0.03), and calibration uncertainties, added in quadrature.

The Luminosity Function
The size distribution of the KBOs can be determined directly from their luminosity function. The relation between the luminosity function and the size distribution has been derived previously (see, for example, Gladman et al. (2001)). We reproduce the relation and discuss complications due to possible variations of KBO albedos, using a slightly different approach than that in Gladman et al. (2001).
Consider the number of objects N in a survey field, given by where R (r) and S (D) are the radial and size distribution functions, r and D are object heliocentric distance and diameter, and A is some normalization constant. We model the size distribution as S (D) ∝ D −q . The magnitude m of an object at heliocentric distance r, geocentric distance ∆, with diameter D and albedo p λ is where K λ and p λ are wavelength dependent; in R-band, p R ∼ 4% gives K R ∼ 18.8. The smallest observable object at a given distance r in a survey with limiting magnitude m max , has diameter D min = p −1/2 λ r 2 10 (K λ −mmax)/5 where we have made the approximation r ≡ ∆. Similarly, the largest observed object at distance r will have diameter D max = p −1/2 λ r 2 10 (K λ −m min )/5 . Using these limits and integrating, Eq. 5 becomes if q = 1. Note: this assumes that D min /D max are not the smallest/largest objects in the belt; if this were the case, using D min /D max as the integration limits as we have defined them here is incorrect.
The observed cumulative luminosity function is well represented by a powerlaw of the form which gives the number of KBOs per square-degree, on the belt mid-plane, brighter than magnitude m. Comparing this with Eq. 8, reveals that the luminosity function slope α and the size distribution slope q are related by From Eqs. 8 and 10 we see that, for a size distribution that is independent of distance, the choice in radial distribution is not significant and only affects the interpretation of the normalization of the observed luminosity function, not the inferred size distribution slope.
Here we have assumed that the KBO albedos do not vary with distance, or size. If the distribution of KBO albedos varies only with distance, then this will only affect the normalization of the luminosity function.
The albedo of Pluto and Eris (D 2000km) are p P luto ∼ 0.6 and p Eris ∼ 0.6 − 0.86 (Young et al., 2001;Bertoldi et al., 2006;Brown et al., 2006) while smaller objects (D ∼ 100 km) have been shown to have albedos p ∼ 0.06 (Grundy et al., 2005). These data suggest an increase in albedo with size. We can understand the effects of such a trend on the interpretation of the LF by considering a toy model where KBO albedos vary as p ∼ D −β . This functional form retains the analyticity of Eq. 5, and reveals the effects of albedo variations on the inferred size distribution. Under this assumption, Eq. 8 becomes Thus the slope of the size distribution is We see that in the case of constant albedo (β = 0) we get back Eq. 10.
Current albedo data imply β ∼ −1 (Stansberry et al., 2007). Thus, an estimate of q which assumes β = 0 potentially under-estimates the steepness of the intrinsic size distribution. Our knowledge of the relation between albedo and size however, is currently insufficient to constrain β.
The reader is cautioned that the current determinations of the size distribution from the Kuiper belt luminosity function are based on a few poorly constrained estimates of KBO physical properties and the inferred slope q is probably underestimated.

Maximum Likelihood fits
To determine the optimal values of α and m o , we use a maximum likelihood method. We extend the maximum likelihood analysis of Gladman et al. (2001) to account for observations made in different wave-bands, as well as systematic differences in the normalization, m o , that may affect the results of combining separate data sets, such as systematic calibration errors, sky density variations, and variations in average colour of the observed KBOs compared to the average colour of all KBOs.
Given that the probability of detecting an individual KBO is independent of detecting the next, the likelihood function for a single survey k takes the form where m i is the magnitude of detection i,Ñ k is the number of objects expected to be detected in a given survey, and P i is the probability of having object i given the underlying luminosity function.Ñ k is given bỹ where Ω is the survey area, η(m) = η(m|η max , m * , g) is the detection efficiency for an object with magnitude m, and Σ(m|α, m o ) is the differential density of objects on the sky. P i and Σ(m|α, m o ) are given by and where α is the slope of the luminosity function (LF), m o is the magnitude at which the sky density is 1 KBO per square degree, and ǫ i (m) is a functional representation of the photometric uncertainty for object i. A formal derivation of this likelihood function is presented in Loredo (2004).
We choose to represent the magnitude uncertainties as gaussian. ie.
where ∆m i is the uncertainty in the magnitude measurement of each object. This treatment of uncertainties for faint sources is incorrect, but is a sufficient approximation for detections brightward of the 50% efficiency threshold. Brightward of this threshold, the gaussian approximation will not affect the results of the maximum likelihood inference (Bernstein et al., 2004).
For a group of surveys with calibration uncertainties and colour offset variations much smaller than the uncertainty in the flux measurements of the brightest objects detected, the net likelihood resulting from combining multiple surveys together is the product of the likelihoods of each individual survey. In reality, photometric calibration and colour offset uncertainties are not insignificant, and therefore, need to be considered when combining different surveys.
Additionally, the apparent sky density of KBOs and, hence, m o are strong functions of latitude and longitude (Kavelaars et al., In Press). When combining separate surveys, differences in the flux-limited sky densities will skew the inferred slope of the LF.
We account for these effects with the addition of two new parameters for each survey, the colour parameter C k , and the density parameter ∆m o,k . Eqs. 14, 15, and 16 becomẽ and Σ(m − C k |α, m o , ∆m o,k ) = ln(10) α 10 α(m−C k −(mo−∆m o,k )) .
As can be seen, for the case of a power-law LF, the parameters C k and ∆m o,k are degenerate. For a non power-law model however, C k and ∆m o,k are no longer degenerate, and need to be treated separately. As we consider non power-law LFs in this manuscript, we choose to treat C k and ∆m o,k as different parameters when evaluating the likelihood.
This treatment substantially increases the number of parameters in the fit, while not providing any new information about the shape of the LF. Thus we treat these additional parameters as nuisance parameters, and marginalize

Results
Each detection in our survey is listed in Table 6 along with estimated barycentric distance and inclination determined from fit radec (Bernstein and Khushalani, 2000). The fit radec routine determines a set of orbital parameters that best fit the observations in a least-squares sense. Because of the large degeneracies between orbital parameters for short-arc observations such as those presented here, fit radec initially assumes that the objects are near perihelion, and are on bound, nearly circular orbits, near the ecliptic plane. The routine determines 2-sigma uncertainties, which contain 95% of the orbits which are consistent with the measured positions of the objects. In this survey, we have detected 72 KBOs, 53 of which are brightward of the 50% threshold of our search, and these 53 detections are used in our maximum likelihood analysis.
As the UNE and UNW fields are close on the sky, observed together through the same filter/telescope, and the skies were photometric during those observations, we combine the two fields for the maximum likelihood fits. Similarly, we treat the N10033 and N10032W3 and the NEP0815NW3 and NEP0815NE3 fields in the same fashion. We refer to these combined fields as the UN, CTIO01, and CTIO02 fields. The CTIO01 and CTIO02 fields were not combined together, as they were observed in separate years.
To extend our results, we include the observations from other surveys that have well-measured efficiency functions for each field in the survey. We define the F07 sample as those objects detected in the surveys from Jewitt et al. (1998), Gladman et al. (1998), Chiang and Brown (1999), Gladman et al. (2001), Allen et al. (2002), and Petit et al. (2006), as well as those detected in the UN, CTIO01, and CTIO02 fields.
In the F07 sample we also include all on ecliptic detections in  (T01) in fields with detection efficiencies classified as 'good' or 'medium'. Because the T01 images cover a large range of latitudes, longitudes, and detection efficiencies, we divided the survey into smaller "sub-fields" to avoid large density variations from field to field. Each sub-field spans ∼1 -1.5 hrs. RA. We do not include any of the high-latitude data from T01 as the KBO latitude distribution is not sufficiently understood to be able to predict the density of KBOs at high latitudes and we are thus unable to provide a reasonable estimate of ∆m o,k for the high-latitude fields. We avoid complications due to detection efficiency variations by further sub-dividing the data from T01 into separate groups based on the good and medium detection efficiencies defined in that work. We do not consider any T01 fields with detection efficiencies classified as 'poor'. The sub-field divisions listed by field ID (see Table 2 in ) are presented in Appendix B.
We include all fields from the above surveys, even those with few or no detections in the F07 sample. This is the sample on which all our LF analysis was performed.
We chose to exclude the data from Bernstein et al. (2004) as the intent of this work is to determine the shape of the bright end of the luminosity function.
When computing our LF estimate, we ignored all detections fainter than the η(m) = 50% detection threshold of each individual survey, and truncated ("cut") the detection efficiency to zero faintward of that threshold. Correctly determining the efficiency function at low levels is notoriously difficult and prone to error. Effects such as Malmquist-bias further complicate detection efficiency measurements, at low brightnesses.
To test for a bias in the fits caused by an efficiency truncation, we ran a series of Monte-Carlo realizations that simulated the observation of a luminosity function from a single survey. We also performed this analysis using a set of surveys with parameters like those in the F07 sample. We performed 2000 random realizations for efficiency cuts ranging from 0-90%. We found that, for a single survey, the mean of the estimated LF slope is biased to steeper slopes when truncating the detection efficiencies at some non-zero threshold. The bias was found to decrease asymptotically to zero as the efficiency threshold was moved faintward.
We found the bias caused by including an efficiency cut is removed however, when combining multiple surveys if the separate surveys have different detection efficiencies. Most detections in a survey occur where η(m) ∼ 80%. Thus, when two or more surveys are well separated in their limiting magnitudes, so are their detections, which creates a "lever arm" that dominates the slope determination, effectively removing any efficiency cut bias in the fit.
For each field with 5 or more detections which is the minimum number of detections required to provide a reliable measure of the LF, we determined the best-fit parameters for the LF using Eq. 13 (see Table 7). The 1, 2, and 3-sigma credible regions of these fits are presented in Fig. 5. The 1-sigma uncertainties of the best-fit parameters, taken as the extrema of the 1-sigma credible regions, are presented in Table 7. Presenting the fit parameter uncertainties in this way describes the full range of allowed values for each parameter individually. These uncertainties are somewhat misleading as they do not describe the correlation between α and m o (see Fig. 5).
The average (V-R) colour of KBOs is < V −R >= 0.56 (Hainaut and Delsanti, 2002). Using the VR and R magnitude transformation given by Allen et al. (2001), as well as the relations between V, R, r' and g' given by Smith et al. (2002), we determined the expected offsets between the various filters used in the surveys considered here. We found that the average colours of KBOs are < V R − R >= 0.03, < r ′ − R >= 0.26, and < g ′ − R >= 0.95. Using these colours to offset the separate fields to R-band, we find the estimates of m o provided independently by each survey are not consistent with a single value (see Fig. 5). Small shifts in m o are necessary to make all fields agree at the 1 -2 sigma level. This is acceptable when considering possible sky density variations, and photometric calibration errors (see Section 5) and justifies the more complicated likelihood equation given by Eq. 21. We therefore shifted all data to R-band using the expected average KBO colours, and used Eq. 21 when determining the best-fit LF parameters.
To determine the possible range in m o values between surveys, we consider a toy model of the Kuiper belt. As the majority of the objects observed in the surveys considered in this work are Plutinos and classical belt objects (CKBOs), we only consider these objects in the model. We represent the LFs of both populations as power-laws with the same slope, but with different observed sky densities. Then the observed cumulative LF (assuming constant detection efficiencies) is a power-law, and is given by where m o,c , m o,p , and m o are the magnitudes at which one object per square degree is observed for the CKBOs, the Plutinos, and all populations respectively and f is the ratio of the number Plutinos to CKBOs.
The mean eccentricity of Plutinos isē = 0.15, and they are ∼ 30% as numerous as the classical belt objects of the intrinsic population. ie. f = 0.3 (Kavelaars et al. (In press)). The magnitude of an object as a function of heliocentric distance r is given by m ∼ K + 10.0 log(r). Thus the variation in m o,p for Plutinos at perihelion versus aphelion is ∼ 1 mag. If all CKBOs are observed at r = 43 AU, then the difference in m o between fields that observe Plutinos at perihelion versus aphelion is ∼ 0.8 mags. This is an upper estimate of the true offset because not all Plutinos come to perihelion at once, but provides a reasonable range of integration for each m o,k parameter.
Latitude differences between surveys will cause variations in m o requiring knowledge of the latitude distribution of the Kuiper belt. While the location of the mid-plane is unclear, it is secure that the mid-plane is inconsistent with the ecliptic (Brown and Pan, 2004;Elliot et al., 2005); there is disagreement about whether the mid-plane is the invariable plane. Brown (2001) has shown that the KBO latitude distribution above the plane has a Full-Width-Half-Maximum of ∼ 7 o . Most of the observations considered in the F07 sample are made near the ecliptic. If the true plane of the belt is similar to the invariable plane, or the Kuiper plane proposed by Brown and Pan (2004), then the density variation due to latitude between ecliptic fields is at most 5% and would affect the m o,k values by ∼ 0.05 mags. This variation is small compared to that expected from changes in longitude, and is ignored in our analysis. We bound the integration of Eq. 21 over the m o,k parameters between ±0.4 mags.
The standard deviation of the < V − R > colour from the MBOSS data set is σ(V − R) ∼ 0.15 (Hainaut and Delsanti, 2002). Hence, the expected variation in average object colour between separate surveys due solely to KBO colour variation is ∼ 0.15/ √ N where N is the number of objects in a survey. For the fields considered in the F07 sample, N ∼ 1 − 15. Zeropoint calibration errors between separate surveys have been as high as 0.1 magnitudes. Thus a reasonable range of offsets due to calibrations and colour variation is ∼ 0.2 mags. This offset can occur in either a positive or negative sense on the measured magnitude. Thus we bound the integration of Eq. 21 over the C k parameters between ±0.2 mags.
Maximizing Eq. 21 using the F07 sample, we find a best fit slope of α = 0.65 ± 0.05 and normalization m o = 23.42 ± 0.13. The likelihood contours of the fit are shown in Fig. 6.
The factor e −Ñ in Eq. 21 weights the fit towards the lowest possible number of detections. The m o that maximizes the likelihood equation given by Eq. 21 is the maximal value within the range allowed by the ∆m o,k that best describes each individual field considered. Thus, the best-fit m o is not applicable to any one field, but rather is a value typical of all data considered in the fit.
To test if the fit is consistent with the observations, we employed a series of Monte-Carlo simulations. The simulations involved random realizations of a number of objects drawn from the best-fit power-law model (α = 0.65, m o = 23.42) equal to the number of objects in the F07 sample with the random magnitudes scattered according to our uncertainty model (see Eqs. 3 and 17). A best-fit power-law of each realization was determined using our maximum likelihood technique, and the distribution of maximum likelihoods was determined from 1000 realizations. We found that the probability of getting a maximum likelihood less than or equal to the maximum likelihood computed from the F07 sample was 43.2%: the model is fully consistent with these data.
In Fig. 7 (top) we present a histogram of the differential LF of the F07 sample corrected for detection efficiencies. These data are well described by our best fit. The reader is cautioned however, from drawing conclusions about the LF from this diagram alone. All source magnitudes have been adjusted to R-band using typical KBO colours (see above). The observational data from different surveys however, contain calibration and colour offsets, and the observations have not been adjusted to reflect variations in sky densities, as a standard model of the density variations is not known. The minor discrepancies between the observed LF and the model apparent in Fig. 7 maybe caused by these effects.
In Fig. 7 (middle) we show the net effective area for all fields as well as a differential histogram of the object magnitudes shifted to R-band with completeness corrections. The different fields which are included in the F07 sample all have different 50% thresholds, resulting in the broad fall off in effective area between m(R) ∼ 21 to 26.
In Fig. 7 (bottom) we show the logarithm of the ratio of the observed and fit differential LFs, as a function of magnitude. This figure shows that the LF from m(R) = 21 to 26 is well described by a single power-law. Bernstein et al. (2004) conclude that the slope of the luminosity function rolls over to shallower values for fainter magnitudes and that the luminosity function is well described by a rolling power-law given by the functional form

Broken Power-law Models
where Σ(m) is the differential surface density of objects, Σ 23 is the surface density of objects at m(R) = 23, α is the bright end slope, and α ′ is the slope derivative. They found that a rolling power-law was a better fit to their observations than the single power-law used here. This suggests that a deviation in the form of a flattening of the power-law might be visible in our data.
To look for evidence of a roll-over in the KBO LF we fit Eq. 23 to the F07 sample. As the density parameter for Eq. 23, Σ 23 , is different from m o in Eq. 9, we introduce a density offset parameter ∆Σ k as a multiplicative factor in front of Eq. 23. We maintain the same range in colour and density offsets as used in the power-law fit (see above) and marginalize C k and ∆Σ k over the range ±0.2 and ±0.25 respectively. Our maximum likelihood method gives a best-fit A = 0.79 ± 0.14, α = 0.74 ± 0.09, and α ′ = −0.03 ± 0.04. Bernstein determined a best-fit of A = 1.07, α = 0.66 ±0.03 α ′ = −0.05 ±0.015 (1-sigma uncertainties were extracted from Fig. 4 of Bernstein et al. (2004)), consistent with our results. The increase in uncertainty of α and α ′ in our fit versus that from Bernstein et al. (2004) is caused by the marginalization in Eq. 21. We feel that our approach provides a more realistic estimate of the true uncertainties.
Additionally we fit a simple broken power-law model by Eq. 21. The model LF has a sudden change from the bright-side slope α 1 to the faint-side slope α 2 at a break magnitude m B , and is given by The best-fit parameters from maximizing the likelihood for the broken powerlaw LF while marginalizing over nuisance parameters C k and ∆m o,k as per Eq. 21 are m o = 23.2 ± 0.5, α 1 = 0.69 ± 0.08, α 2 = 0.57 ± 0.2, and m B = 24.4 ± 0.7.
The improvement of the maximum likelihood value when considering the rolling and broken power-law LFs given by Eqs. 23 and 24 fits is only a few percent over that of the single slope model. Hence the additional degrees of freedom included in the more complicated LFs compared to a single slope power-law does not substantially improve the fit and the higher order models are not statistically warranted in the range of magnitudes considered in the F07 sample. Indeed the best fits of both equations are consistent at the 1-sigma level with no break at all (α ′ = 0 and α 2 = α 1 ).
To test at what magnitude a roll over to shallower slopes would be inconsistent with the F07 sample, we made use of Monte-Carlo simulations and the Kolmogorov-Smirnov (KS) test. In these simulations, we simulated the observation of a broken power-law LF of the form of Eq. 24, with α 1 = 0.65 breaking to some faint-end slope α 2 at break magnitude m B . The simulated surveys and total number of detections were chosen to match the surveys and number of detections in the F07 sample. We included the effects of variable m o and colour offsets between surveys by randomly selecting these offsets for each individual survey.
For a given break magnitude and faint-end slope, a parent population of ∼ 10, 000 objects was generated. This parent population was simulated with a set of random C k and ∆m o,k offsets linearly sampled from the marginalization range described above (±0.2 for C k and ±0.4 for ∆m o,k ). Calibration and colour offsets occur in the observations in the F07 sample. The random sampling we implement provides a simple means of generating offsets consistent with the way in which those observations were made.
From the parent population, a sub-sample of objects was bootstrapped and the KS-statistic of this sub-sample compared to the parent population was calculated. This procedure was repeated for 1000 random sub-samples of the parent population. In this way, the KS statistic distribution was bootstrapped from the parent population.
The KS statistic of the F07 sample when compared to the parent population was compared to the distribution of KS-statistics generated from the bootstrapped samples. The roll-over model (α 2 , m B ) was rejected if 99.7% of the bootstrapped KS statistics were smaller than the KS statistic of the F07 sample. This was repeated 25 times using different random realizations of the C k and ∆m o,k offsets for each choice of faint-end slope and break magnitude. The repetition was necessitated by the randomness included with variable C k and ∆m o,k values; in certain circumstances, a particular set of colour and m o offsets will enhance the chances of observing a break. A model could not be rejected if the break was not rejected in 2 or more of the 25 simulation repetitions for that model. Fig. 8 is the 95% probability contour that a particular break model is rejected by the F07 sample. The contour has been smoothed to remove the effects of course sampling in α 2 -m B space. As expected, brightward of a particular magnitude most break models would likely have been detected in the observations. From these simulations, we conclude that the observations cannot reject the possibility of a break in the Kuiper belt LF with break magnitudes fainter than m B (R) ∼ 24.3. Bernstein et al. (2004) find that the LF is well described by the harmonic mean of a steep power-law for bright objects and a shallow power-law for faint objects with both power-laws contributing equally at magnitude R eq . They find R eq ∼ 22.8 − 23.6. This is not however inconsistent with the results of our KS test as m B and R eq are not equivalent parameters between the two LF models. The reader is cautioned about drawing conclusions from comparison of these two parameters.

Presented in
While no break is apparent in the observations, the best-fit single sloped model is inconsistent with the results from Bernstein et al. (2004) at m(R) 28 as the number of objects detected in that survey are too few by a factor of ∼ 6 from that predicted by the model. To account for this, the LF must roll-over to shallower slopes at some faint magnitude not present in the F07 sample.

The Size Distribution
Our best-fit luminosity function slope is α = 0.65±0.05. Under the assumption that Eq. 10 holds, we find that the KBO differential size distribution has a slope of q = 4.25 with a 1-sigma uncertainty of 0.25. If the belt is in a state of collisional equilibrium, we would expect the slope to be q ≃ 3.5. This is inconsistent with the inferred size distribution at the 3-sigma level. We conclude that, for objects larger than r ≈ 50 km, the size-distribution of the belt is inconsistent with a system in collisional equilibrium.

Discussion
To understand what the KBO size distribution tells us about the history of objects in the outer solar system, we must interpret our observations in terms of models that account for the size distribution of a belt of planetesimals. Generally, there are two broad types of models that attempt this, fragmentation models and accretion models.
Analytic fragmentation models are those in which a series of equations accounting for accretion and collisional disruption are solved, producing a steadystate collisional-cascade equilibrium. These models vary in complexity. Some assume that each disruption produces a number of equal size collision remnants (Dohnanyi, 1969;Pan and Sari, 2005), while some model object-disruptions with a distribution of collision remnant sizes (O'Brien and Greenberg, 2003). Some calculations also include various models of KBO physical strengths (O'Brien and Greenberg, 2003). The detailed results depend on the calculations, but generally the outcome is a differential size distribution in a quasiequilibrium steady-state with a slope equal to or less than Dohnanyi slope, q 3.5. Pan and Sari (2005), assumed that the strengths of large Kuiper belt objects are gravity dominated and modeled a collisional cascade in which all collisions were purely disruptive (no accretion). They found that, for objects smaller than some break size, the size distribution slope was q = 3 and could not account for the steep slope observed for large KBOs with this model. (2003) accounted for accretion and fragmentation, and a variation of object strength as a function of size. They found that the equilibrium size distribution slope was related to the power-law slope describing the variation in object disruption energy per unit mass with size. From this model, the Kuiper belt size distribution inferred from the F07 sample implies a unit disruption energy that decreases rapidly with increasing size. This is consistent with objects whose physical strength is dominated by tensile forces, but is inconsistent with objects whose strength is dominated by gravity; these objects have an increasing disruption energy with size. Objects in the tensile strength dominated regime are typically smaller than 1 km in size. The objects we observe here are likely to be gravity dominated and the scaling of KBO disruption energy versus size implied by the model from O'Brien and Greenberg (2003) seems unlikely. We conclude from these results that the observed LF and the inferred size distribution is inconsistent with analytic equilibrium models.

O'Brien and Greenberg
The large diversity of fragmentation models have similar outcomes; they predict a shallow size distribution slope of q ∼ 3.5 which is incompatible with the large object size distribution inferred from the luminosity function of the Kuiper belt.
Numerical accretion models are those in which the size distribution, and orbits of a population of bodies is calculated. Unlike analytic collisional cascade models, the size distribution calculated from accretion models is not assumed to be in steady state, but evolves in time along with the orbital distribution. Models which account for accretion and fragmentation of planetesimals in the region of the Kuiper belt, predict a broken power-law size distribution with a steep slope for large objects (D 10 km) that rolls over to a shallower slope for smaller objects (D ∼ 2 km) (Kenyon and Bromley, 2001;Kenyon, 2002;Kenyon and Bromley, 2004).
These models have two general evolutionary phases. The first phase is planetesimal accretion, in which a large reservoir of planetesimals on nearly circular orbits accrete to form larger bodies via low encounter-velocity collisions. This process rapidly produces a steep size distribution for large objects (q 5). Models from Kenyon and Bromley (2001) and Kenyon (2002) calculate planet growth in the 40-47 AU zone. These models have a relatively low-mass intial Kuiper belt, and do not include stirring from Neptune. They calculate planet growth before Neptune attains its current mass and orbit, and start from bodies with radius 1km. Initially, the size distribution of large objects is very steep (q 5). After ∼ 100 − 300 Myr, the slope of the largest objects flattens to q ∼ 4.1 − 4.4.
The second phase is started when the belt members are stirred up onto dynamically excited orbits (via interactions with Neptune, oligarchs, etc.) such that mutual collisions are mainly catastrophic. Collisional grinding rapidly reduced the slope of the size distribution for small bodies (D 1 km). The longer the collisional evolution continues, the larger the radius at which the break to shallower slopes occurs. In the models from Kenyon and Bromley (2001) and Kenyon (2002) the size distribution breaks to q 3.5 at D ∼ 2 − 20km after ∼ 1 Gyr.
More complicated models that start with a more massive initial Kuiper belt, and include effects of gas-drag and a more detailed model of object strength, produce large object slopes of q ∼ 3 and after 4.5 Gyr, and a break at D ∼ 2 − 40 km to very flat slopes of q ∼ 0 − 0.5 for small objects. (Kenyon and Bromley, 2004) The observed size distribution exhibits a steeper slope than those determined by Kenyon and Bromley (2004) and is more consistent with models that start with a less massive Kuiper belt. In these models, the break radius grows larger and the large object slope flattens as accretion and collisional evolution continues. An accurate measure of the break size in addition to the large object slope determined here, combined with modeling, may further constrain the duration of accretion in the region of the Kuiper belt.
Accretion models which include migration of Neptune can account for some of the features of the KBO orbit distribution, but come at the cost of complicating the accretion process (Kenyon et al. 2007 and references there in). In these models, most of the current Kuiper belt consists of objects that are scattered by Neptune to their current positions, and originate from regions interior to Neptune's current orbit. Accretion of the KBOs in their original locations proceeds until stirring and scattering by Neptune disrupts the accretion process.
In these migration models the KBOs are initially closer to the Sun than there current locations and accrete in regions of high density, and therefore both the accretion and collision process occur on much faster time-scales than in-situ formation models predict. The migration time-scale becomes an important parameter in these models; too fast a migration and the largest KBOs are not formed, while too slow, and the break size can evolve to be large enough such that it would be detected in the available observations.
Numerical accretion models of KBO formation, are generally in better agreement with the observed large size distribution than collisional cascade equilibrium models. Accretion models suggest that knowledge of the size distribution slope, and the radius at which the break to shallower slopes occurs will provide strong constraint on the accretion dynamics, and KBO formation time-scales. Due to the large degree of complexity in these models however, there is much work needed to interpret the observed size distribution. Models which achieve collisional-cascade equilibrium for large objects cannot account for the observed steep size distribution.

Conclusions
We have performed a survey of the Kuiper belt with an sky coverage of 3.0 square degrees with our deepest field having a depth of m(g ′ ) = 26.43 magnitudes. An analysis of current survey data confirms that the luminosity function of the belt is well described by a single slope power-law between m(R) = 21−26 with a slope of α = 0.65 ± 0.05 and normalization m o = 23.42 ± 0.13 which is typical of all the fields considered in the F07 sample.
We have shown the necessity of considering calibration, and density effects when inferring the luminosity function from the combination of different surveys that have observed separate regions of the sky, and are not directly calibrated to one-another. Such effects can skew the inferred LF away from the true form, and cause uncertainties in the inference to be highly underestimated.
We conclude that: (1) The Kuiper belt is not in a state of collisional equilibrium for objects larger than D ≈ 50 km.
(2) There is no evidence for a change in the slope of the luminosity function for magnitudes brighter than m(R) ∼ 24.3 corresponding to objects with diameters D 110 km at 40 AU. The sample of observations considered here is consistent with a best-fit model that breaks at m(R) ∼ 24.4 to a slope of α 2 ∼ 0.6, but the more complicated fit is not warranted by the slight improvement in the maximum likelihood value. (3) The observed slope of the luminosity function implies that the size distribution of the Kuiper belt is consistent with a power-law with slope q = 4.25 ± 0.25 for objects with diameters larger than D ∼ 50km.

Acknowledgements
We thank the reviewers for providing excellent feedback and suggestions for improvement of this manuscript. We also thank Tom Loredo for his many discussions in the use of statistics.
This project was funded by the National Science and Engineering Research Council and the National Research Council of Canada. This research used the facilities of the Canadian Astronomy Data Centre operated by the National Research Council of Canada with the support of the Canadian Space Agency.

A Psfmatch3
Psfmatch3 is a subtraction routine that compensates for variations between the image quality of the template image and that of individual images. A description of the basic algorithm follows.
Consider an observed image O ij (i, j refer to pixel row and column), and a suitably chosen reference image I ij . We define an n k × n k convolution kernel K ij such that the convolution K * I provides some "best" (in a least-squares sense) estimator,Õ, of O; that is, the optimal kernel K is the one that minimizes i,j (O ij −Õ ij ) 2 . Differentiating with respect to each kernel element K ij yields a system of M = n 2 k linear equations which can be solved for the n 2 k kernel elements. This basic procedure can be improved considerably by adding an extra term for sky differences, viz.Õ = K * I + ∆s, where ∆s is a constant denoting the difference in sky between the two images. Weighting by errors is easily included by minimizing Perhaps most important is allowing for spatial variations in the kernel (and sky background); this is accomplished simply by solving for polynomial coefficients representing the spatial variation of each kernel coefficient.
There are a number of features and advantages to this method of matching the point spread function on two images.
• The kernel K has arbitrary form; it does not need to be symmetric, and can handle PSF variations that may not always be handled by methods involving a basis set of functions to represent the kernel. The kernel can even be solved to perform deconvolution (which is the case when the reference image is erroneously chosen to have worse seeing), though this results in noise amplification. • K is not necessarily normalized to unity; this automatically takes out transparency fluctuations. • The method automatically removes small shifts (even spatially variable shifts) between two images. • As noted above, the method can easily be adapted to include background differences, spatially variable backgrounds, and spatially variable kernels.   TE1G 476727o, 476728o, 476729o, 476848o, 476849o, 476850o, 476851o, 476852o, 476853o, 476854o, 476855o, 476856o, 476857o, 476858o, 476859o, 476984o, 476985o, 476986o, 476987o, 476988o, 476989o, 476990o, 476991o TE1M 476717o, 476718o, 476719o, 476720o, 476721o, 476992o, 476993o, 476994o, 476995o, 476996o TE2G 476885o, 476886o, 476887o, 476888o, 476889o, 476890o, 476891o, 476892o, 476893o, 476894o, 476895o, 476896o,476924o TE3G 476758o, 476759o, 476760o, 476761o, 476762o, 476763o, 476764o, 476765o, 476766o, 476767o, 476768o, 476769o,476795o, 476796o, 476797o, 527174o, 527305o, 527458o, 527461o TE3M 476798o, 476799o, 527175o, 527306o, 527307o, 527459o, 527460o TE4G 502047o, 502048o, 502049o, 502183o, 502184o, 502215o, 502217o, 502218o TE4M 502050o, 502051o, 502052o, 502182o, 502185o, 502186o, 502214o, 502216o, 502374o TE5G 502102o, 502139o TE5M 502098o, 502099o, 502100o, 502101o, 502103o, 502136o, 502137o, 502138o, 502140o, 502248o, 502249o, 502250o, 502426o   Table 7: Results of the maximum likelihood fits for all fields fit individually. Uncertainties are taken from the extrema of the 1-sigma likelihood contours. *-Provides a 3-sigma lower limit only. N -Number of objects detected in the field. N 50 -number of objects brightward of 50 % threshold.   Rates are in arcseconds per hour. Angles are in degrees below the horizontal. The source visible in the images is a real 23.8 mag. object with a rate of motion close to 2 arcsec. hr −1 at 20.8 • , consistent with the ecliptic. The subtraction well and characteristic trailing can be seen at other rates of motion. Fig. 3 Distribution of inclinations and distances for the artificial objects planted in the observations, taken from the planted objects lists of 10 chips in the MEGAPrime field. Solid circles mark those planted objects that were found during the object search. Open circles mark those that were not found. The same magnitude distribution of artificial objects was used at all distances. Thus, for a given magnitude, the detection efficiency did not depend on object distance. Fig. 4 Net efficiency for each field with 50% efficiency marked with a vertical line. Magnitudes are converted to R-band using average colours (see Section 6). Points are the binned efficiencies determined from the image search. Errorbars are 1-sigma Poisson confidence regions. The solid curve is the best fit efficiency using Eq. 2. The dotted curve is a cumulative histogram of visually rejected false candidates. Each field is labeled; a: MEGAPrime N1, b: UNE, c: UNW, d: N10032W3, e: N11033, f: NEP0813NW3, g: NEP0815NE3.   Fig. 8 Probability of rejection of a broken power-law from the data. Plotted is faint-end slope α 2 versus break magnitude m B . Contour represents the 95% rejection confidence region. Contour has been smoothed to remove effects of course simulation sampling.