text:
# Beware the Non-uniqueness of Einstein Rings

## 1 Introduction

It seems likely now that all multiple-image lens systems have Einstein rings, if imaged deep enough. Radio rings, formed out of multiply-imaged radio lobes and jets, have been known for more than a decade (Hewitt et al. 1988, Langston et al. 1990, Lehar et al. 1993). But more recently, QSO lenses when deeply imaged—especially in the infrared—have shown rings formed out of multiply-imaged parts of the QSO host galaxy: Falco et al. (1997) examined HST images of MG0414+0534 and found an arc connecting three of the four QSO images; Courbin et al. (1997) found an incipient ring in 1115+080; Impey et al. (1998) uncovered the full ring. Blandford, Surpi, & Kundić (2000, hereafter BSK) identified ring components in 1608+656.

Such a rich source of new image data is very exciting, because the paucity of constraints is a very serious problem in lens modeling, leading to large uncertainties in current lens reconstructions (see e.g., Bernstein & Fischer 1999, Williams & Saha 2000, hereafter WS). BSK and Kochanek, Keeton, & McLeod (2001, hereafter KKM) present new theoretical ideas for using ring observations to constrain lenses. In particular KKM make two very interesting points (see especially their Figure 1): (i) in a double (quad), the ring is the image of a two (four) lobed figure on the source; (ii) where a ring has brightness minima, it intersects a critical curve.

This paper presents some further theoretical work, in which we address the following question: for a lens which already has accurate image positions and time-delays, how much does a ring further constrain the lens? Our approach is very different from previous ones, and leads to some new insights: the down side is that the same ring can arise from a large range of lenses and $`h`$ values, and thus the prediction by KKM that “deep, detailed observations of Einstein rings will be revolutionary for constraining mass models and determining the Hubble constant from time delay measurements” is, unfortunately over-optimistic; the up side is that rings turn out to be related to time delay measurements in an unexpected but testable way.

## 2 Rings and the arrival-time surface

While working on WS we noticed that a contour map of the arrival-time surface in a multiple-image system resembles a ring if the contours are very close. Figure 14 of WS shows an example. In that paper we gave a handwaving argument for why this should happen, but returning to the question afterwards led to interesting developments, which we now discuss.

The arrival-time surface produced by a source at $`\theta _\mathrm{s}`$ can be expressed as

$$\tau (\theta )=2_\theta ^2\left[1\kappa (\theta )\right]\theta \theta _\mathrm{s}$$
(1)
with $`\tau (\theta )`$ being related to the physical arrival time111We will use ‘arrival-time surface’ to refer to both $`\tau (\theta )`$ and $`t(\theta )`$, in places where this causes no ambiguity.

$$t(\theta )=H_0^1(1+z_\mathrm{L})\frac{d_\mathrm{L}d_\mathrm{S}}{d_{\mathrm{LS}}}\tau (\theta ),$$
(2)
where the $`d`$’s are dimensionless angular diameter distances, $`z_\mathrm{L}`$ is the lens-redshift, and $`\kappa (\theta )`$ the dimensionless mass profile (see e.g., equations 9 and 10 in Saha 2000). If the source is much further than the lens, then

$$t(\theta )h^1z_\mathrm{L}(1+z_\mathrm{L})\times 80\mathrm{days}\mathrm{arcsec}^2\tau (\theta ).$$
(3)
If image positions, QSO time delays, redshifts, cosmology, and lens profile are all known, $`H_0`$ can be read off from equations (1) and (2). Unfortunately, this is not the case; image positions, time delays, and redshifts are accurately measurable, and cosmology produces a readily-quantifiable (and fairly small) uncertainty, lens profile is a large uncertainty: a large range of lens models $`\kappa (\theta )`$ can reproduce QSO image positions and time delay ratios exactly or nearly exactly, resulting in a range of plausible $`H_0`$’s.

An Einstein ring, being an image of an extended source, is not determined solely by the arrival-time surface for a point source, but also by the brightness distribution of the source—see e.g., equation (7) of KKM. KKM and BSK argue that not all the mass models that reproduce the QSO images will be able to reproduce the image of the QSO host galaxy. If so, it would narrow down the range of acceptable $`H_0`$’s.

How does an extended source change things? To study this, let us consider an extended source with a very simple brightness profile: peak brightness at $`\overline{\theta _\mathrm{s}}`$ and around it source brightness falling off in a conical profile

$$I_\mathrm{s}(\overline{\theta _\mathrm{s}}+\mathrm{\Delta }\theta _\mathrm{s})=1k\left|\mathrm{\Delta }\theta _\mathrm{s}\right|.$$
(4)
Here $`k`$ is a constant, and it is understood that where $`1k\left|\mathrm{}\right|<0`$, the surface brightness is zero. Say one of the images of $`\overline{\theta _\mathrm{s}}`$ is at $`\overline{\theta }`$ on the lens plane, and that the image-plane region $`\overline{\theta }+\mathrm{\Delta }\theta `$ corresponds to the source-plane region $`\overline{\theta _\mathrm{s}}+\mathrm{\Delta }\theta _\mathrm{s}`$. Since surface brightness is preserved, we have for the image brightness

$$I(\overline{\theta }+\mathrm{\Delta }\theta )=I_\mathrm{s}(\overline{\theta _\mathrm{s}}+\mathrm{\Delta }\theta _\mathrm{s})=1k\left|\mathrm{\Delta }\theta _\mathrm{s}\right|.$$
(5)
A small displacement in the source plane, $`\mathrm{\Delta }\theta _\mathrm{s}`$, is related to the corresponding displacement in the lens plane, $`\mathrm{\Delta }\theta `$, through $`\mathrm{\Delta }\theta _\mathrm{s}=𝐌^1\mathrm{\Delta }\theta `$, where $`𝐌`$ is the inverse magnification. With this, the illumination in the lens plane can be rewritten as

$$I(\overline{\theta }+\mathrm{\Delta }\theta )=1k\left|𝐌^1\mathrm{\Delta }\theta \right|.$$
(6)
Now $`\tau (\theta )`$, the arrival time for the point source at $`\theta _\mathrm{s}`$, has zero derivative at $`\overline{\theta }`$. In the neighborhood of $`\overline{\theta _\mathrm{s}}`$, the second derivatives of $`\tau `$ are simply $`𝐌^\mathrm{𝟏}`$. Hence in this region, $`\tau `$ equals $`𝐌^1\mathrm{\Delta }\theta `$. So the brightness of the ring is

$$I(\theta )=1k\left|\tau \right|.$$
(7)
Now picture a sheet of paper with contours of arrival time drawn on it, viewed so that nearby contour lines blur into a grayscale. Equation 7 then also happens to represent the grayscale pattern on this sheet of paper, where

$$k\frac{\text{thickness of contour lines}}{\tau \text{-spacing between contour lines}},$$
(8)
with ‘1’ representing blank paper, and $`k\left|\tau \right|`$ the darkening due to the contour lines. We thus have a quantitative explanation for why Figure 14 of WS looks like a ring.

More importantly, we see that given a lens model for QSO images, generating a ring image for a conical source is very simple: one has only to plot a closely spaced contour map of the arrival-time surface; the only adjustable parameter is the contour interval, which behaves as a surrogate for the source size $`k^1`$.

The relation (7) between the arrival-time surface (of the QSO or other compact source) and the ring (formed from the host galaxy) provides two insights, which are quite generic and not limited to conical sources.

1. Consider time delays versus image separations between different pairs of images in a quad. Where the ratio of time delay to image separation is small, the density of contour lines will be small, therefore the corresponding segment of the ring will be brighter. Thus, if time delay ratios in a quad have been measured then the relative brightness of different ring segments can be predicted to some extent. Conversely, if QSO images and the image of the ring are observed then time delay ratios can be estimated.
2. Altering the lens mass model, and hence $`\tau (\theta )`$ in such a way that there is little or no change in the density of contour lines in the region of the ring will produce little or no change in the ring (other than its overall brightness), but may cause a large change in $`H_0`$. One example is evident from equations (1) and (2): if $`[1\kappa (\theta )]`$, $`\theta _\mathrm{s}`$, and $`H_0`$ are all multiplied by the same constant factor, $`t(\theta )`$ is unaffected. This is the mass-disk degeneracy, first noted by Falco et al. (1985) and discussed in several recent papers, e.g., Saha (2000); basically it means that steeper lens profile give higher $`H_0`$ from the same data—see the “radial-index versus $`h`$” plots in WS for examples. This particular degeneracy is exact, but numerically one can cook up any number of near-degeneracies that scale $`H_0`$ significantly while producing only slight, observationally negligible, changes in rings.

In the following section we follow up the implication of these two points in detail for two quads.

## 3 Case studies: 1115+080 and 1608+656

These are among the most-studied lenses: 1115+080 was discovered by Weymann et al. (1980) and its time delays were determined by Schechter et al. (1997); 1608+656 was discovered by Myers et al. (1995) and Fassnacht et al. (1999) measured its time delays.

Figure 1 shows the image configurations for both lenses. We have plotted the saddle-point contours, which make up a sort of skeleton for the arrival-time surface, and renamed the images by arrival time. Upon doing this, it becomes evident how similar these two lenses are.

We can now use the ideas of the previous section to predict the morphology of the rings. The contour-density argument show that

$$\frac{(\mathrm{\Delta }\tau _{13}/\mathrm{\Delta }\theta _{13})}{(\mathrm{\Delta }\tau _{24}/\mathrm{\Delta }\theta _{24})}\frac{\text{brightness of 2–4 ring segment}}{\text{brightness of 1–3 ring segment}},$$
(9)
and similarly for any pair of ring segments. Now, arrival-time contours will be furthest apart in the region 2–3, closer together in 1–3 and 2–4, and closest in 1–4. So the brightest part of the ring will be 2–3, while 1–4 will be a faint arc or a gap in the ring. The same conclusions follow from a more complicated argument based on caustics in Blandford & Kovner (1988). And they may be verified by inspecting Figure 1 of Impey et al. (1998) and Figure 1 of BSK.

We now show the predicted rings from models that precisely fit the QSO image and time-delay data. We generate ensembles of such models using the method developed in WS, with one change: instead of letting $`h`$ vary in the ensemble and then examining its range, we fix it at some trial value; in this way we can examine how the predicted ring varies with the trial $`h`$. All models we show in this paper are averages of ensembles of 100. We have not attempted to fit for the contour step, we have just chosen a suitable value empirically.

Figure 2 shows a model of PG1115, for $`h`$ fixed at 0.5. The model is constrained by image positions and Barkana’s (1997) estimate of the time delay from the Schechter et al. (1997) data. It predicts the gross features of the observed ring. Figure 3 is similar ($`h`$ is also set to 0.5), but uses the earlier preliminary time delay estimate by Schechter et al. (1997); the only significant difference is that $`\mathrm{\Delta }t_{2,4}/\mathrm{\Delta }t_{1,3}`$ is larger, and in response the ring segment 2–4 has become fainter than 1–3. In the observed ring the 2–4 segment is actually slightly brighter than 1–3. (See Figure 1 of Impey et al. 1998, and also Figure 3 of KKM). Thus, the observed ring favors the Barkana estimate of $`\mathrm{\Delta }t_{2,4}/\mathrm{\Delta }t_{1,3}`$. Moreover, since the model is still slightly under-predicting the brightness of 2–4 versus 1–3, we further predict that as the time-delay measurements improve, $`\mathrm{\Delta }t_{2,4}/\mathrm{\Delta }t_{1,3}`$ will get somewhat smaller. This illustrates point 1 from the end of Section 2, that QSO image positions, time delay ratios, and the ring are not all independent.

The above conclusions do not depend on the trial value chosen for $`h`$. Figure 4 shows a model using the Barkana values and $`h=0.75`$. The ring is similar to that in Figure 2 (Barkana values and $`h=0.5`$), and demonstrates that QSO image positions and time delays determine the ring, while the mass distribution can vary substantially: compare the top right mass map in Figure 4 to that in Figure 2. This illustrates point 2 from the end of Section 2. It is also contrary to the assertion by KKM that “The Einstein ring in PG1115+080 rules out the centrally concentrated mass distributions leading to a high Hubble constant ($`H_0>60\mathrm{km}\mathrm{s}^1\mathrm{Mpc}`$) given the measured time delays”.222KKM’s limit actually uses the preliminary Schechter et al. value for $`\mathrm{\Delta }t_{2,4}/\mathrm{\Delta }t_{1,3}`$, rather than the later Barkana value. Our argument would not change if we used the other value. Incidentally, note that the $`h=0.75`$ model in Figure 4 is much more elliptical than the $`h=0.5`$ model in Figure 2, but it is not much more centrally concentrated. It is easier to produce $`h>0.6`$ models for this system if they are more centrally concentrated, but it is not necessary (see Figure 13 of WS).

Figures 5 and 6 are for 1608+656, and are analogous to Figures 2 and 4. The data are from Fassnacht et al. (1999) and the predicted rings may be compared with Figure 1 of BSK. Again, we find that changing $`h`$ makes little difference to the ring.

Why do we, though we have no disagreement with the basic theoretical ideas in KKM, come to the opposite conclusion regarding constraints on $`h`$? The reason lies in the modeling strategy. KKM assume a fixed functional form for the lens described by a few parameters; these models have little freedom to vary radial profile, and do not allow ellipticity to vary or to twist with radius. Pixelated models, such as the ones in this paper, do not have such restrictions and thus are able to expose the degeneracies left by rings. Of course it is not necessary to use pixelated models for this purpose, parametric models with some additional parameters would also be effective. What these additional parameters need to be is not obvious, but from examining the pixelated mass maps we guess that letting the index of the radial profile vary and letting the ellipticity vary with radius are probably the most important ingredients.

## 4 Conclusions

We show that the arrival-time surface for a multiply-imaged QSO (or other compact source), when plotted as a contour map with densely packed contours, is approximately the lensed image of an extended source with a conical surface-brightness profile. An observed ring will depend on the source’s detailed brightness profile, so our result does not furnish a method for making detailed models of rings. However, our result does imply that the gross features of a ring are determined by the information encoded in the QSO point images alone, which leads to two useful findings, as follows.

1. It becomes possible to translate statements about time delays into statements about the ring (eq. 9). In particular, the general morphology of rings in quads becomes very easy to predict, from just the configuration of the QSO images. For the case of 1115+080, the ring can even be used to predict the sign of the current errors in the ratio of two time delays.
2. It becomes evident that the same ring can arise from a large range of mass distributions, and hence from a large range of $`h`$. The main (but not the only) contributor to this non-uniqueness is the mass disk degeneracy, which makes the lens mass profile steeper while making the Universe smaller but keeps image and time-delay data the same. We illustrate non-uniqueness in models of 1115+080 and 1608+656.

We thank J.C.B. Papaloizou for elucidating the connection between contour maps and ring images.


format: text

text:
# An Interface View of Directed Sandpile Dynamics

## Abstract

We present a directed unloading sand box type avalanche model, driven by slowly lowering the retaining wall at the bottom of the slope. The avalanche propagation in the two dimensional surface is related to the space-time configurations of one dimensional Kardar-Parisi-Zhang (KPZ) type interface growth dynamics. We express the scaling exponents for the avalanche cluster distributions into that framework. The numerical results agree closely with KPZ scaling, but not perfectly.

Avalanche phenomena are common in nature . They are characterized by fast relaxation dynamics under a slow driving force. Models that describe such dynamics, e.g., so-called sandpile models, have been studied extensively for more than a decade following the work of Bak *et al* . Directed sandpile models are a special sub-class in which relaxation follows a directional rule, that is, the propagation of active sites occurs only in one direction and never backfires. The central issue in this type of research is whether the dynamics is critical, such that the avalanche distribution functions are scale invariant (power-law decay), and if so, whether these scaling properties are universal in the same sense as equilibrium critical phenomena. Our understanding of these issues is still restricted. There are only a few exactly soluble models, e.g., some deterministic Abelian directed sandpile models , but most insight is still limited to numerical simulation data. The evidence for scaling and universality in other types of non-equilibrium dynamics is less ambiguous: In surface growth (another example of intrinsic critical behavior) several universality classes are well established; e.g., Edwards-Wilkinson (EW)  and Kardar-Parisi-Zhang (KPZ)  type surface growth; population and catalysis type dynamics undergo absorbing state type phase transitions with distinct universality classes like directed percolation and directed Ising .

Efforts are under way to link avalanche dynamics to these better understood dynamic phenomena. This ranges from mappings to driven interfaces  and directed percolation  to using the concepts of renormalization  and universality classes . In this letter we introduce a two dimensional (2D) directed avalanche model in which the scaling exponents of the avalanche dynamics are directly related to the exponents of KPZ growth of a 1D interface.

The physical system we have in mind is a sand box with a movable retaining wall to let out sand from the bottom of the slope, see Fig.(1). The retaining wall is lowered very slowly, such that grains tumble out sporadically forming distinct avalanches instead of a continuous flow. We model the sand surface by continuous height variables $`h(x,y)`$, with respect to a 2D square lattice, which is rotated over $`45^{}`$, meaning that in the even (odd) $`y`$ rows $`x`$ takes only even (odd) integer values. The retaining wall is placed at the $`y=0`$ boundary. The slope is stabilized by the following constraint. The surface particle in column $`(x,y)`$ is supported by the two columns $`(x\pm 1,y1)`$ just below it, and must be lower than the lowest of the two increased by an amount $`s_c`$,

$$h(x,y)\mathrm{min}[h(x+1,y1),h(x1,y1)]+s_c.$$
(1)
An avalanche is triggered by selecting the highest site, $`(x_i,0)`$, at the $`y=0`$ wall boundary (it is the $`i`$-th avalanche) and reducing its height by a random amount, $`0<\eta _i<s_c`$. This represents the lowering of the retaining wall. Next, all sites that violate the stability condition topple according to the rule,

$`h(x,y)`$ $``$ $`\mathrm{min}[h(x+1,y1),h(x1,y1)]`$ (2)
$`+`$ $`\eta _i(x,y),`$ (3)
where $`0<\eta _i(x,y)s_c`$ are uncorrelated random numbers. This toppling continues until the whole system is stable again. Since the toppling of a site in row $`y`$ can only effect the stability of two sites immediately above it in row $`y+1`$, the sites can be updated row by row starting from the $`y=0`$ boundary.

This process is idealized compared to a real unloading sand box in the sense that the toppled grains drop out without disturbing the already stabilized lower regions of the surface. It is possible to justify this as the low gravity or strong bond limit where the falling sand does not gain enough momentum to disturb the stabilized surface on its way out. Without this idealization we loose the directed nature of the dynamics and the following surface growth interpretation.

The row-by-row toppling sequence (3) can be reinterpreted as a dynamic rule for a 1D growing interface, in which the $`y`$ coordinate plays the role of time. Each stable surface configuration represents a space-time configuration of the 1D interface. Imagine creating an initial stable surface configuration, before the retaining wall starts to drop: choose an arbitrary configuration with all $`0<h(x,0)s_c`$ in row $`y=0`$. Next, apply Eq.(3) to all sites in the next row, $`y=1`$, to create the next slice of the surface. Repeat this for all higher rows. The configuration in every row is like a 1D interface evolving in time $`t=y`$. Fig.2 illustrates how this interface propagates during each time step, $`yy+1`$. The upper panel shows the first half of the update (the drawn to the dashed line). This is the deterministic part of the propagation (the $`\mathrm{min}[]`$ operator) in Eq.(3). The lower panel illustrates the second half of the update, where the heights increase by a random amount $`0<\eta s_c`$. The first step removes material, and the second step deposits particles.

This type of interface dynamics almost certainly belongs to the KPZ universality class. Eq.(3) can be rewritten as

$`h(x,t)`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left[h(x+1,t1)+h(x1,t1)\right]`$ (4)
$``$ $`{\displaystyle \frac{1}{2}}\left|h(x+1,t1)h(x1,t1)\right|`$ (5)
$`+`$ $`\eta (x,t),`$ (6)
which is a discrete version of the KPZ Langevin equation ,

$$\frac{h}{t}=^2h\frac{\lambda }{2}\left(h\right)^2+\eta ,$$
(7)
To be absolutely sure, and also to make sure that $`\lambda `$ is large enough such that corrections to scaling from the EW point (at $`\lambda =0`$) do not interfere, we checked numerically the behavior of the surface width $`\mathrm{\Delta }(L,t)`$, defined as $`\mathrm{\Delta }^2\left(hh\right)^2`$. The scaling properties of KPZ growth are known exactly in 1+1 dimensions. Starting from a flat surface at $`y=0`$, the width increases as $`\mathrm{\Delta }t^\beta `$ for $`0tL^z`$ and saturates at $`\mathrm{\Delta }L^\alpha `$ for $`tL^z`$, with $`L`$ the lattice size in the $`x`$ direction. The exponents are exactly equal to $`\alpha =\frac{1}{2}`$, $`\beta =\frac{1}{3}`$ and $`z=\alpha /\beta =\frac{3}{2}`$. Our numerical results are shown in Fig.(3) and are consistent with the KPZ values.

Throughout this paper we present numerical results in terms of finite size scaling plots of effective exponents, like $`\alpha (L)`$ in Fig.(3)(a). Global straight line fits to log-log plots of, e.g., $`\mathrm{\Delta }L^\alpha `$ are inaccurate. $`\alpha (L)`$ is obtained from fitting the form $`\mathrm{\Delta }L^\alpha `$ to two nearby values of $`L`$. The approach to $`L\mathrm{}`$ in Fig.(3) is consistent with a leading corrections to scaling exponent $`y=\frac{1}{2}`$.

The characteristic feature of self-organized criticality (SOC) is the lack of typical avalanche length, width, depth, or mass scales. The probability distributions follow power-laws, like $`P_ww^{\tau _w}`$ characterized by the scaling exponents, $`\tau _l`$, $`\tau _w`$, $`\tau _\delta `$ and $`\tau _m`$. Our numerical simulation results confirm the existence of scale invariance. The critical exponents converge well, see Fig.(4). The length $`l`$ is the maximum $`y`$ coordinate the avalanche reaches. The width $`w`$ is the maximum departure in the $`x`$ direction, $`\left|xx_i\right|`$, from the trigger point $`x_i`$. The depth $`\delta `$ is the maximum height change, $`h_ih_{i1}`$, caused by the avalanche, and the mass is the total amount of sand carried off by the avalanche.

The meta-distribution function, $`P(l,w,\delta )`$, should obey the scaling form

$$P(l,w,\delta )=b^\sigma P(b^zl,b^1w,b^\alpha \delta ).$$
(8)
with $`b`$ an arbitrary scale parameter. The exponents, $`\sigma `$, $`z`$ and $`\alpha `$, are expected to be robust with respect to details of the dynamic rule, and thus define the universality class. Single parameter distributions, like $`P_ww^{\tau _w}`$, follow by integrating out the other two parameters. This implies the following expressions for the $`\tau `$ exponents

$$\tau _l=\frac{\sigma 1\alpha }{z},\tau _w=\sigma z\alpha ,\tau _\delta =\frac{\sigma 1z}{\alpha }.$$
(9)
or inverted,

$$z=\frac{\tau _w1}{\tau _l1},\alpha =\frac{\tau _w1}{\tau _\delta 1},\sigma =\tau _w+z+\alpha .$$
(10)
In Fig.(5) we replot the numerical finite size scaling estimates of the $`\tau `$ exponents in terms of $`\alpha `$, $`z`$, and $`\sigma `$.

The values $`z=1.52\pm 0.02`$ and $`\alpha =0.46\pm 0.01`$ are very close to those of 1D KPZ growth. Every stable slope configuration represents a possible world sheet of a 1D KPZ type interface, and every avalanche the difference between two such world sheets. Therefore it is natural to expect that the length (depth) of the avalanche scales with the KPZ value for $`z`$ ($`\alpha `$).

The distribution of avalanche cluster sizes is the most commonly studied and experimentally the most accessible property of SOC. Its exponent is linked to $`\alpha `$ and $`z`$ in the following manner. At the start of the avalanche, the height of a boundary site ($`y=0`$) is lowered on average by $`\frac{1}{2}s_c`$. Thus, for a sandbox of width, $`L_x`$, the boundary row is lowered by $`\frac{1}{2}s_c`$ after $`L_x`$ avalanches. In the stationary state, the entire surface moves down on average by $`\frac{1}{2}s_c`$ and the average amount of removed sand is equal to $`L_y\times \frac{1}{2}L_xs_c`$. Thus, the average mass of an avalanche must be equal to

$$m=𝑑mmP_m(m)=\frac{1}{2}s_cL_y.$$
(11)
This is analogous to conservation of current in conventional deposition type avalanche systems (see, e.g., ). Assume that the avalanche clusters are compact, i.e., that the sizes of holes of unaffected regions inside an avalanche do not scale with the avalanche size. In that case, the mass scales as $`ml\times w\times \delta `$, and we can use the meta-distribution function to evaluate Eq.(11), as

$`m`$ $``$ $`{\displaystyle _0^{L_y}}𝑑l{\displaystyle _0^{\mathrm{}}}𝑑w{\displaystyle _0^{\mathrm{}}}𝑑\delta lw\delta P(l,w,\delta )`$ (12)
$`+`$ $`m_{L_y}{\displaystyle _{L_y}^{\mathrm{}}}𝑑l{\displaystyle _0^{\mathrm{}}}𝑑w{\displaystyle _0^{\mathrm{}}}𝑑\delta P(l,w,\delta )`$ (13)
This applies when the box is wide and deep enough that $`L_y`$ is the only limiting factor to the avalanches. The first term accounts for all avalanches that fit inside the box, and the second term for the ones that reach the $`L_y`$ edge and thus are prematurely terminated. The first integral scales as $`L_y^{(\sigma +2+2z+2\alpha )/z}`$ for large $`L_y`$, while the second one scales as $`L_y^{(\sigma +1+z+\alpha )/z}`$. We have assumed $`ml\times w\times \delta `$ so $`m_{L_y}L_y^{(1+z+\alpha )/z}`$. The two terms in Eq.(13) scale in the same way, as

$$mL_y^{(\sigma +2+2z+2\alpha )/z}.$$
(14)
Solving Eqs.(11) and (14) for $`\sigma `$ gives

$$\sigma =2+z+2\alpha $$
(15)
Numerically we find $`\sigma =4.43\pm 0.05`$, see Fig.(5), in agreement with the numerical values for $`\alpha `$ and $`z`$, and also with their KPZ values (which yields $`\sigma =9/2`$).

In conclusion, we introduced a directed unloading sand box model, in which the stable slope configurations obey KPZ type scaling, and the avalanches represent the difference between two such KPZ world sheets. The avalanche distribution exponents can be reformulated in KPZ language, Eq.(10), and the numerical results agree with the KPZ values.

An intriguing issue for further study is whether the slight systematic deviations from the KPZ values in Fig.(5) are for real or just an artifact of, e.g., finite size scaling (the avalanches cover only a small part of the surface). In KPZ dynamics the scaling exponents follow from an ensemble average over independent MC runs, i.e., over a large set of independent world sheets. The avalanche dynamics performs this ensemble average in a correlated manner. All subsequent world sheets are identical except for the avalanche area. Correlated MC runs like this might well change the scaling exponents (but only slightly apparently), and this might prove a key feature of avalanche type SOC. The dashed curve in Fig.3(b) shows initial results for the global surface roughness in such avalanche correlated MC runs. The surface is noticeably rougher in amplitude but the finite size scaling values of the exponent $`\beta `$ are somewhat smaller (consistent with the values for $`\alpha `$ and $`z`$ in Fig.5). Whether they will converge in the end to a value $`\beta <\frac{1}{3}`$ remains to be seen.

In this context it is noteworthy that Eq.(10) applies also to recent solutions to the stochastic directed sandpile model (SDSM) of Paczuski et.al.  and Kloster et.al.  where the avalanche dynamics is related to Edwards-Wilkinson (EW)  surface growth, with $`z=2`$ and $`\alpha =1/2`$; and also that the exactly soluble directed sandpile model (DSM) of Dhar and Ramaswamy satisfies Eq.(10) with $`z=2`$ and $`\alpha =0`$. However, that does not resolve the issue, because both models are intrinsically simpler than ours, and the numerical accuracy of SDSM in, e.g., Ref. is not sufficient to detect a small effect like this.

This research is supported by the National Science Foundation under grant DMR-9985806.


text:
# PROPERTIES OF THE TRANS-NEPTUNIAN BELT: STATISTICS FROM THE CFHT SURVEY1footnote 11footnote 1Based on observations collected at Canada-France-Hawaii Telescope, which is operated by the National Research Council of Canada, the Centre National de la Recherche Scientifique de France, and the University of Hawaii.

## 1 Introduction

The rate of discovery of Kuiper Belt Objects (KBOs) has increased dramatically since the first member (1992 $`\mathrm{QB}_1`$) was found (Jewitt & Luu 1993). As of Dec 2000, $`400`$ KBOs are known. These bodies exist in three dynamical classes (Jewitt, Luu & Trujillo 1998): (1) the Classical KBOs (CKBOs) occupy nearly circular (eccentricities $`e<0.25`$) orbits with semimajor axes $`41\text{ AU}a46\text{ AU}`$, and they constitute $`70`$ % of the observed population; (2) the Resonant KBOs occupy mean-motion resonances with Neptune, such as the 3:2 (the Plutinos, $`a39.4`$ AU) and 2:1 ($`a47.8`$ AU), and comprise $`20`$ % of the known objects; (3) the Scattered KBOs represent only $`10`$ % of the known KBOs, but possess the most extreme orbits, with median semimajor axis $`a90`$ AU and eccentricity $`e0.6`$, presumably due to a weak interaction with Neptune (Duncan & Levison 1997, Luu et al. 1997, and Trujillo, Jewitt & Luu 2000). Although these classes are now well established, only rudimentary information has been collected about their populations. One reason is that only a fraction of the known KBOs were discovered in well-parametrized surveys that have been published in the open literature (principally Jewitt & Luu 1993 (1 KBO); Jewitt & Luu 1995 (17 KBOs); Irwin, Tremaine & Żytkow 1995 (2 KBOs); Jewitt, Luu & Chen 1996 (15 KBOs); Gladman et al. 1998 (5 KBOs); Jewitt, Luu & Trujillo 1998 (13 KBOs); and Chiang & Brown 1999 (2 KBOs)). In this work, we characterize the fundamental parameters of the CKBOs: the size distribution, inclination distribution and radial distribution using a large sample (86 KBOs) discovered in a well-characterized survey.

The quintessential measurement of the size distribution relies on the Cumulative Luminosity Function (CLF). The CLF describes the number of KBOs $`\text{deg}^2`$ ($`\mathrm{\Sigma }`$) near the ecliptic as a function of apparent red magnitude ($`m_R`$). It is fitted with the relation $`\mathrm{log}(\mathrm{\Sigma })=\alpha (m_Rm_0)`$, where $`m_0`$ is the red magnitude at which $`\mathrm{\Sigma }=1`$ KBO $`\text{deg}^2`$. The slope ($`\alpha `$) is related to the size distribution (described later). Although many different works have considered the CLF, two papers are responsible for discovering the majority of KBOs found in published surveys: Jewitt, Luu & Chen (1996) and Jewitt, Luu & Trujillo (1998). The former constrained the CLF over a 1.6 magnitude range ($`23.2<m_R<24.8`$) with 15 discovered KBOs while the latter covered a complimentary 2.5 magnitude range ($`20.5<m_R<23.0`$), discovering 13 objects. Jewitt, Luu & Trujillo (1998) measured the CLF produced from these two data sets and found $`\alpha =0.58\pm 0.05`$ and $`m_0=23.27\pm 0.11`$. Gladman et al. (1998) criticized this work on 2 main counts: (1) they believed that Jewitt, Luu & Chen (1996) underestimated the number of KBOs and (2) the fit in the Jewitt, Luu & Trujillo (1998) survey used a least-squares approach that assumed Gaussian errors rather than Poissonian errors. Gladman et al. (1998) found 5 additional KBOs and re-analysed the CLF using a Poissonian maximum likelihood method to refit the CLF to (1) the Jewitt, Luu & Trujillo (1998) data without the Jewitt, Luu & Chen (1996) data and (2) a fit to the 6 different surveys available at the time except for Tombaugh (1961), Kowal (1989) and Jewitt, Luu & Chen (1996). Both of these fits were steeper but formally consistent with the original Jewitt, Luu & Trujillo (1998) data at the $`1.5\sigma `$ level: (1) $`\alpha =0.72_{0.26}^{+0.30}`$ and $`m_0=23.3_{0.4}^{+0.2}`$ and (2) $`\alpha =0.76_{0.11}^{+0.10}`$ and $`m_0=23.40_{0.18}^{+0.20}`$. Chiang & Brown (1999) find a flatter size distribution of $`\alpha =0.52\pm 0.05`$ and $`m_0=23.5\pm 0.06`$ much closer to the Jewitt, Luu & Trujillo (1998) result. They observed that the steep size distribution reported by Gladman et al. (1998) was an artifact of their selective exclusion of part of the available survey data, not of their use of a different fitting method. The first goal of this work is to measure the CLF and additionally constrain the power-law slope of the size distribution using a single well-characterized survey and a maximum likelihood simulation which allows for the correction of observational biases.

An accurate characterization of the inclination distribution of the KBOs is critical to understanding the dynamical history of the outer Solar System since the era of planetesimal formation. We expect that the KBOs formed by accretion in a very thin disk of particles with a small internal velocity dispersion (e.g. Kenyon & Luu 1998 and Hahn & Malhotra 1999) and a correspondingly small inclination distribution. However, the velocity dispersion indicated by the inclination distribution in the present-day Kuiper Belt is large. Jewitt, Luu & Chen (1996) measured the apparent Half-Width of the Kuiper Belt inclination distribution to be $`5^{}`$. They noted a strong bias against observing high inclination objects in ecliptic surveys, and they estimated the true distribution to be much thicker, with an inclination distribution Half-Width of $`15^{}`$, corresponding to a vertical velocity dispersion of $`1`$ km/s. Several conjectures have been advanced to explain the thickness of the Kuiper Belt: Earth-mass planetesimals may have been scattered through the belt in the late stages of the planet-formation era, exciting the Kuiper Belt (Morbidelli & Valsecchi 1997 and Petit et al. 1999); stellar encounters may have enhanced the velocity dispersion of the distant KBOs (Ida, Larwood & Burkert, 2000); and the dispersion velocity of small bodies tends to grow to roughly equal the escape speed of the bodies contributing the most mass (the large bodies for size distributions with $`q<4`$) in the belt (Aarseth, Lin & Palmer 1993). As there is much speculation about the origin of the large velocity dispersion of the Kuiper Belt, but only one published measurement (Jewitt, Luu & Chen 1996), the second goal of this work is to accurately quantify the inclination distribution from our large sample of objects.

The radial extent of the Classical Kuiper Belt has not been well constrained. None of the CKBOs have been discovered beyond $`R50`$ AU. This trend was first noted by Dones (1997) who suggested that the 50–75 AU region may be depleted; he found the results of a Monte-Carlo simulation of CKBOs drawn from a rather flat differential size distribution (power-law index $`q=3`$) to be inconsistent with the observations of the 6 CKBOs discovered by Jewitt, Luu & Chen (1996). Jewitt, Luu & Trujillo (1998) discovered all of their KBOs at heliocentric distances $`R<46`$ AU. In the absence of other effects, one should expect to find fewer bodies with $`R>50`$ AU than with $`R40`$ AU, as the former are about a magnitude fainter than the latter. However, through the use of a Monte-Carlo model they demonstrated that the bias against objects beyond 50 AU is not strong enough to explain the distribution of discovery distances. They speculated that the lack of bodies discovered beyond 50 AU could be caused by a combination of (1) a decrease in the maximum KBO size (and reduction in the brightest and most detectable objects) beyond 50 AU or (2) the size distribution might steepen beyond 50 AU, putting more of the mass in the smaller, less-detectable bodies. They also suggested that the lack of $`R>50`$ AU objects could be explained by an outer edge to the Classical Kuiper Belt at 50 AU.

Two later papers questioned the existence of an edge to the Kuiper Belt near 50 AU. Gladman et al. (1998) suggested that the number of objects expected to be discovered beyond 50 AU is highly dependent on the size distribution because steep size distributions reduce the number of large (bright) bodies relative to small (faint) bodies. Gladman et al. (1998) adopted a relatively steep distribution ($`q=4.65`$), and found no significant evidence of a truncated belt. Chiang & Brown (1999) found that 8%–13% of the $`100`$ objects known at the time should have been found beyond 50 AU, and suggested that this precludes the presence of a density enhancement beyond 50 AU, but could not definitively rule out a density deficit. Allen, Bernstein & Malhotra (2001) have also recently reported the detection of an outer edge to the Kuiper Belt. The third goal of the present work is to test the distribution of the discovery distances for the presence of an outer edge to the Kuiper Belt.

## 2 Survey Data

Observations were made at the 3.6 m diameter Canada-France-Hawaii Telescope using the 12288 x 8192, 15 µm pixel mosaic CCD (CFHT 12k; Cuillandre et al. 2000). Built at the University of Hawaii (UH) the CFHT 12k comprises 12 edge-abutted, thinned, high quantum efficiency (QE $`0.75`$), 4096 x 2048 pixel Lincoln Laboratory CCDs. It is currently the largest close-packed CCD camera in the world. When mounted at the CFHT f/4 prime focus the camera yields a plate scale of 0.206 arc sec/pixel, corresponding to a 0.330 sq deg field of view in each 200 Mbyte image. Images were taken through a Mould $`R`$ filter, with a central wavelength of 6581 Å and a bandwidth of 1251 Å. Instrumental parameters of the survey are summarized in Table 1.

Observations were taken within a few days of new moon under photometric conditions during three periods: Feb 10 – 15 1999, Sep 5–8 1999, and Mar 31 – Apr 3 2000. Fields were imaged at airmasses $`<1.7`$ and were within 1.5 hours of opposition. We chose to use short 180 s exposures at the CFHT to maximize area coverage and detection statistics. All discovered objects were accessible for recovery at the UH 2.2 m telescope during comparable seeing conditions with exposure times of $`<600`$ seconds. Each field was imaged three times (a “field triplet”), with about 1 hour timebase between exposures. Fields imaged appear in Figures 1, 2 and 3, and in Table 3. The CFHT observations were taken at three ecliptic latitudes $`\beta =0^{}`$, $`10^{}`$, and $`20^{}`$ to probe the inclination distribution of the KBOs (see §4).

Photometric calibrations were obtained from Landolt (1992) standard stars imaged several times on each chip. Three CFHT 12k chips of poor quality were replaced between the Feb 1999 and Sep 1999 runs. The positions of four other CFHT 12k chips within the focal plane array were changed to move the cosmetically superior chips towards the center of the camera. The photometric calibration accounts for these changes, as shown in Table 2, containing the measured photometric zero points of the chips. In addition, chip 6 was not used in Feb 1999 because of its extremely poor cosmetic quality. The area covered in the fields from Feb 1999 was corrected for this 8% reduction in field-of-view. The area imaged in Mar 2000 included some small field overlap (6%), resulting in a minor correction applied to the reported total area surveyed.

Each of the 12 CCDs in the CFHT 12k functions as an individual detector, with its own characteristic bias level, flat field, gain level, and orientation (at the $`1^{}`$ level). The bias level for each chip was estimated using the row-by-row median of the overscan region. Flat fields were constructed from a combination of (1) the median of normalized bias-subtracted twilight flat fields and (2) a median of bias-subtracted data frames, with a clipping algorithm used to remove excess counts due to bright stars. Fields were analysed by subtracting the overscan region, dividing by the composite flats and searching for moving objects using our Moving Object Detection Software (MODS, Trujillo & Jewitt 1998). We rejected bad pixels through the use of a bad pixel mask.

Artificial moving objects were added to the data to quantify the sensitivity of the moving object detection procedure (Trujillo & Jewitt 1998). The seeing during the survey typically varied from 0.7 arc sec to 1.1 arc sec (FWHM). Accordingly, we subdivided and analysed the data in 3 groups based on the seeing. Artificial moving objects were added to bias-subtracted twilight sky-flattened images, with profiles matched to the characteristic point-spread function for each image group. These images were then passed through the data analysis pipeline. The detection efficiency was found to be uniform with respect to sky-plane speed in the 1 – 10 arc sec/hr range. At opposition, the apparent speed in arc sec/hr, $`\dot{\theta }`$, of an object is dominated by the parallactic motion of the Earth, and follows

$$\dot{\theta }148\left(\frac{1R^{0.5}}{R1}\right),$$
(1)
where $`R`$ is heliocentric distance in AU (Luu & Jewitt 1988). From Equation 1, our speed limit criterion for the survey, $`1\text{ ″}\text{ hr}^1<\dot{\theta }<10\text{ ″}\text{ hr}^1`$, corresponds to opposition heliocentric distances $`10\text{ AU}R140\text{ AU}`$, with efficiency variations within this range due only to object brightness and seeing.

The magnitude-dependent efficiency function was fitted by

$$\epsilon =\frac{\epsilon _{\mathrm{max}}}{2}\left(\mathrm{tanh}\left(\frac{m_{R50}m_R}{\sigma }\right)+1\right),$$
(2)
where $`0<\epsilon <1`$ is the efficiency with which objects of red magnitude $`m_R`$ are detected, $`\epsilon _{\mathrm{max}}`$ is the maximum efficiency, $`m_{R50}`$ is the magnitude at which $`\epsilon =\epsilon _{\mathrm{max}}/2`$, and $`\sigma `$ magnitudes is the characteristic range over which the efficiency drops from $`\epsilon _{\mathrm{max}}`$ to zero. Table 4 shows the efficiency function derived for each seeing category, along with an average of the seeing cases, weighted by sky area imaged, applicable to the entire data set. The efficiency function is known to greater precision than the $`0.1`$ magnitude uncertainty on our discovery photometry. Changes to the efficiency function of $`<0.1`$ magnitudes produce no significant variation in our results for the size or inclination distributions.

The MODS software, running on two Ultra 10 computers, was fast enough to efficiently search for the KBOs in near real-time, so that newly detected objects could be quickly discovered and re-imaged. We imaged $`35`$ field triplets each night at the CFHT, corresponding to $`20`$ Gbytes of raw data collected per night, plus several more Gbytes for flat fields and standard stars. Eighty-six KBOs were found in the CFHT survey, 2 of which were serendipitous re-detections of known objects. The discovery conditions of the detected objects appear in Table 5. Photometry was performed using a 2.5 arc second diameter synthetic aperture for discovery data, resulting in median photometric error of 0.15 magnitudes and a maximum photometric error of 0.3 for the faintest objects. Our results are unaffected by this error; randomly introducing $`\pm 0.15`$ magnitude errors in our simulations (described later), and $`\pm 0.3`$ magnitude errors in the faintest objects produced no statistically significant change. Trailing loss was insignificant as the KBOs moved only 0.15 arc sec during our integration.

### 2.1 Recovery Observations and Orbits

Extensive efforts were made to recover all objects using the UH 2.2 m telescope. Attempts were made to recover the objects one week after discovery, then one, two and three months after discovery. Most of these attempts were successful, as demonstrated by the fact that 79 of the 86 CFHT objects were recovered. The loss of 7 objects is the result of unusually poor weather during the Mar–May 1999 recovery period. Only 6 of the 79 recovered objects have arc-lengths shorter than 30 days as of Dec 1, 2000. Second opposition observations have been acquired for 36 of the 78 KBOs found in Feb 1999 and Sep 1999.

Orbits derived from the discovery and recovery data appear in Table 6. The listed elements are those computed by Brian Marsden of the Minor Planet Center. We also benefited from orbital element calculations by David Tholen (Univ. of Hawaii). Both sources produced comparable orbital solutions to the astrometric data.

With only first opposition observations, the inclination and heliocentric distance at discovery can be well determined for nearly all KBOs, as depicted in Figure 4. We find that the semimajor axis and eccentricity determinations are less reliable but are usually good enough to classify the objects as either Classical, Resonance or Scattered KBOs, as depicted in Figure 5. We find that 6 out of 36 (17%) of the objects exhibit orbital changes large enough for their dynamical classification to change from the first opposition to the second opposition. Randomly rejecting 17% of our sample (to simulate misclassification) does not significantly change the results. In addition, rejection of all but the multi-opposition objects does not significantly change our results; as expected, the total number of KBOs estimated decreased by a factor $`2`$ and error bars increased by a factor $`\sqrt{2}`$ due to the sample size reduction. The eccentricity and semimajor axes of all objects with $`a<50`$ AU (this includes all Classical KBOs) appear in Figure 6.

In the next two sections, we use our observations to constrain three fundamental quantities of the Classical KBOs: (1) the size distribution index, $`q`$, (2) the half-width of the inclination distribution, $`i_{1/2}`$, and (3) the total number of CKBOs larger than 100 km in diameter, $`N_{\mathrm{CKBOs}}(D>100\text{ km})`$. The quantities $`i_{1/2}`$ and $`q`$ are uncorrelated, as the observable constraining $`i_{1/2}`$ is the inclination distribution and the observable constraining $`q`$ is the absolute magnitude distribution. However, $`N_{\mathrm{CKBOs}}(D>100\text{ km})`$ is a function of both $`q`$ and $`i_{1/2}`$ as a steeper size distribution or thicker inclination distribution will each allow more bodies to be present. In the maximum likelihood simulations that follow, the ideal case would be to constrain $`q`$, $`i_{1/2}`$ and $`N_{\mathrm{CKBOs}}(D>100\text{ km})`$ and estimate errors in one simulation, however, this is difficult computationally. Therefore, we find the best-fit values of the three parameters in a single simulation, but estimate the errors on the parameters in two simulations, one that estimates the $`q`$-$`N_{\mathrm{CKBOs}}(D>100\text{ km})`$ joint errors and one that estimates the $`i_{1/2}`$-$`N_{\mathrm{CKBOs}}(D>100\text{ km})`$ joint errors. We then combine the two simulation results in quadrature to determine the errors on $`N_{\mathrm{CKBOs}}(D>100\text{ km})`$.

## 3 The Size Distribution of the Classical KBOs

We estimate the size distribution of the KBOs from our data in two ways. The first is a simple estimate made directly from the distribution of ecliptic KBO apparent magnitudes (CLF). The second is a model which simulates the discovery characteristics of our survey through the use of a maximum likelihood model constrained by the absolute magnitude of the Classical KBOs.

### 3.1 Cumulative Luminosity Function

We model the CLF with a power-law relation, $`\mathrm{log}\mathrm{\Sigma }=\alpha (m_Rm_0)`$1). The KBOs are assumed to follow a differential power-law size distribution of the form $`n(r)drr^qdr`$, where $`n(r)dr`$ is the number of objects having radii between $`r`$ and $`r+dr`$, and $`q`$ is the index of the size distribution. Assuming albedo and heliocentric distance distributions that are independent of KBO size, the simple transformation between the slope of the CLF ($`\alpha `$) and the exponent of the size distribution ($`q`$) is given by (Irwin et al. 1995),

$$q=5\alpha +1.$$
(3)
Under these assumptions, the size distribution can be estimated directly from the CLF.

We estimated the CLF by multiplying the detection statistics from the observed distribution of object brightnesses by the inverse of the detection efficiency. We assumed Poisson detection statistics, with error bars indicating the interval over which the integrated Poisson probability distribution for the observed number of objects contains 68.27% of the total probability (identical to the errors derived by Kraft, Burrows & Nousek 1991). This is nearly equal to the Gaussian case for all data points resulting from more than a few detections. We have included all 74 KBOs discovered in our 37.2 sq deg of ecliptic fields in the estimate of the CLF. This includes the lost objects, as the CLF is simply a count of the number of bodies discovered at a given apparent magnitude. Our results appear in Figure 7, with other published KBO surveys. All observations were converted to $`R`$-band if necessary assuming $`VR=0.5`$ for KBOs (Luu & Jewitt 1996), and error bars were computed assuming Poisson detection statistics. The data point of Cochran et al. (1995) near $`m_R=28`$ was omitted because of major uncertainties about its reliability (Brown, Kulkarni & Liggett 1997, cf. Cochran et al. 1998). Early photographic plate surveys (Tombaugh 1961, portions of Luu & Jewitt 1988, and Kowal 1989) have unproven reliability at detecting faint slow-moving objects, and plate emulsion variations and defects make accurate photometric calibration difficult. The photographic plate survey data were not used in our analysis.

The CLF points are highly correlated with one another, resulting in a heavy weighting of the bright object data points. Thus, we fitted the Differential Luminosity Function (DLF) instead. We plot the DLF data points at the faint end of the bin, representing the modal value in that bin. Very small bin sizes were chosen (0.1 magnitudes) to negate binning effects incurred from averaging the detection efficiency (Equation 2) over a large magnitude range. For any non-zero CLF slope, $`\alpha `$, the DLF and CLF slopes are equal due to the exponential nature of the CLF. The DLF was modelled by evaluating the Poisson probability of detecting the observed DLF given a range of $`m_0`$ and $`\alpha `$, with the maximum probability corresponding to our best-fit values. Error bars were determined by finding the contours of constant joint probability for $`m_0`$ and $`\alpha `$ enclosing 68.27% of the total probability, a procedure similar to that used below for the maximum likelihood simulation. Computations from this procedure are summarized in Table 7. We find that the slope of the CLF is $`\alpha =0.64_{0.10}^{+0.11}`$ with $`m_0=23.23_{0.20}^{+0.15}`$, which corresponds to $`q=4.2\pm 0.5`$ from Equation 3. We also fitted the CLF by applying the maximum likelihood method described by Gladman et al. (1998) to our data, which yields statistically identical results to the binned DLF procedure: $`\alpha =0.63\pm 0.06`$ and $`m_0=23.04_{0.09}^{+0.08}`$, corresponding to $`q=4.2\pm 0.3`$. The maximum-likelihood method provides slightly better signal-to-noise, and is independent of binning effects. However, it does not provide a visualization of the data, as presented for the DLF fit. We adopt the maximum-likelihood procedure as our formal estimate of the CLF slope. Both methods estimating the size distribution are in statistical agreement with the more detailed analysis presented in the next section.

The best-fit $`\alpha =0.63`$ magnitude distribution was compared with the observed magnitude distribution using a Kolmogorov-Smirnov test (Press et al. 1992), producing a value of $`D=0.13`$. If the model and the data distributions were identical, a deviation greater than this would occur by chance 12% of the time. Thus, our linear model is not a perfect fit, but it is statistically acceptable.

### 3.2 Maximum Likelihood Simulation

We now present more detailed analysis of the size distribution. Since we model the detection statistics of an assumed population, we choose to model the 49 Classical KBOs discovered on the ecliptic as they are numerically dominant in the observations and their orbital parameters are more easily modelled than other KBO classes. Our selection criteria for CKBOs are perihelion $`q^{}>37\text{ AU}`$ and $`40.5\text{ AU}<a<46\text{ AU}`$. Given the size of an object and its orbital parameters, we can compute its position, velocity, and brightness, allowing a full “Monte-Carlo” style analysis of the bias effects of our data collection procedures. The apparent brightness was computed from:

$$m=m_{}2.5\mathrm{log}(p\mathrm{\Phi }(\alpha ^{})r^2)+2.5\mathrm{log}(2.25\times 10^{16}R^2\mathrm{\Delta }^2),$$
(4)
where $`\alpha ^{}`$ is the phase angle of the object, $`\mathrm{\Phi }(\alpha ^{})`$ is the Bowell et al. (1989) phase function, geometric albedo is given by $`p`$, $`r`$ is the object radius in kilometers, $`R`$ is the heliocentric distance, and $`\mathrm{\Delta }`$ is the geocentric distance, both in AU (Jewitt & Luu 1995). The apparent red magnitude of the Sun was taken to be $`m_{}=27.1`$. For this work, we assume $`p=0.04`$, consistent with a Centaur-like albedo (Jewitt & Luu 2000). We neglect phase effects (setting $`\mathrm{\Phi }(\alpha ^{})1`$) since the maximum phase angle of an object at $`R=40`$ AU within 1.5 hours of opposition is $`\alpha ^{}=0.55\mathrm{°}`$. This corresponds to $`\mathrm{\Phi }(\alpha ^{})=0.91`$, a change in brightness of only 0.09 magnitudes, which is less than other uncertainties in the data.

This apparent brightness is used in a biasing-correction procedure (Trujillo, Jewitt & Luu 2000 and Trujillo 2000), summarized here:

1. A model distribution of KBOs is assumed (described in Table 8).
2. KBOs are drawn randomly from the model distribution.
3. For each KBO, the apparent speed and ecliptic coordinates are computed from the equations of Sykes & Moynihan (1996, a sign error was found in Equation 2 of their text and corrected), and compared to the observed fields and speed criteria.
4. The apparent magnitude is computed from Equation 4.
5. The efficiency function (Equation 2) and our field area covered are used to determine if the simulated object would be “detected” in our survey.
6. A histogram of the detection statistics for the simulated objects is constructed, logarithmically binned by object size for the size distribution model and binned by inclination for the inclination-distribution model. Binning effects were negligible due to small bin choice.
7. Steps 1-6 are repeated until the number of detected simulated objects is at least a factor 10 greater than the number of observed objects in each histogram bin (typically requiring a sample of $`10^6<N<10^8`$ simulated objects, depending on the observed distribution).
8. The likelihood of producing the observed population from the model is estimated by assuming that Poisson detection statistics ($`P=\frac{\mu ^n}{n!}\mathrm{exp}\left(\mu \right)`$) apply to each histogram bin, where $`\mu `$ represents the expected number of simulated objects “discovered” given the number of objects simulated and $`n`$ represents the true number of KBOs observed. Thus, the observed size distribution, calculated from Equation 4, is used to constrain the $`q`$ model, and the observed inclination distribution is used to constrain the $`i`$ model (§ 4).

These steps are repeated for each set of model parameters in order to estimate the likelihood of producing the observations for a variety of models.

For the size distribution analysis, we take our best-fit model of the width of the inclination distribution (Half-Width $`i_{1/2}=20^{}`$, as estimated in the next section), and vary the size distribution index $`q`$, and the total number of objects $`N_{\mathrm{CKBOs}}(D>100\text{ km})`$. Model parameters are summarized in Table 8 and results appear in Figure 8. Our best-fit values are

$$\begin{array}{cccc}\hfill q& =& 4.0_{0.5}^{+0.6}\hfill & (1\sigma )\text{ and}\hfill \\ & =& 4.0_{2.1}^{+1.3}\hfill & (3\sigma ),\hfill \\ & \text{and}& & \\ \hfill N_{\mathrm{CKBOs}}(D>100\text{ km})& =& 3.8_{1.5}^{+2.0}\times 10^4\hfill & (1\sigma )\text{ and}\hfill \\ & =& 3.8_{2.7}^{+5.4}\times 10^4\hfill & (3\sigma ),\hfill \end{array}$$

where the errors for $`N_{\mathrm{CKBOs}}(D>100\text{ km})`$ have been combined in quadrature from the results of the $`q`$ and $`i_{1/2}`$ fits, as described at the end of § 2.1. The values for $`q`$ are consistent with previously published works (Table 9) and the $`q`$ derived from the CLF data in the simple model (Equation 3). The results are consistent with the distribution of large ($`D>150`$ km) main-belt asteroids ($`q=4.0`$, Cellino, Zappalá, & Farinella 1991) and rock crushed by hypervelocity impacts ($`q=3.4`$, Dohnanyi 1969). In addition, the scenario where the cross-sectional area (and thus optical scattered light and thermal emission) is concentrated in the largest objects ($`q<3`$, Dohnanyi 1969) is ruled out at the $`>2\sigma `$ ($`>95.4\%`$ confidence) level. Our results are also consistent with Kenyon & Luu (1999) who simulate the growth and velocity evolution of the Kuiper Belt during the formation era in the Solar System. They find several plausible models for the resulting size distribution, all of which have $`q4`$. In Figure 9 we plot the best-fit model CKBO distribution with the observed DLF to demonstrate the expected results from different size distributions.

The magnitude distribution expected from the maximum likelihood model was compared to the observed magnitude distribution, as was done for the CLF-derived magnitude distribution in § 3.1. The Kolmogorov-Smirnov test produced $`D=0.17`$; a greater deviation would occur by chance 11% of the time.

In our Classical KBO maximum likelihood simulation, we have ignored possible contributions of the 7 lost KBOs, since their orbital classes are not known. However, including them in the simulations by assuming circular orbits at the heliocentric distance of discovery results in statistically identical resulted for $`q`$, and the expected 7/49 rise in $`N_{\mathrm{CKBOs}}(D>100\text{ km})`$.

## 4 Inclination Distribution of the Classical KBOs

The dynamical excitation of the Kuiper Belt is directly related to the inclination distribution of the KBOs. We present the inclinations of the CKBOs found in the CFHT survey in Figure 10. Assuming heliocentric observations, a KBO in circular orbit follows

$$\mathrm{sin}\beta =\mathrm{sin}i\mathrm{sin}f$$
(5)
where $`\beta `$ is the heliocentric ecliptic latitude, $`0<i<90^{}`$ is the inclination, and $`0<f<360^{}`$ represents the true anomaly of the object’s orbit with $`f=0`$ and $`180\mathrm{°}`$ representing the ecliptic plane crossing (the longitude of perihelion is defined as 0 in this case). Using Equation 5, we plot the fraction of each orbit spent at various ecliptic latitudes as a function of $`i`$ (Figure 11). This plot demonstrates two trends concerning the ecliptic latitude of observations $`\beta _{\mathrm{obs}}`$. First, high-inclination objects are a factor 3–4 times more likely to be discovered when $`\beta _{\mathrm{obs}}i`$ than when observing at low ecliptic latitudes ($`\beta _{\mathrm{obs}}<i`$). Second, the number of expected high-inclination objects drops precipitously, roughly as $`1/i`$, once $`i>1.5\beta _{\mathrm{obs}}`$ (Jewitt, Luu & Chen 1996).

These facts led us to observe at three different ecliptic latitudes ($`0^{}`$, $`10^{}`$ and $`20^{}`$) to better sample the high-inclination objects. During two observation periods (Sep 1999 and Mar 2000) care was made to interleave the ecliptic fields with the off-ecliptic fields on timescales of $`30`$ minutes. This technique provides immunity to drift in the limiting magnitude which might otherwise occur in response to typical slow changes in the seeing through the night. The results for the robust, interleaved fields matched those for the seeing-corrected Feb 1999 fields where fields were interleaved on much longer timescales of $`3`$ hours. Accordingly, we combined the data sets from all epochs to improve signal-to-noise. In the next sections, we analyse the inclination distribution using two techniques to demonstrate the robustness of our method.

### 4.1 Simple Inclination Model

First, since fields were imaged at three different ecliptic latitudes, the surface density of objects at each latitude band ($`\mathrm{\Sigma }(0\mathrm{°})`$, $`\mathrm{\Sigma }(10\mathrm{°})`$ and $`\mathrm{\Sigma }(20\mathrm{°})`$) can directly yield the underlying inclination distribution. In our simple model, we generate an ensemble of inclined, circular orbits drawn from a Gaussian distribution centered on the ecliptic, and having a characteristic Half-Width of $`i_{1/2}`$. The probability of drawing a KBO with inclination between $`i`$ and $`i+di`$ is given by

$$P(i)di=\frac{1}{\sigma \sqrt{2\pi }}\mathrm{exp}\left(\frac{i^2}{2\sigma ^2}\right)di,$$
(6)
where $`\sigma =i_{1/2}/\sqrt{2\mathrm{ln}2}`$. Using this relation, and Equation 5, we simulate the expected values of $`\mathrm{\Sigma }(0\mathrm{°})`$, $`\mathrm{\Sigma }(10\mathrm{°})`$ and $`\mathrm{\Sigma }(20\mathrm{°})`$ for various $`i_{1/2}`$. These are compared to two ratios measured from our observations, $`R(10^{},0^{})\mathrm{\Sigma }(10\mathrm{°})/\mathrm{\Sigma }(0\mathrm{°})`$ and $`R(20^{},0^{})\mathrm{\Sigma }(20\mathrm{°})/\mathrm{\Sigma }(0\mathrm{°})`$. Results appear in Table 10, and demonstrate that the characteristic half-width of the inclination distribution in the Kuiper Belt is $`i_{1/2}17\mathrm{°}_{4\mathrm{°}}^{+10\mathrm{°}}`$ ($`1\sigma `$ = 68.27% confidence). This simple model does not use the observed inclination distribution of the individual objects, merely the surface density of objects found at each ecliptic latitude, thus we have combined all objects from all KBO classes into this estimate.

### 4.2 Full Maximum Likelihood Inclination Model

Second, we use the maximum likelihood model described in §3.1. We list the parameters of the model in Table 11. This model encompasses the additional constraint of the observed inclination distribution, as well as the parallactic motion of the Earth and KBO orbital motion to produce more realistic results. Results appear in Figure 12, with $`N_{\mathrm{CKBOs}}(D>100\text{ km})`$ representing the number of CKBOs with diameters greater than 100 km. The maximum likelihood occurs at

$$\begin{array}{cccc}\hfill i_{1/2}& =& 20\mathrm{°}_{4\mathrm{°}}^{+6\mathrm{°}}\hfill & (1\sigma )\text{ and}\hfill \\ & =& 20\mathrm{°}_{8\mathrm{°}}^{+26\mathrm{°}}\hfill & (3\sigma ),\hfill \\ & \text{and}& & \\ \hfill N_{\mathrm{CKBOs}}(D>100\text{ km})& =& 3.8_{1.5}^{+2.0}\times 10^4\hfill & (1\sigma )\text{ and}\hfill \\ & =& 3.8_{2.7}^{+5.4}\times 10^4\hfill & (3\sigma ),\hfill \end{array}$$

where the errors for $`N_{\mathrm{CKBOs}}(D>100\text{ km})`$ have been estimated from the $`i_{1/2}`$ and $`q`$ fits, combined in quadrature, as described at the end of § 2.1. This maximum likelihood model is consistent with the simple model described in §4.1. In Figure 13, we plot the observed surface density of objects as a function of ecliptic latitude and compare these data to our best-fit models. This illustrates the fundamental fact that even though the true inclination distribution of the KBOs is very thick ($`i_{1/2}20\mathrm{°}`$), the surface density drops off quickly with ecliptic latitude, reaching half the ecliptic value at an ecliptic latitude of $`\beta 3\mathrm{°}`$ ($`\mathrm{\Sigma }(3\mathrm{°})/\mathrm{\Sigma }(0\mathrm{°})<0.5`$).

The functional form of the inclination distribution cannot be well constrained by our data. However, the best-fit Gaussian distribution was compared to a flat-top (“top-hat”) inclination distribution, with a uniform number of objects in the $`0\mathrm{°}<i<30\mathrm{°}`$ range. The Gaussian and flat-top models were equally likely to produce the observed distribution in the 65% confidence limit ($`<1\sigma `$). A Gaussian model multiplied by $`\mathrm{sin}(i)`$ was also tried but could be rejected at the $`>3\sigma `$ level because it produced too few low-inclination objects. We also tested the best-fit model presented by Brown (2001), consisting of two Gaussians multiplied by $`\mathrm{sin}(i)`$,

$$[a\mathrm{exp}(\frac{i^2}{2\sigma _1^2})+(1a)\mathrm{exp}(\frac{i^2}{2\sigma _2^2})]\mathrm{sin}i,$$
(7)
where $`a=0.93`$, $`\sigma _1=2.2\mathrm{°}`$, and $`\sigma _2=18\mathrm{°}`$, and found it equally compatible with our single Gaussian model (Equation 6). Because the Gaussian model was the simplest model that fit the observed data well, we chose it to derive the following velocity dispersion results.

We first find the mean velocity vector of all the simulated best-fit CKBOs, $`\stackrel{}{\overline{v}}`$, in cylindrical coordinates (normal vectors $`\widehat{r}`$, $`\widehat{\theta }`$, and $`\widehat{z}`$ representing the radial, longitudinal and vertical components respectively). The mean velocity vector $`\stackrel{}{\overline{v}}`$ is consistent with a simple Keplerian rotation model at $`R46`$ AU. We then compute the relative velocity of each KBO from this via $`|\stackrel{}{\overline{v}}\stackrel{}{v_i}|`$, where $`\stackrel{}{v_i}`$ is the velocity dispersion contribution of the $`i`$th KBO. We find the resulting root-mean-square (RMS) velocity dispersion of the $`\widehat{r}`$, $`\widehat{\theta }`$, and $`\widehat{z}`$ components to be equal to $`\mathrm{\Delta }v_r=0.51`$ km/s, $`\mathrm{\Delta }v_\theta =0.50`$ km/s, and $`\mathrm{\Delta }v_z=0.91`$ km/s, combining in quadrature for a total velocity dispersion of $`\mathrm{\Delta }v=\sqrt{\mathrm{\Delta }v_r^2+\mathrm{\Delta }v_\theta ^2+\mathrm{\Delta }v_z^2}=1.16`$ km/s. An error estimate of the velocity dispersion can be found by following a similar procedure for the $`i_{1/2}`$ = 16°and 26°($`\pm 1\sigma `$) models, yielding $`\mathrm{\Delta }v=1.16_{0.16}^{+0.25}`$ km/s.

### 4.3 Inferred Mass

The Kuiper Belt mass inferred from these results can be directly calculated from the size distribution and the number of bodies present. For the best-fit $`q=4.0`$ size distribution, the mass of CKBOs in bodies with diameters $`D_{\mathrm{min}}<D<D_{\mathrm{max}}`$ is

$$M(D_{\mathrm{min}},D_{\mathrm{max}})=\frac{4}{3}\pi \rho \mathrm{\Gamma }\mathrm{ln}(D_{\mathrm{max}}/D_{\mathrm{min}}),$$
(8)
where $`\rho `$ is the bulk density of the object. The normalization constant $`\mathrm{\Gamma }`$ is calculated from the results of our simulation,

$$\mathrm{\Gamma }3.0\times 10^{12}\text{ m}^3p_R^{1.5}N(D>100\text{ km}),$$
(9)
where $`N(D>100\text{ km})=3.8\times 10^4`$4.2), yielding $`\mathrm{\Gamma }=1.4\times 10^{19}\text{ m}^3`$ assuming $`p_R0.04`$. The mass for $`100\text{ km}<D<2000\text{ km}`$ then becomes

$$M(100\text{ km},2000\text{ km})0.03M_{\mathrm{}}\left(\frac{\rho }{1000\text{ kg}\text{ m}^3}\right)\left(\frac{0.04}{p_R}\right)^{1.5},$$
(10)
where $`M_{\mathrm{}}=6.0\times 10^{24}`$ kg is the mass of the earth. The uncertainties on this value are considerable as the characteristic albedo and density of the CKBOs are unknown.

### 4.4 Comparison of the Classical KBOs to Other Dynamical Classes

We found that the total number of CKBOs is given by $`N_{\mathrm{CKBOs}}(D>100\text{ km})=3.8_{1.5}^{+2.0}\times 10^4`$. This can be compared to the other main dynamical populations (the Resonant and Scattered KBOs) from our data. Observational biases favor the detection of the Plutinos over the Classical KBOs due to their closer perihelion distance. We found only 7 Plutinos (4 ecliptic and 3 off-ecliptic) so we can make only crude (factor $`2`$ statements) about the true size of the population. Thus, we use the results of Jewitt, Luu & Trujillo (1998) who estimate that the apparent fraction of Plutinos ($`P_a`$) in the Kuiper Belt is enhanced relative to the intrinsic fraction ($`P_i`$) by a factor $`P_a/P_i2`$ for $`q=4.0`$ and $`r_{\mathrm{max}}=1000`$ km. Applying this correction to our ecliptic observations (4 Plutinos and 49 Classical KBOs) indicates that the total number of Plutinos larger than 100 km in diameter is quite small,

$$N_{\mathrm{Plutinos}}(D>100\text{ km})\frac{4}{4+49}\frac{P_i}{P_a}N_{\mathrm{CKBOs}}1400.$$
(11)
The populations of the Plutinos and the 2:1 Resonant objects are important measures of the resonance sweeping hypothesis (Malhotra 1995), which predicts equal numbers of objects in each resonance. Since the 2:1 objects are systematically farther from the sun than the Plutinos, the true Plutino/2:1 ratio is higher than the observed ratio. Jewitt, Luu & Trujillo (2000) estimate the observed/true bias correction factor to be $`0.310`$ for a survey similar to ours ($`q=4`$ and $`m_{R50}=24.0`$). Only 2 of our objects (both found on the ecliptic) are $`<0.5`$ AU from the 2:1 Resonance, so we find the Plutino/2:1 fraction is given by $`(4/2)0.3100.6`$. Due to the small number of bodies involved, this is only an order of magnitude estimate. Within the uncertainties, our observations are consistent with the hypothesis that the 3:2 and 2:1 resonances are equally populated.

The observational biases against the Scattered KBOs are considerable. Trujillo, Jewitt & Luu (2000) estimate the total population of the Scattered KBOs to be $`N_{\mathrm{SKBOs}}(D>100\text{ km})=3.1_{1.3}^{+1.9}\times 10^4`$, approximately equal to the population of Classical KBOs derived from our data. We summarize the relative populations by presenting their number ratios:

$$\text{Classical}:\text{Scattered}:\text{Plutino}:\text{Resonant 2:1}=1.0:0.8:0.04:0.07.$$
(12)
## 5 The Edge of the Classical Kuiper Belt

We found no objects beyond heliocentric distance $`R_{\mathrm{obs}}=48.9`$ AU. There are two possibilities to explain this observation: (1) this is an observational bias effect and the bodies beyond $`R_{\mathrm{obs}}`$ cannot be detected in our survey, or (2) there is a real change in the physical or dynamical properties of the KBOs beyond $`R_{\mathrm{obs}}`$. In order to test these two possibilities, we compare the expected discovery distance of an untruncated Classical Kuiper Belt to the observations, as depicted in Figure 14. This untruncated CKBO distribution is identical to our best-fit model from § 4.2, except $`40.5\text{ AU}<a<200\text{ AU}`$, instead of $`40.5\text{ AU}<a<46\text{ AU}`$. The total number of bodies produced was considered a free parameter in this model. Inspecting Figure 14, the absence of detections beyond 50 AU is inconsistent with an untruncated model with $`R^2`$ radial power to the ecliptic plane surface density. Assuming Poisson statistics apply to our null detection beyond $`R_{\mathrm{max}}`$, the 99.73% ($`3\sigma `$) upper limit to the number of bodies ($`\mu `$) expected beyond $`R_{\mathrm{max}}`$ can be calculated from $`10.9973=\mathrm{exp}(\mu )`$, yielding $`\mu =5.9`$ KBOs. We found 49 ecliptic Classical KBOs inside the $`R_{\mathrm{max}}`$ limit, so the $`3\sigma `$ upper limit to the number density of KBOs beyond $`R_{\mathrm{max}}`$ is $`49/5.98`$ times less than the number density of Classical KBOs. Although we have constrained the outer edge by the heliocentric distance at discovery $`R`$, which is a directly observable quantity, a dynamical edge would be set by the semimajor axes ($`a`$) of the object orbits. This difference has little effect on our findings as the known CKBOs occupy nearly circular orbits with median eccentricity $`e=0.08`$ (the calculated median is conservative as it includes only bodies with $`e>0`$ to protect against short-arc orbits, which typically assume $`e=0`$). Since an untruncated distribution (1) is incompatible with our data, we must conclude that scenario (2) applies — there must be a physical or dynamical change in the KBOs beyond $`R_{\mathrm{max}}`$.

There are several possible physical and dynamical mechanisms that could produce the observed truncation of the belt beyond $`R_{\mathrm{max}}=50`$ AU (Jewitt, Luu & Trujillo 1998): (1) the size distribution of the belt might become much steeper beyond $`R_{\mathrm{max}}`$, putting most of the mass of the belt in the smallest, undetectable objects; (2) the size distribution could be the same ($`q=4`$), but there might be a dearth of large (i.e. bright) objects beyond $`R_{\mathrm{max}}`$, suggesting prematurely arrested growth; (3) the objects beyond $`R_{\mathrm{max}}`$ may be much darker and therefore remain undetected; (4) the eccentricity distribution could be lower in the outer belt, resulting in the detection of fewer bodies; (5) the ecliptic plane surface density variation with radial distance may be steeper than our assumed $`p=2`$; and (6) there is a real drop in the number density of objects beyond $`R_{\mathrm{max}}`$. We consider each of these cases in turn, and their possible causes.

Detailed simulations of the growth of planetesimals in the outer Solar System have not estimated the radial dependence of the formation timescale (e.g. Kenyon & Luu 1999). However, it is expected that growth timescales should increase rapidly with heliocentric distance, perhaps as $`tR^3`$ (Wetherill 1989). One could then expect a reduction in the number of large objects beyond 50 AU, as per (2) above, and a correspondingly steeper size distribution, as in (1), at larger heliocentric distances. However, with $`tR^3`$, the timescales for growth at $`R=41`$ AU (inner edge) and $`R=50`$ AU (outer edge) are only in the ratio 1.8:1. In addition, we observe no correlation between size and semimajor axis among the Classical KBOs.

To test scenario (1), we took our untruncated best-fit model and varied the size distribution index $`q_{\mathrm{out}}`$ for bodies with semimajor axes $`a>R_{\mathrm{max}}`$, keeping the KBO mass across the $`R_{\mathrm{max}}`$ boundary constant. We then found the minimum $`q_{\mathrm{out}}`$ consistent with our null detection beyond $`R_{\mathrm{max}}`$. This mass-conservation model is very sensitive to the chosen minimum body radius $`r_{\mathrm{min}}`$, because for any $`q_{\mathrm{out}}>4`$, most of the mass is in the smallest bodies (Dohnanyi 1969). The minimum size-distribution index required as a function of $`r_{\mathrm{min}}`$ appears in Table 12. If mass is conserved for the observable range of bodies, $`r_{\mathrm{min}}=50`$ km, the observed edge cannot be explained by a change in the size distribution unless $`q>10`$ ($`3\sigma `$), an unphysically large value. For the conservative case of $`r_{\mathrm{min}}=6`$ km (roughly the size of cometary nuclei, Jewitt 1997), the observed edge could only be explained by $`q>5.6`$ ($`3\sigma `$) beyond $`R_{\mathrm{max}}`$. We know of no population of bodies with a comparably steep size distribution. Thus, we conclude that the observed edge is unlikely to be solely caused by a change in the size distribution beyond $`R_{\mathrm{max}}`$.

A similar procedure was followed for possibility (2). Here again, we took our best-fit truncated model and extended it to large heliocentric distances. Then, $`r_{\mathrm{max}}`$ was varied to find the largest value that could explain our null detection beyond $`R_{\mathrm{max}}`$, keeping the total number density of objects with radii $`r<r_{\mathrm{max}}`$ constant. We found that $`r_{\mathrm{max}}<75`$ km ($`3\sigma `$) was required beyond $`R_{\mathrm{max}}`$ to explain the observed edge. This is a factor $`5`$ smaller radius and a factor $`150`$ less volume than our largest object found within $`R_{\mathrm{max}}`$ (1999 $`\mathrm{CD}_{158}`$, $`400`$ km in radius). Such a severe change in the maximum object size beyond $`R_{\mathrm{obs}}`$ would have to occur despite the fact that growth timescales vary by less than a factor of $`2`$ over the observed Classical KBO range, as explained above.

One might also expect (3) to be true, as KBO surfaces could darken over time with occasional resurfacing by collisions (Luu & Jewitt 1996), and long growth timescales indicate long collision timescales as well. However, the geometric albedo would have to be $`p<0.008`$, a factor 5 lower than that of the CKBOs in our model, assuming a constant number density of objects across the transition region. We are not aware of natural planetary materials with such low albedos.

The dynamical cases (4), a drop in the eccentricity distribution, and (5), a steeper ecliptic plane density index, can also be rejected. Even an extreme change in the eccentricity distribution cannot explain our observations. Lowering eccentricity from $`e=0.15`$ (a high value for the Classical KBOs) to $`e=0`$ results in a perihelion change from 42.5 AU to 50 AU for an object with semimajor axis 50 AU. Such a change corresponds to a 0.7 magnitude change in perihelion brightness, and to a factor 2.8 change in the surface density of objects expected from our $`\alpha =0.63`$ CLF. This model is rejected by our observations at the $`>5\sigma `$ level. The variation in ecliptic plane surface density with respect to heliocentric distance was assumed to follow a power law with index $`p=2`$ in our model. However, even a large increase to $`p=5`$ would result in a reduction in surface density of a factor 2.7 in the 41 AU to 50 AU range, which can also be rejected as the cause of our observed edge at the $`>5\sigma `$ level.

Since scenarios (1) through (5) seem implausible at best, we conclude that the most probable explanation for the lack of objects discovered beyond $`R_{\mathrm{max}}`$ is (6), the existence of a real, physical decrease in object number density. There have been few works considering mechanisms for such truncation. The 2:1 mean-motion Neptune resonance at $`a47.8`$ AU is quite close to the observed outer edge of the belt. However, given the Neptune resonance sweeping model (Malhotra 1995), the resonance could not cause an edge. The sweeping theory predicts that the 2:1 resonance should have passed through the Classical Kuiper Belt as Neptune’s orbit migrated outwards to its present semimajor axis. Thus, the KBOs interior to the current 2:1 resonance ($`a47.8`$ AU) could have been affected by this process, but an edge at $`R_{\mathrm{max}}`$ cannot be explained by such a model. Ida, Larwood & Burkert (2000) simulate the effect of a close stellar encounter on the Kuiper Belt, suggesting that KBO orbits beyond 0.25–0.3 times the stellar perihelion distance would be disrupted and ejected for a variety of encounter inclinations. Thus, an encounter with a solar mass star with perihelion at $`200`$ AU might explain the observed edge. Such encounters are implausible in the present solar environment but might have been more common if the sun formed with other stars in a dense cluster.

## 6 Constraints on a Distant Primordial Kuiper Belt

While our observations indicate a dearth of objects beyond 50 AU, it is also possible that a “wall” of enhanced number density exists at some large $`R100`$ AU distance, as suggested by Stern (1995). We know that the Kuiper Belt has lost much mass since formation because the present mass is too small to allow the observed objects to grow in the age of the solar system. Kenyon & Luu (1999) found that the primordial Kuiper Belt mass in the $`30\text{ AU}<R<50\text{ AU}`$ region could have been some $`10M_{}`$, compared to the $`0.1M_{}`$ we see today. Stern (1995) also speculated that the primordial surface density may be present at large heliocentric distances. We model this primordial belt as analogous to the CKBOs in terms of eccentricity, inclination and size distribution, but containing a factor of 100 more objects and mass per unit volume of space. These objects would be readily distinguishable from the rest of the objects in our sample as they would have low eccentricities characteristic of the CKBOs ($`e<0.25`$) yet would have very large semimajor axes ($`a>90`$ AU). Since we have discovered no such “primordial” objects, Poisson statistics state that the $`3\sigma `$ upper limit to the sky-area number density of primordial KBOs is 5.9 in 37.2 sq deg, or 0.16 primordial KBOs $`\text{deg}^2`$. We constrain the primordial KBOs by allowing the inner edge of the population, $`a_{\mathrm{min}}`$, to vary outwards, while keeping the outer edge fixed at 250 AU. We find that $`a_{\mathrm{min}}=130`$ AU coincides with the $`3\sigma `$ limit on the inner edge of the belt, nearly at the extreme distance limit of our survey. An object discovered at our survey magnitude limit $`m_{R50}=23.7`$ at this distance would have diameter $`D1800`$ km (approximately 25% smaller than Pluto) assuming a 4% albedo.

## 7 Summary

New measurements of the Kuiper Belt using the world’s largest CCD mosaic array provide the following results in the context of our Classical KBO model.

(1) The slope of the differential size distribution, assumed to be a power law, is $`q=4.0_{0.5}^{+0.6}`$ ($`1\sigma `$). This is consistent with accretion models of the Kuiper Belt (Kenyon & Luu 1999). This distribution implies that the surface area, the corresponding optical reflected light and thermal emission are dominated by the smallest bodies.

(2) The Classical KBOs inhabit a thick disk with Half-Width $`20\mathrm{°}_{4\mathrm{°}}^{+6\mathrm{°}}`$ ($`1\sigma `$).

(3) The Classical KBOs have a velocity dispersion of $`1.16_{0.16}^{+0.25}`$ km/s.

(4) The population of Classical KBOs larger than 100 km in diameter $`N_{\mathrm{CKBOs}}(D>100\text{ km})=3.8_{1.5}^{+2.0}\times 10^4`$ ($`1\sigma `$). The corresponding total mass of bodies with diameters between 100 km and 2000 km is $`M(100\text{ km},2000\text{ km})0.03M_{\mathrm{}}`$ , assuming geometric red albedo $`p_R0.04`$ and bulk density $`\rho 1000\text{ kg}\text{ m}^3`$.

(5) The approximate population ratios of the Classical, Scattered, 3:2 Resonant (Plutinos) and 2:1 Resonant KBOs are 1.0:0.8:0.04:0.07.

(6) The Classical Kuiper Belt has an outer edge at $`R=50`$ AU. This edge is unlikely to be due to a change in the physical properties of the CKBOs (albedo, maximum object size, or size distribution). The edge is more likely a real, physical depletion in the number of bodies beyond $`50`$ AU.

(7) There is no evidence of a primordial (factor 100 density increase) Kuiper Belt out to heliocentric distance $`R=130`$ AU.

We thank Dave Tholen and Brian Marsden for providing orbits and ephemerides. We appreciate the vital observational assistance provided by Scott Sheppard and his help with astrometric measurements. We thank David Woodworth, Ken Barton, Lisa Wells and Christian Viellet for help at the CFH telescope. We are grateful for the assistance of John Dvorak, Chris Merrick, Lance Amano, Paul DeGrood and Farren Herron-Thorpe at the University of Hawaii 2.2 m telescope. A NASA grant to DCJ provided financial support for this project.