COAGULATION CALCULATIONS OF ICY PLANET FORMATION AT 15–150 AU: A CORRELATION BETWEEN THE MAXIMUM RADIUS AND THE SLOPE OF THE SIZE DISTRIBUTION FOR TRANS-NEPTUNIAN OBJECTS

Figure 4 illustrates the time evolution of r_max for the calculation described in Figure 2. In each annulus, the largest objects first grow slowly; r_max grows roughly linearly with time. Once runaway growth begins, large objects grow exponentially; r_max reaches ∼1000 km in 3 Myr at 15–18 AU, 20 Myr at 22–27 AU, 50 Myr at 33–40 AU, and 200 Myr at 48–58 AU. After planets achieve these radii, the collisional cascade removes small planetesimals more rapidly than planets accrete them. Thus, planets grow much more slowly. Despite the different timescales, planets in all annuli reach the same maximum size of roughly 3000 km.

Zoom In Zoom Out Reset image size

Figure 5 shows how q₄ evolves with time. When the first objects with r > 100 km form, there are many 100 km objects and only a few with larger sizes. Thus, the size distribution at large sizes is very steep and q₄ is large. As the evolution proceeds, runaway growth produces more and more large objects from the ensemble of 100 km objects. With more objects at r > 500 km and fewer at 100 km, the power-law size distribution grows shallower and shallower (Figure 2); q₄ decreases. In the late stages of oligarchic growth, planets reach their maximum radius. Because large objects mostly accrete material from much smaller objects, the slope then converges on a characteristic value q₄ ≈ 3 (Goldreich et al. 2004; Schlichting & Sari 2011). Although the growth of the largest objects is monotonic in time (Figure 4), small number statistics produces the fluctuations around the steady decrease of q₄ with time shown in the figure.

Zoom In Zoom Out Reset image size

To examine the relationship between r_max and q₄, we consider calculations for each set of initial conditions (x_m, r₀, q_init, f_i). For each set of calculations, we collect all pairs (r_max, q₄) in all annuli at all times. We then bin these data in increments of 0.1 in log r_max and calculate the median q₄ in each r_max bin. This procedure yields the variation of median q₄ as a function of r_max and the initial conditions. This median provides a good representation of the typical slope: the inter-quartile range is roughly 2 when r_max is small (∼100 km) and roughly 0.5–1.0 when r_max is large (≳1000 km).

Figure 6 shows results for the complete set of calculations at a = 15–75 AU with r₀ = 1 km and the f_w fragmentation parameters. For the ensemble of massive disks with x_m = 3 (violet points in the figure), the median slope is large, q₄ ≈ 10, when r_max ≈ 100–300 km. For r_max ≈ 200–2000 km, the median q₄ declines roughly monotonically with increasing r_max. As r_max grows from ∼2000 km to ∼5000 km, the median q₄ is roughly constant.

Zoom In Zoom Out Reset image size

The variation of q₄ with r_max is sensitive to the initial disk mass. As the initial disk mass declines from x_m = 3 to x_m = 0.01, growth times are longer and longer. Due to the impact of gas drag and long-range stirring, small planetesimals in lower mass disks tend to have larger eccentricities at the onset of runaway growth than those in higher mass disks. Larger eccentricities reduce gravitational focusing factors, slowing growth and enabling stirring to grow faster. Because slower growth is more orderly, large objects have shallower size distributions. When r_max ≈ 100–300 km, lower mass disks thus have a smaller median q₄ than higher mass disks. Because growth is very slow for the smallest disk masses, all calculations with x_m ≲ 0.1 have roughly the same median q₄ ≈ 4.

For disks with x_m ≳ 0.3, gas drag (long-range stirring) is more (less) effective. With larger gravitational focusing factors in more massive disks, growth is relatively fast. Once a few 100 km objects form, a stronger runaway leads to a more pronounced decline of q₄ with r_max. As r_max increases from roughly 100 km to ∼2000 km, the median q₄ declines roughly linearly with r_max. Once r_max ≳ 2000 km, the median q₄ is roughly constant. In lower mass disks, the median q₄ declines more slowly with increasing r_max. Because low-mass disks never produce objects with r_max ≳ 2000 km, the median q₄ is never constant with increasing r_max.

In addition to the clear dependence on x_m, the median q₄ is also sensitive to the initial maximum planetesimal size r₀ (Figure 7). For r_max ≈ 100–300 km and x_m ≈ 1–3, the median q₄ steadily increases with larger r₀, from q₄ ≈ 8–9 for r₀ = 1 km to q₄ ≈ 10–11 for r₀ = 10 km to q₄ ≈ 13 for r₀ = 100 km. Independent of r₀, the median q₄ declines steadily with larger r_max. At the largest r_max attained in our calculations, the median q₄ converges on an equilibrium slope ∼ 3–4. Calculations with smaller r₀ yield smaller values for this equilibrium slope.

Zoom In Zoom Out Reset image size

The behavior of the q₄–r_max relation is remarkably independent of semimajor axis (Table 2, Figure 8). In calculations at 30–150 AU, our 10 Gyr evolution times yield larger r_max than the 1 Gyr evolution times for calculations at 15–75 AU. Otherwise, the median slope at small r_max, the monotonic decline in slope at r_max ≈ 300–2000 km, and the roughly constant slope for r_max ≳ 2000 km closely follow the relations established at 15–75 AU. Our results show some indication for a smaller inter-quartile range in the median q₄ at larger semimajor axis. Confirming this trend requires a larger ensemble of calculations.

Zoom In Zoom Out Reset image size

Despite the insensitivity to semimajor axis, the median slope depends strongly on the initial size distribution (Table 2, Figure 9). In calculations where the initial mass is distributed among planetesimals with sizes of 1 m to 1 km, the median q₄ is systematically smaller (∼1–2) throughout the evolution. This behavior is independent of r₀. All of these calculations begin with a substantial amount of mass in objects with sizes of 0.1 km or smaller. During runaway growth, collisional grinding removes some of this material from the grid. Compared to calculations where all of the initial mass is concentrated in the largest planetesimals, the largest objects in these calculations accrete material from a lower mass reservoir of small planetesimals which are easier to erode. Thus, the growth of the large objects is relatively slower and the median slope is smaller. Because calculations with larger r₀ have a smaller fraction of the initial mass in small planetesimals, this difference is smaller in calculations with larger r₀.

Zoom In Zoom Out Reset image size

To conclude this discussion, Figure 10 shows how the q₄–r_max relation depends on the fragmentation parameters. For all x_m, larger planets grow in disks composed of stronger planetesimals. Once the collisional cascade begins, collisions remove more mass from weak planetesimals than from strong planetesimals. At identical evolution times, disks composed of weak planetesimals therefore have less mass in bins with small planetesimals. With oligarchs accreting most of their mass from small planetesimals, they do not grow as rapidly or as large when planetesimals are weak.

Zoom In Zoom Out Reset image size

Despite the differences in r_max, the median q₄ is fairly independent of the fragmentation parameters (Table 2). In the most massive disks with x_m = 0.33–3, calculations with stronger planetesimals have slightly larger q₄ than calculations with weak planetesimals: q₄ = 2.9–4.1 for f_w, q₄ = 3.2–4.3 for f_s, and q₄ = 3.1–4.2 for f_vs. These differences are smaller than the typical inter-quartile range of 0.5–1.0 when r_max ≈ 1000 km. In lower mass disks with x_m ≲ 0.1, calculations with stronger planetesimals tend to have shallower slopes, but the inter-quartile ranges are still much larger than the differences.

To summarize, Figures 6–10 demonstrate a clear relation between the size of the largest object and the slope of the cumulative size distribution for r ≳ 100 km. For a broad range of initial conditions and fragmentation parameters, q₄ is large (6–12) when r_max ≈ 300–500 km and small (2–5) when r_max ≈ 1000–3000 km. For any r_max, disks with smaller masses have smaller q₄ than more massive disks.

To understand how the rest of the size distribution evolves with r_max, we now consider the other three power-law exponents. We begin with the size distribution of 0.1–1 km objects in Section 3.2 and then discuss the evolution of 1–10 km (10–100 km) objects in Section 3.3 (Section 3.4).

The slope of the cumulative size distribution for 0.1–1 km objects, q₁, is a proxy for the evolution of objects smaller than r₀. As the largest objects grow, they remove planetesimals from this size range. During runaway growth and the early stages of oligarchic growth, oligarchs accrete only a small fraction of the mass in 0.1–1 km objects (KB08, KB10). Thus, the slope of this size distribution (q₁) should remain roughly constant with increasing r_max. Once the collisional cascade begins, destructive collisions among small objects remove material from this range of sizes; destructive collisions among larger objects add material. Although the evolution of the slope then depends on the relative importance of destructive collisions among large and small objects, the collisional cascade eventually removes all of the mass in 0.1–1 km objects from every annulus (KB08, KB10). Thus, q₁ should evolve dramatically during the cascade.

Figure 11 shows the variation of the median q₁ with r_max for all calculations with the f_w fragmentation parameters. Each panel illustrates results for a different value of r₀; the lower (upper) set of points corresponds to calculations with initial n_c∝r^−0.5 (n_c∝r⁻³). Colors indicate initial disk mass ranging from x_m = 3 (violet) to x_m = 0.01 (maroon) as in Figure 6. The remarkable overlap of colored points demonstrates that the evolution of the median q₁ is independent of initial disk mass.

Zoom In Zoom Out Reset image size

For calculations with r₀ = 1 km, there are two main features of the evolution of q₁ with r_max. As long as r_max ≲ 300–500 km, q₁ is constant. In the lower panel of Figure 11, q₁ remains fixed at its initial value as the largest objects grow from sizes of a few km up to roughly 300 km. Thus, the size distribution of 0.1–1 km objects "remembers" the initial size distribution throughout runaway and oligarchic growth.

Once the collisional cascade begins, the general evolution of q₁ is nearly independent of the initial size distribution and the initial size of the largest planetesimal. As the largest objects grow to sizes of ∼1000 km, all size distributions begin to converge on q₁ ≈ 2.0–2.5 (Table 3). After maintaining this limit for a small range in r_max, q₁ approaches a limiting value of 0.0–0.5. As the cascade removes nearly all of the mass from an annulus, q₁ remains at this limiting value.

The intermediate value of q₁ ≈ 2.0–2.5 results from collisional erosion (see the Appendix for a simple derivation). When the cascade begins, the ratio of the collision energy to the binding energy is Q_c/Q*_D ≈ 0.1–1. In this regime, collisions produce a modest net reduction in the mass of the larger object. Together with the mass of the smaller object, this material is lost as debris and accumulates in much smaller mass bins. Smaller objects are easier to erode. Thus, the size distribution evolves from the growth limit of roughly 3–4 to an intermediate limit of roughly 2.0–2.5.

The final limiting value of q₁ results from complete shattering (see the Appendix). As the cascade continues, the largest objects grow. Collision velocities among the smaller objects also grow. Once Q_c/Q*_D ≳ 2–3, nearly all of the mass of two colliding objects goes into debris. Within the 0.1–1 km size range, all objects are susceptible to shattering. When shattering of the small objects becomes important, q₁ approaches zero (Kenyon & Bromley 2004c).

For calculations with r₀ = 10–100 km, 0.1–1 km objects collide and merge into larger objects before the collisional cascade begins. Based on the evolution of the largest objects in Figures 6–10, growth tends to produce size distributions with q ≈ 2–4. Thus, when the initial size distribution is steep (n_c∝r⁻³), there is little evolution in q₁ for r_max ≲ 300–500 km (Figure 11, upper sets of points in the middle and upper panels). When the initial size distribution is shallow (n_c∝r^−0.5), slow growth leads to a gradual evolution from q₁ ≈ 0 to q₁ ≈ 2–3 (Figure 11, lower sets of points in the middle and upper panels).

Aside from this difference in the early evolution of q₁, these calculations also approach the two limiting values earlier in their evolution (at smaller r_max). When r₀ = 10–100 km, there is initially less mass in the 0.1–1 km range; thus, the collisional cascade modifies q₁ earlier. At the onset of the collisional cascade, calculations with r₀ = 1 km have most of their mass in 0.1–1 km planetesimals. Because it takes more collisional debris to change q₁ in these calculations, q₁ begins to change when r_max is larger. In these calculations, q₁ also has a larger dispersion as a function of initial disk mass and r_max.

Analyzing calculations with other fragmentation parameters yields nearly identical results (Table 3). The trend of the evolution is always independent of initial disk mass. For r₀ = 1 km, q₁ always remains close to its initial value. When r₀ = 10–100 km, q₁ remains roughly constant (when the initial q₁ is larger than 2) or slowly evolves to q₁ ≈ 2–3 (when the initial q₁ is smaller than 2). After the cascade commences, q₁ first evolves to the canonical slope for erosive collisions of 2.0–2.5 (at r_max ≈ 1000 km) and then to the canonical shattering slope of roughly 0.0–0.5 (at r_max ≈ 3000–5000 km). Although it always occurs at r_max ≈ 500–1500 km, the transition from the initial q₁ to the canonical q₁ depends on the amount of mass initially in 0.1–1 km objects. Disks with less (more) mass in these objects make the transition at smaller (larger) r_max and earlier (later) in evolution time.

Depending on the initial size of the largest object, the slope of the cumulative size distribution for 1–10 km objects, q₂, serves two distinct roles. For calculations with r₀ = 1 km, the median q₂ first measures changes in the slope as objects grow to sizes of 10 km and then tracks the changes in the size distribution of these objects throughout the rest of runaway/oligarchic growth and the collisional cascade. In these models, we expect the median q₂ first to decline with increasing r_max (as in the evolution of q₄ with r_max) and then to evolve to a slope characteristic of the collisional cascade. For calculations with larger objects, r₀ = 10–100 km, changes in the size distribution for 1–10 km objects should follow the changes in q₁ observed for the 0.1–1 km objects: the median q₂ should evolve slowly to ∼2–3 during runaway/oligarchic growth and then evolve more rapidly to a smaller value during the collisional cascade.

Figure 12 shows the variation of the median q₂ with r_max for all calculations with the f_w fragmentation parameters. The format is identical to Figure 11, with separate panels for each r₀, colors indicating the initial disk mass, and an upper and lower set of points corresponding to our two adopted initial size distributions. As in Figure 11, the evolution of the median q₂ as a function of r_max is independent of the initial disk mass.

Zoom In Zoom Out Reset image size

In each panel, the evolution of q₂ follows expectation. In the lower panel (r₀ = 1 km), the median q₂ declines from ∼5–6 when r_max ≈ 20 km to ∼4 for r_max ≳ 100 km. This evolution is independent of the initial size distribution. Once the collisional cascade begins, q₂ declines to ∼2.0–2.5 as r_max grows from 1000 km to 3000 km. At early times, the slope of the 1–10 km size distribution is similar to the slope of the large object size distribution (q₄); at late times, the slope evolves into one of the canonical slopes for the collisional cascade.

For larger r₀, the evolution of q₂ tracks features of the evolution of q₁. In the middle panel of Figure 12, q₂ follows much of the q₁ evolution for 1 km objects (Figure 11, lower panel). During runaway and oligarchic growth, q₂ remains constant with r_max. Once the cascade begins, q₂ evolves to the canonical value for the early stages of the collisional cascade q₂ ≈ 2.0–2.5. As r_max approaches 10⁴ km, q₂ remains constant at roughly 2.5 instead of dropping to ∼0.

In the upper panel, the evolution of q₂ for r₀ = 100 km similarly tracks the q₁ evolution for r₀ = 10 km. As the largest objects grow from ∼100 km to ∼1000 km, 1–10 km objects collide and merge into larger objects. Growth tends to evolve the slope to 2–4. For calculations with n_c∝r⁻³ (upper set of points), q₂ does not change. When the initial size distribution is concentrated into large objects (n_c∝r^−0.5), q₂ slowly approaches 2. Eventually, the collisional cascade begins; q₂ evolves to 2–2.5 and remains in this range as r_max reaches roughly 10⁴ km.

The ratio of the collision energy Q_c to the binding energy Q*_D allows us to understand the difference between the q₁ evolution for r₀ = 1 km and the q₂ evolution for r₀ = 10–100 km. When the collisional cascade begins, erosive collisions remove mass from the 0.1–1 km and the 1–10 km size ranges. During this evolution, both q₁ and q₂ evolve to 2.0–2.5. Once 3000 km ≲ r_max ≲ 10⁴ km, collisions among smaller planetesimals produce shattering. However, the larger planetesimals remain in the erosive regime. As a result, objects are completely removed from the ensemble of 0.1–1 km objects (producing q₁ = 0.0–0.5) and slowly eroded among the ensemble of 1–10 km objects (leaving q₂ at 2.0–2.5).

Results for calculations with the f_s and f_vs fragmentation parameters are nearly identical to those with the f_w fragmentation parameters (Table 4). When collisions produce mergers instead of debris, q₂ evolves from its initial value to the canonical "merger" limit of 2–4. When r₀ is 1 km, q₂ has no memory of the initial size distribution. Calculations with larger r₀ remember the initial q₂ and maintain this slope throughout much of the evolution. Once the cascade begins, q₂ always evolves to ∼2. Because 1–10 km objects are not completely shattered at any point in the collisional cascade, q₂ always remains larger than ∼2 and never drops to 0–1.

Based on our results for q₁ and q₂, the evolution of q₃, the slope for 10–100 km objects, should follow one of two paths. When the initial size of the largest planetesimal is 10 km or smaller, planetesimals must grow to fill the 10–100 km size bins. Growth produces a characteristic evolution of the slope from large values to ∼2–4 (Figure 6). As the evolution makes the transition from growth to the collisional cascade, the slope evolves to a canonical value of roughly 2. Because 10–100 km objects are harder to shatter than 1–10 km objects, the slope produced by the cascade should never fall below 2. When planetesimals have initial sizes larger than 10 km, collisional growth requires a long time to erase the initial size distribution. Eventually, growth and the cascade modify the size distribution, which should approach the canonical value of 2.

Figure 13 confirms these expectations for calculations with the f_w fragmentation parameters. Following the format of Figures 11 and 12, there is a panel for each r₀, colors indicating the range in x_m, and sets of points for each initial size distribution. Once again, the evolution is nearly independent of the initial disk mass.

Zoom In Zoom Out Reset image size

For calculations with r₀ = 1 km, growth initially produces large q₃ ≈ 4–6. Shallower initial size distributions and more massive disks tend to yield larger q₃. As objects grow, all calculations converge on a smaller range of q₃, ∼ 2.5–3.5, with smaller initial disk masses favoring smaller q₃. For r_max ≳ 1000–2000 km, q₃ is independent of initial disk mass.

When r₀ = 10 km, growth leads to a smaller range of initial median slopes, q₃ ≈ 5–6 for r_max ≈ 200 km. These median slopes are independent of the initial size distribution and the initial disk mass. As the largest objects grow to r_max ≈ 1000 km, q₃ reaches an approximate equilibrium at ∼4. Once the cascade begins, this slope gradually declines to q₃ ≈ 2.

This evolution for q₃ follows the evolution of q₂ for calculations with r₀ = 1 km. Independent of the initial size distribution, growth produces a characteristic size distribution which slowly evolves to shallower slopes as the largest objects grow. Eventually, the cascade changes the size distribution and the median q₃ becomes smaller. In the middle panel of Figure 13, this change occurs when r_max ≳ 3000 km. A similar, but less obvious, change in slope occurs for the r₀ = 1 km calculations in the lower panel.

When r₀ = 100 km, q₃ remembers the initial size distribution throughout most of the evolution (Figure 13, upper panel). In these calculations, the initial size distribution is unaffected by runaway/oligarchic growth and the early stages of the collisional cascade. Throughout these phases of the evolution, growth is slow and collisions have little impact on large objects. Thus, the size distribution is unchanged. During the late stages of the cascade, occasional collisions among large objects produce some debris which fills all mass bins in the 10–100 km size range. Independent of the initial size distribution, this debris drives the slope of the size distribution to the canonical value of q₃ = 2.0–2.5.

The evolution of q₃ is independent of the fragmentation parameters (Table 5). When the planetesimals are initially small, q₃ first evolves from large values (∼4–6) to small values (∼3). Once the cascade begins, q₃ evolves to 2.0–2.5. In this size range, the collisional cascade has less impact than at smaller sizes. Thus, the transition from the q₃ ≈ 3–4 characteristic of runaway/oligarchic growth to q₃ ≈ 2.0–2.5 occurs at larger r_max (and later evolution time).

At 15–150 AU, the evolution of planetesimals follows a similar path at every distance from the central star. Initially, collisional damping and gas drag reduce the velocities of small objects. Among larger objects, low-velocity collisions produce mergers instead of debris. Along with dynamical friction and viscous stirring, damping produces large gravitational focusing factors; this evolution leads to runaway growth of the largest planetesimals. As these protoplanets grow, they stir up the leftover smaller planetesimals. Eventually, collisions among the leftovers produce debris instead of mergers. A collisional cascade ensues, which grinds planetesimals with sizes of 1–3 km and smaller to dust.

Despite the similarity in the evolutionary path, there is considerable diversity in the overall outcome. For identical starting conditions, the stochastic nature of coagulation leads to a broad range in the rate protoplanets grow, in the maximum radius of protoplanets, in the distribution of sizes as a function of time, and in the amount of mass lost through the collisional cascade. Although this diversity does not mask clear trends in the median outcome as a function of initial conditions (Kenyon & Bromley 2008, 2010), the chaotic nature of the evolution prevents associating a single outcome with a single set of starting conditions.

To quantify the overall evolution of the size distribution of icy planetesimals at 15–150 AU, we define four slopes, q₁–q₄, which measure the best-fitting slope of a power-law size distribution between 0.1–1 km (q₁), 1–10 km (q₂), 10–100 km (q₃), and 100 km to r_max (q₄). For all r_max ≳ 200 km, we derive the median slopes and the inter-quartile range as a function of r_max. We use median q₁–q₄ and the inter-quartile range to track the evolution of the size distribution and to measure the range of outcomes as a function of initial conditions.

Based on these slopes, we infer several general features of the evolution (see also Tables 2–5).

• 1.  
  For all distances from the central star, the growth of icy planets leads to size distributions with q ≈ 2–4 (see also Kenyon & Luu 1999b; Kenyon 2002; Goldreich et al. 2004; Schlichting & Sari 2011). This size distribution results from large objects accreting primarily small objects rather than other large objects. As the largest objects grow, the slope evolves from q₄ ≈ 6–14 when r_max ≈ 200 km to q₄ ≈ 2–4 when r_max ≳ 1000 km.
• 2.  
  The collisional cascade produces shallower slopes in the size distribution, q ≲ 2.0–2.5. During the cascade, several processes set the size distribution for objects with r ≳ 0.1 km (see also O'Brien & Greenberg 2003; Kenyon & Bromley 2004c; O'Brien et al. 2005; Pan & Sari 2005): (1) growth by mergers (q ≈ 3), (2) debris production (q ≈ 2.5–3, e.g., Dohnanyi 1969; O'Brien & Greenberg 2003), (3) erosion (q ≈ 2.0–2.5; see the Appendix), and (4) shattering (q ≈ 0.0–0.5, e.g., Kenyon & Bromley 2004c, and the Appendix). At the smallest sizes (r ≲ 0.01–0.1 km), debris production dominates the size distribution and q ≈ 3.5–4 (O'Brien & Greenberg 2003; Kenyon & Bromley 2004c, 2005, 2008, 2010; O'Brien et al. 2005; Krivov et al. 2005; Benavidez & Campo Bagatin 2009; Fraser 2009). At larger sizes, erosion slowly removes mass from all objects and produces q ≈ 2.0–2.5. Once the collisional cascade reaches its peak (r_max ≳ 3000 km), shattering of 0.1–1 km objects leads to q ≈ 0.0–1.0.
• 3.  
  The evolution of the large object size distribution is independent of radial distance from the central star and is generally independent of the fragmentation parameters. However, the evolution of q₄ is sensitive to the initial size of the largest planetesimal (r₀) and the initial mass of the disk (x_m). Lower mass disks evolve over a smaller range of q₄. Calculations with larger planetesimals tend to produce larger median q₄.
• 4.  
  For any r_max, the range in q₄ is fairly large. At small sizes (r_max ≈ 200–500 km), the inter-quartile range about the median q₄ is roughly 2. This range drops to roughly 1 at r_max ≳ 1000 km. For a typical q₄ ≈ 3, roughly 50% of calculations with the same initial conditions yield q₄ = 2.5–3.5. The rest have q₄ = 2–2.5 and q₄ = 3.5–5.
• 5.  
  Among the intermediate mass objects with r = 1–100 km, the evolution of the slope of the size distribution depends on the initial size of the largest planetesimal.
• 6.  
  When r₀ is small, growth produces objects with r = 1–100 km. Growth erases knowledge of the initial size distribution; the slope evolves to the standard "growth slope," q₂ ≈ q₃ ≈ 2–4. Once the collisional cascade begins, the slopes evolve to the canonical value of 2 for erosive collisions.
• 7.  
  When r₀ is large, growth cannot change the initial slope of the size distribution. Thus, q₂ and q₃ remember the initial q until the onset of the collisional cascade. Once the cascade begins, erosive collisions lead to q₂ ≈ q₃ ≈ 2.0–2.5.

Current theories for the origin of the solar system broadly explain the nature of TNOs as a several step process. During the early evolution of the protosolar disk, dust grains within the disk grow into icy planetesimals (e.g., Weidenschilling 2006; Chiang & Youdin 2010) and then into icy protoplanets (e.g., Lissauer 1987; Kenyon & Luu 1999a; Kenyon & Bromley 2008, 2010). Four protoplanets become gas giants (Pollack et al. 1996; Alibert et al. 2005; Ida & Lin 2008; Lissauer et al. 2009; Bromley & Kenyon 2011a). Sometime after the gas giants reach roughly their final masses, dynamical interactions among the gas giants, the much less massive icy protoplanets, and any leftover planetesimals result in the ejection or migration of icy objects into the trans-Neptunian region and the Oort Cloud (e.g., Malhotra 1993; Hahn & Malhotra 1999; Morbidelli et al. 2004; Tsiganis et al. 2005; Gomes et al. 2005). Once dynamical evolution begins, lower surface densities and larger orbital velocities among the TNOs reduce growth rates and promote collisional erosion (e.g., Stern & Colwell 1997b; Kenyon & Luu 1999a; Kenyon & Bromley 2004c, 2008, 2010; O'Brien et al. 2005; Fraser 2009). Thus, dynamical evolution ends the growth of TNOs.

Despite the uncertainty in the theory, this picture constrains the initial mass of the protosolar disk and the likely range of formation times for TNOs. The protosolar disk probably has properties similar to the disks observed in nearby star-forming regions such as Ophiuchus, Orion, and Taurus-Auriga (e.g., Williams & Cieza 2011). Infrared and radio observations demonstrate that these systems lose their gas and dust on timescales of 1–10 Myr (Haisch et al. 2001; Currie et al. 2009; Kennedy & Kenyon 2009; Williams & Cieza 2011). For the solar system, these data constrain the initial disk mass, 0.001–0.1 M_☉, and the timescale, 1–10 Myr, for the formation of gas giant planets.

Dynamical calculations of the solar system place additional limits on the initial disk mass and the formation time. Several results suggest some TNOs formed at 20–25 AU and then migrated outward with the giant planets (e.g., Malhotra 1993, 1995; Hahn & Malhotra 1999; Levison & Morbidelli 2003; Morbidelli et al. 2008). To explain the orbital architecture of the giant planets, the TNOs and the Oort Cloud, these calculations rely on a massive disk of solids with x_m ≳ 1 at 15–30 AU (Hahn & Malhotra 1999; Levison & Morbidelli 2003; Morbidelli et al. 2008). Although these calculations place little direct constraint on the formation time, Figures 6–10 demonstrate that more massive disks have larger q₄ at fixed r_max. Thus, relying on a massive disk to enable outward migration places indirect con-straints on any combination of formation time, r_max and q₄.

Various lines of evidence associate migrating gas giant planets with an increased rate of cratering during the Late Heavy Bombardment (hereafter LHB; e.g., Levison et al. 2001; Gomes et al. 2005; Strom et al. 2005). Although constraints on the cratering history prior to this epoch are weak (e.g., Chapman et al. 2007), detailed analyses suggest the intense period of bombardment ended ∼500–700 Myr after the formation of the solar system (e.g., Tera et al. 1974; Hartmann et al. 2000; Stöffler & Ryder 2001). It is unlikely that the gas giants migrated after the LHB. Thus, the end of the LHB sets an upper limit on the formation time for TNOs.

For our analysis, we concentrate on the formation of TNOs before the gas giants dynamically rearrange the solar system. The initial conditions chosen for our calculations broadly match the range of initial disk masses and semimajor axes applicable to TNO formation. By following the evolution of TNOs for 1–10 Gyr, we cover the most likely range of formation times, ∼10–700 Myr. Thus, our calculations can provide clear tests of the coagulation picture for TNO formation.

Recent observations place excellent limits on the masses and radii of the largest TNOs. All of the major dynamical classes (e.g., Gladman et al. 2008) have at least one large object with radius r ≈ 450–1300 km (e.g., Sheppard et al. 2011; Petit et al. 2011). All-sky surveys are now nearly complete to an R-band magnitude R ≈ 21, the approximate brightness of an object with radius r ≈ 350 km and an albedo of 10% at 42 AU (e.g., Schwamb et al. 2010; Sheppard et al. 2011). Thus, the sample of largest TNOs is probably complete at 40–50 AU. For large TNOs with binary companions, period measurements yield reliable masses of roughly 1.5 × 10²⁴ g for Quaoar (Fraser & Brown 2010) to roughly 1.5 × 10²⁵ g for Pluto (Buie et al. 2006) and 1.6 × 10²⁵ g for Eris (Brown & Schaller 2007). In our calculations, this mass range corresponds to r ≈ 600–1400 km. Thus, a typical r_max = 1000 ± 400 km matches the current population of the largest TNOs.

For deep surveys with sufficient observations, the luminosity function—written as σ(m), the surface number density of objects brighter than magnitude m—is reasonably well fit by a double power law,

$$\begin{equation} \sigma (m) = \sigma _{m_0} {1 + 10^{(\alpha _L - \alpha _S)(m_{{\rm eq}} - m_0)} \over 10^{\alpha _L(m - m_0)} + 10^{(\alpha _L-\alpha _S)(m_{{\rm eq}} - m_0) - \alpha _S (m - m_0)}}, \end{equation} \tag{ 9 }$$

where $\sigma _{m_0}$ is a normalization factor, m₀ is a reference magnitude, m_eq is the brightness where the two power laws merge, α_L is the slope of the size distribution at large sizes, and α_S is the slope of the size distribution at small sizes (e.g., Petit et al. 2006, 2008, 2011; Fuentes & Holman 2008; Fuentes et al. 2009; Fraser et al. 2008, 2010; Fraser & Kavelaars 2009; Schwamb et al. 2010; Sheppard et al. 2011).

Converting observations of the surface number density to constraints on the size distribution requires additional information (e.g., Fraser et al. 2008; Petit et al. 2008). If the cumulative size distribution and the distribution of albedos follow power laws with exponents q (for the size distribution) and β (for the albedo),

$$\begin{equation} q = 5 (\alpha - \beta / 2) + 1. \end{equation} \tag{ 10 }$$

Although TNOs display a broad range of albedos (e.g., Stansberry et al. 2008; Brucker et al. 2009), current data suggest a roughly constant albedo of 2%–20% for r ≈ 30–500 km and a larger albedo of 20% to nearly 80% for r ≳ 700 km (Figure 3 of Stansberry et al. 2008). Most objects comprising the luminosity function are faint; adopting β ≈ 0 is therefore a reasonable assumption. The slope of the size distribution is then q = 5α + 1.

Recent analyses indicate a range in the slope of the luminosity function at large sizes. Several deep surveys yield α_L ≈ 0.70–0.76, implying q_L ≈ 4.5–4.8 for β = 0 (Fuentes & Holman 2008; Fuentes et al. 2009; Fraser & Kavelaars 2009). Fraser et al. (2010) divided their sample of the classical Kuiper Belt into a set of cold objects with inclination i ⩽ 5° and a set of hot objects with i ⩾ 5°. Their results suggest α_L ≈ 0.80 (q_L ≈ 5) for the cold objects and α_L ≈ 0.40 (q_L ≈ 3) for the hot objects (see also Levison & Stern 2001; Bernstein et al. 2004; Fuentes & Holman 2008; Fuentes et al. 2011). Sheppard et al. (2011) also prefer a shallower slope for the large objects.

Other observations also distinguish the hot and cold populations within the Kuiper Belt. The largest cold KBOs are smaller (Levison & Stern 2001) and have a larger fraction of binary companions (Noll et al. 2008) than the largest hot KBOs. The cold KBOs also have larger albedos (Brucker et al. 2009). Although there is no consensus, some results suggest that the cold KBOs are considerably redder than the hot KBOs (e.g., Tegler & Romanishin 2000; Trujillo & Brown 2002; Peixinho et al. 2004, 2008; Gulbis et al. 2006).

This diversity suggests that the cold and hot KBOs formed in different locations or have had different collisional histories (Levison & Stern 2001; Malhotra 1993; Gomes 2003; Levison & Morbidelli 2003; Levison et al. 2008). In current dynamical models, the observed (a, e, i) distributions of the hot KBOs are a result of formation at 15–25 AU and subsequent migration and scattering during the migration of the giant planets. In this picture, the cold population formed at a larger distance, perhaps in situ (see Batygin et al. 2012, and references therein).

To assign r_max for the hot and cold populations, we use results from direct mass estimates of individual KBOs (e.g., Grundy et al. 2009, 2011; Fraser & Brown 2010) and deep probes of the KBO luminosity function (e.g., Fraser & Kavelaars 2009; Fraser et al. 2010; Petit et al. 2011). These analyses consistently demonstrate that the largest hot KBOs are systematically larger than the largest cold KBOs (see also Levison & Stern 2001). We infer a typical "size ratio" ≳1.5–2, suggesting r_max ≈ 1000 km for the hot objects and r_max ≲ 600 km for the cold objects. Although the true r_max for the cold population might be less than 600 km, results for (q₄, r_max) = (5, 600 km) are a reasonable proxy for any size in the 400–600 km ((0.4–1.3) × 10²⁴ g) range. After deriving the constraints these observations place on the models, we will consider alternative models for a smaller r_max = 100–200 km.

To use the observational and theoretical constraints to test our calculations, we assume that our q₄ is equivalent to the observed q_L. For our purposes, the q₄ derived for the cold KBOs is identical to the q₄ inferred from most deep surveys. Thus, we consider the implications of q₄ = 3 and q₄ = 5 for r_max = 600–1400 km. After a basic application of our models to these data in Section 4.3.1, we derive a set of "success fractions," which measure the fraction of models that match the observations of q₄ and r_max (Section 4.3.2). This analysis yields two plausible models for the formation of the cold and hot TNOs: (1) ensembles of 1 km planetesimals in a massive disk with x_m = 0.3–3 and (2) ensembles of 10 km planetesimals in a less massive disk with x_m = 0.03–0.3. To select the "best" model, we consider the range of formation times at 20–50 AU (Section 4.3.3). These results favor ensembles of 1 km planetesimals as the most likely starting point for the current populations of hot and cold TNOs.

4.3.1. Basic Tests

4.3.2. Success Fractions

Having established rough limits on the initial disk mass and planetesimal size, we now quantify the range of calculations which match the observational constraints on r_max and q₄. To make these choices, we consider ensembles of results for each set of initial conditions and fragmentation parameters. From Table 1, there are three ensembles with r₀ = 10 km and 100 km, and five ensembles with r₀ = 1 km. For each of these ensembles, we collect all pairs of (r_max, q₄) and bin these data in increments of 0.1 in log r_max. Within each r_max bin, stochastic effects produce a broad range of q₄. Thus, we identify the fraction of calculations within each bin that match an adopted target q₄ as a function of r_max. To span the observed range, we consider two target slopes, q₄ = 3 ± 0.25 and 5 ± 0.25, for r_max = 600 km, 1000 km, and 1400 km.

To judge the uncertainties of the derived "success fractions," we examine the variation of the success fraction as a function of r_max. For a single set of initial conditions (x_m, r₀, q_init, f_i), the success fraction usually has an obvious trend (from larger to smaller success fraction or from smaller to larger success fraction) as a function of r_max. We set the uncertainty in a single success fraction as the dispersion around this trend. Typically, this dispersion (and thus the adopted uncertainty) is 0.05–0.10.

Figure 14 shows results for a target q₄ = 5 ± 0.25. The three sets of vertical panels plot the fraction of successful models for r_max = 600 km (left panels), 1000 km (middle panels), and 1400 km (right panels). The three sets of horizontal panels plot the success fraction for calculations with r₀ = 1 km (lower panels), 10 km (middle panels), and 100 km (upper panels). Within each panel, colored lines show the variation of the success fraction as a function of the initial disk mass. Different colors indicate different initial mass distributions and fragmentation parameters as listed in the figure caption (see also Table 1).

Zoom In Zoom Out Reset image size

Clearly, a broad range of calculations match the target slope at rates of roughly 20%–50%. For r₀ = 1 km and r_max = 600 km (lower left panel), calculations with any set of initial conditions and low disk masses (x_m ≲ 0.1) fail to match the target slope. Models with larger initial disk masses (x_m ≳ 0.1) have larger success fractions. Within this range of x_m, calculations with an initially large range of planetesimal sizes (green, orange, and magenta curves) are more successful than calculations with a single planetesimal size (violet and blue curves). As r_max increases to 1000 km (lower middle panel) and to 1400 km (lower right panel), the success fraction systematically decreases. In these panels, calculations with a single planetesimal size are somewhat more successful than those with a range of planetesimal sizes. Once r_max = 1400 km, however, few calculations successfully match the target slope.

Calculations with larger r₀ also match the target slope. For r₀ = 10 km and 100 km, success fractions of 20%–40% are common. Among the r₀ = 10 km panels, there is a clear trend of larger success fraction with increasing initial disk mass. Although there is a small trend of more success for calculations with a single initial planetesimal size, the trend is much weaker than the trend with r_max. Among the models with r₀ = 100 km, less (more) massive disks yield larger success ratios when r_max is small (large). When r_max is small, calculations with a broad range of initial planetesimal sizes often match the data; results for a single initial planetesimal size rarely match the data. This trend reverses at large r_max, when calculations with a single planetesimal size are often more successful.

Figure 15 repeats Figure 14 for a target q₄ = 3 ± 0.25. In these examples, calculations with 1 km planetesimals are much more successful than those with large planetesimals. For any r_max (any vertical column of panels in the figure), the success fractions decline monotonically from 20%–40% at r₀ = 1 km to 0% for r₀ = 100 km. Although there is little trend in the success fraction with r_max for calculations with r₀ =10–100 km, there is a modest trend for r₀ = 1 km: smaller r_max requires disks with smaller masses to match the target slope. In all panels for r₀ = 1 km, there is little trend with the initial size distribution of planetesimals or the fragmentation parameters: all of the initial conditions yield success ratios of 20%–40% for x_m = 0.03–3.

Zoom In Zoom Out Reset image size

To refine these tests, we now include other constraints on the hot and cold populations. In the dynamical models, hot KBOs formed closer to the Sun and then migrated out to their current orbits. In the context of our coagulation models, objects closer to the Sun grow faster (Figure 4) and produce shallower size distributions (Figure 5) over the same evolution time. Both of these features of collisional evolution agree with the observations: the hot KBOs are larger and have a shallower size distribution than the cold KBOs.

To find a combination of starting conditions that match the observations, we consider models that match (q₄, r_max) = (3, 1000 km) for the hot objects and (5, 600 km) for the cold objects. Using results in Figures 14 and 15, we rule out any models starting with 100 km planetesimals. These calculations never match the shallow slope for the hot objects. However, two options with 1 km and 10 km planetesimals are equally likely.

• 1.  
  10 km planetesimals in a low-mass disk (x_m = 0.03–0.3).
• 2.  
  1 km planetesimals in a massive disk (x_m = 0.3–3).

In both cases, models with a single planetesimal size or a range of planetesimal sizes match the constraints reasonably well. Models with a single planetesimal match more often when r₀ = 10 km; models with a range of sizes match more often when r₀ = 1 km.

4.3.3. The Formation Time

To select the "better" alternative among these two options, we add a final constraint—the formation time—to our analysis. The dynamical models require TNO formation at 20–40 AU sometime between 1–10 Myr and 500–700 Myr. The smaller of these two constraints is not useful. Only a few models produce TNOs in 10 Myr at 20 AU; none of our calculations produce TNOs in 10 Myr at 40–50 AU. The longer constraint is very useful. Many calculations require more than 1 Gyr to produce 600–1000 km objects at 40–50 AU. Thus, the dynamical models rule out these coagulation models.

To select between the allowed combinations of r₀ and x_m, we consider the evolution time required to form 1000 km objects as a function of semimajor axis. Figure 16 shows results for calculations with the f_w fragmentation parameters. Results for other calculations are similar. Individual panels plot results for r₀ = 1 km (lower panels), r₀ = 10 km (middle panels), and r₀ = 100 km (upper panels) and two initial size distributions (left and right panels). Colors indicate the initial disk mass, ranging from x_m = 3 (violet) to x_m = 0.01 (maroon) as in Figure 6.

Zoom In Zoom Out Reset image size

These results follow standard expectations for coagulation calculations. Formation times always scale with the initial radius of the largest planetesimal (e.g., Lissauer 1987; Kenyon & Luu 1998; Kenyon & Bromley 2010). As r₀ increases from 1 km to 100 km, the clusters of points shift to longer formation times. In a disk with a surface density Σ∝x_ma^−3/2, the formation timescales as t_f∝x⁻³_ra³ (Goldreich et al. 2004; Kenyon & Bromley 2008). As a increases, formation times in each panel grow rapidly. Similarly, models with smaller x_m have longer formation times.

The solid horizontal line in each panel indicates the approximate end of the LHB, 500–700 Myr. Applying the additional constraint that TNOs form at a ≈ 20–50 AU yields constraints on all of our calculations:

$$\begin{equation} x_m \gtrsim \left\lbrace \begin{array}{@{}cl}0.01\, (r_0 / {\rm 1\, km})^{1/2} (a / {\rm 20\, AU})^3 & {\rm all\ planetesimal\ sizes}, \\ 0.03\, (r_0 / {\rm 1\, km})^{1/2} (a / {\rm 20\, AU})^3 & {\rm one\ planetesimal\ size}. \end{array} \right. \end{equation} \tag{ 11 }$$

When r₀ = 100 km, these constraints are severe: in situ formation of TNOs at 40 AU requires x_m ≳ 1–3. Thus, calculations with very large planetesimals are viable only if the initial mass of the disk exceeds the minimum mass solar nebula. For smaller r₀, calculations with smaller x_m yield 1000 km TNOs in 500 Myr at 40 AU.

Requiring TNOs to form in 10 Myr places more stringent limits on the initial disk mass and the range of semimajor axes. In this case, TNOs can only form at a ≈ 20 AU for massive disks (x_m = 1–3) composed of 1–10 km planetesimals. Rapid formation times exclude all combinations of models with r₀ larger than 10 km.

These conclusions are independent of the fragmentation parameters. When other initial conditions are held fixed, formation times for icy planets with r_max ≈ 1000 km are similar for all three sets of fragmentation parameters. Although the size of the largest object varies with f_i (Figure 10), the variation of q₄ with r_max is similar for calculations with weak, strong, and very strong planetesimals. Thus, the formation time places strong constraints on the initial disk mass and the initial sizes of planetesimals.

From Equation (11), formation at 40–50 AU places the most severe constraints on the initial disk mass. For 1–10 km planetesimals, we adopt the less restrictive model with all planetesimal sizes. Then

$$\begin{equation} x_m \gtrsim \left\lbrace \begin{array}{@{}cl}0.1 - 0.2 & {\rm 1\ km\ planetesimals}, \\ 0.3 - 0.6 & {\rm 10\ km\ planetesimals}. \end{array} \right. \end{equation} \tag{ 12 }$$

Because coagulation models with 10 km planetesimals match the observed range in q₄ and r_max only for x_m ≈ 0.03–0.3, the combined constraints of the formation time, the size of the largest object, and the slope of the large object size distribution strongly favor calculations with 1 km planetesimals.

Some observations suggest that the largest cold KBOs have r_max ≈ 100–200 km instead of the r_max = 400–600 km used in our analysis (e.g., Brucker et al. 2009; Grundy et al. 2009, 2011). To investigate how a smaller r_max impacts our results, we now consider the range of coagulation models which yields an observed q = 5 ± 0.25 for r_max = 100–200 km. For this analysis, we use q₃ as a proxy for the slope of the size distribution for cold KBOs smaller than r_max. For r_max = 100–200 km, q₄ provides a poor measure of the slope below 100 km; q₃ provides an accurate measure of the slope over 20–200 km.

To identify coagulation models which match the observed q₃ for a smaller r_max, we derive success fractions for all sets of initial conditions in Table 1 assuming a target q₃ = 5 ± 0.25 and r_max = 100–200 km. For calculations with r₀ = 100 km and r_max = 100–200 km, q₃ is identical to the initial slope (see Figure 13). Because our adopted initial slope is q₃ = 0.17 or 3.0, our calculations never yield the observed q₃ = 5. Thus, all of our calculations with r₀ = 100 km fail to match the target slope.

Calculations with r₀ = 100 km and a steeper initial size distribution probably cannot match the target slope. When the initial q₃ ≳ 3, most of the mass is in the smallest objects. For modest initial surface densities, x_m ≳ 0.1, growth of these small objects should produce shallower slopes with q₃ ≈ 3–4.

For all other initial conditions, we derive success fractions of 10%–50%. When r_max is much smaller than 500 km, the collisional cascade has no impact on the slopes of the size distributions. Thus, success fractions for a target (q₃, r_max) = (5, 100–200 km) are independent of the fragmentation parameters. For calculations with all sizes or one size of planetesimals, there are weak trends of increasing success fractions in more massive disks. Lower mass disks produce slower growth, which tends to yield shallower size distributions. However, this trend is weak for small r_max.

Despite this lack of trends with f_i or x_m, there are trends with the initial size and size distribution of planetesimals. For r₀ = 10 km, a larger fraction of our models match the target q₃ for r_max = 150–200 km (40%–50%) than for r_max = 100–150 km (10%–30%). Calculations with r₀ = 1 km yield smaller success fractions for r_max = 150–200 km (10%–35%) than for r_max = 100–150 km (35%–50%). In both cases, models with a single planetesimal size have somewhat smaller success fractions, but the differences are not significant.

Thus, coagulation models match the observed q₃ = 5 for the cold KBOs when r_max = 100–200 km. Calculations with r₀ = 1 km (10 km) yield the target slope more often when r_max = 100–150 km (150–200 km). With little preference for the disk mass or the fragmentation parameters, this analysis does not identify a clear preference for r₀ or x_m.

The most basic conclusion of our analysis is that the robust relation between r_max and q₄ allows coagulation models to make clear predictions of (q₄, r_max) for any deep survey of TNOs. Unambiguous relations between q₃ and r_max also enable clear predictions for small TNOs. In principle, observations of (q₄, r_max) or (q₃, r_max) can falsify current coagulation models. In practice, current observations agree with the predictions.

Using observational results from deep surveys of TNOs together with conclusions derived from dynamical models of the solar system, we place good constraints on coagulation models for the origin of TNOs throughout the solar system. Figure 17 summarizes the main results of our analysis. As we add more observational and theoretical constraints to the analysis, the set of coagulation calculations which "match" the constraints narrows. Thus, a complete set of observations for all TNO classes provides a rigorous test of coagulation models.

Zoom In Zoom Out Reset image size

To summarize current limits on the models, we step through the observational and theoretical constraints. For a single size distribution of TNOs, we derive constraints on coagulation models from q₄ and r_max (top two rows of Figure 17).

• 1.  
  With a success rate of 20%–50%, calculations with a broad range of initial conditions and fragmentation parameters can match a target combination (q₄, r_max) which is consistent with observations of TNOs.
• 2.  
  For a specific target (q₄, r_max), the range of success rates allows us to rule out some ranges of initial conditions and fragmentation parameters.
• 3.  
  When the target (q₄, r_max) = (5, 1400 km), calculations with large planetesimals are more successful than those with small planetesimals (Figure 14).
• 4.  
  Smaller target slopes, q₄ ≈ 3, appear to rule out calculations with large planetesimals (Figure 15).
• 5.  
  In all cases, models with large disk masses are more successful than those with small disk masses. Thus, models with massive disks composed of large planetesimals are much more likely to match the target than low-mass disks composed of small planetesimals.

Adding in the dynamical constraint that the hot and cold populations originate in different locations of the protosolar disk yields more severe limits on the models. Requiring (q₄, r_max) = (3,1000 km) for the hot KBOs and (q₄, r_max) = (5,400–600 km) for the cold KBOs (the third row of Figure 17) leads to these conclusions.

• 1.  
  Models with 100 km planetesimals cannot match the data for any set of initial conditions.
• 2.  
  Models with 10 km planetesimals match the data for low-mass disks, x_m = 0.03–0.3.
• 3.  
  Models with 1 km planetesimals match the data for massive disks, x_m = 0.3–3.

Finally, using the constraint on formation time eliminates models with 10 km planetesimals. Forming TNOs in 500–700 Myr at 20–50 AU requires massive disks with x_m = 0.3–3 (the fourth row of Figure 17). Thus, we strongly favor models starting with 1 km planetesimals.

Adopting a different observational constraint, r_max = 100–200 km, for the cold KBOs results in similar conclusions. The observed slope for the hot population, q₄ = 3, drives us to models with small planetesimals in massive disks.

Our analysis is the first attempt to apply multiannulus coagulation calculations to modern observations of TNOs. The results indicate that TNOs form within a massive disk composed of small, 1 km planetesimals. This conclusion is similar to recent results on the origin of asteroids, where an initial population of 0.1 km planetesimals seems sufficient to explain the size distribution of asteroids for sizes of 10–100 km (Weidenschilling 2011, for a different conclusion, see Morbidelli et al. 2009). Although we did not consider formation of TNOs from 0.1 km planetesimals, calculations with smaller planetesimals should yield similar results for the size distributions at 10–1000 km. Thus, current data are also consistent with TNO formation from 0.1 km planetesimals.

The similarity of the size distributions for the TNOs and the trojans of Jupiter (Jewitt et al. 2000; Jedicke et al. 2002) and Neptune (Sheppard & Trujillo 2010) suggest a common origin in a massive disk composed of 1–10 km planetesimals (for a discussion of the dynamical origin of these populations, see Morbidelli et al. 2009). Unless the initial surface density of the disk is small (x_m ≲ 0.03–0.1; see Equation (11)), collisional evolution erases the signatures of many initial conditions on 10–500 Myr timescales. Dynamical models for the orbital architecture require massive disks (Morbidelli et al. 2008). During the dynamical re-arrangement of the giant planets, the large relative velocities among planetesimals leads to a collisional cascade which modifies the size distribution at a ≈ 15–30 AU considerably for r ≲ 100 km (Figures 11–13).

The main alternative to this picture is that the observed size distributions for the cold KBOs and the trojans of Jupiter and Neptune are "primordial," relatively unchanged since the initial process that formed planetesimals from mm-sized or smaller particles (e.g., Sheppard & Trujillo 2010). To evaluate this possibility, we isolate coagulation calculations where an initial size distribution of 100 km and smaller objects produces 110 km or smaller objects prior to the end of LHB. For models with a single size or a range of sizes, calculations with x_m ≲ 0.0003 (a/40 AU)³ satisfy this constraint. Models with larger x_m yield objects larger than 110 km. For primordial objects formed in situ at a ≈ 20–30 AU, this limit implies an initial mass of 0.003 M_⊕, roughly three times larger than the current mass in cold KBOs (e.g., Brucker et al. 2009). Although coagulation models do not directly preclude that cold KBOs and trojans form in a very low mass region of the disk, it seems unlikely that dynamical mechanisms can populate three distinct volumes of the solar system with so little reduction in the total mass (e.g., Morbidelli et al. 2005; Levison et al. 2008; Batygin et al. 2011; Wolff et al. 2011). Thus, we conclude that the observed size distributions are not primordial.

Comparisons between our calculations and the current size distributions of TNOs make several important assumptions. We assume that the hot and cold populations in the Kuiper Belt are representative of all TNOs. Current data do not contradict this assertion (e.g., Sheppard & Trujillo 2010; Petit et al. 2011). Data from ongoing and planned all-sky surveys will eventually yield robust estimates for (q₄, r_max) for all dynamical classes among the TNOs and will likely improve estimates for the trojans of Jupiter/Neptune. If other dynamical classes of TNOs or the trojans of Jupiter/Neptune have different combinations of (q₄, r_max) than analyzed here, our approach provides a framework to test coagulation models with these new data.

Detailed tests of our coagulation calculations also assume that large-scale rearrangement of the planets in the solar system has little impact on relations between r_max and q₄. We do not expect this evolution to produce significant changes in q₄. Before dynamical interactions with the gas giants become important, gravitational focusing factors are already small. Collision rates are in the orderly regime, where all objects collide at comparable rates. Thus, q₄ should not evolve considerably once large-scale dynamical evolution is important.

Throughout large-scale dynamical evolution, we expect significant changes to q₁–q₃ as a result of collisional evolution. Because collisions remain in the orderly regime, these changes should follow the paths outlined in our discussion to Figures 11–13. Depending on the timing of the start of the dynamical evolution, q₂ and q₃ for dynamically hot populations should evolve to roughly 2 or remain roughly constant at 2. If collision rates are large enough, q₁ should evolve from roughly 2 to roughly 0.

Some change in r_max during the epoch of large-scale dynamical evolution seems likely. Because collisions occur at high velocities, it is unlikely that r_max increases. Low collision rates imply little collisional erosion. However, the occasional high-velocity, "hit-and-run" collision (Asphaug et al. 2006) could remove the outer shell of material from a large object and leave the denser core unchanged. This mechanism is one way to explain the relatively high density of Quaoar (Fraser & Brown 2010). Dynamical evolution can also eject a few of the largest TNOs from the solar system. Because the number of large objects is always small, a few ejections can reduce r_max and lead to variations in r_max among different TNO populations.

Deriving the outcome of the dynamical evolution requires codes which combine statistical calculations of small objects as in this paper with the N-body calculations of gas giant planet formation as in Bromley & Kenyon (2011a) and calculations of the migration and scattering of small objects as in Charnoz & Morbidelli (2007), Levison et al. (2008, 2010), and Bromley & Kenyon (2011b). Following the relative timing of the formation of gas giants and TNOs within the same disk will provide better limits for r_max, q₄, and possible locations for the formation of TNOs within the protosolar disk. Combining collisional growth with scattering will yield better predictions for the break radius and the slope of the size distribution below the break. Our encouraging results from pure coagulation models suggest that a more complete model will enable a better understanding for the origin of TNOs and other large objects throughout the solar system.

As oligarchic growth proceeds, cratering collisions gradually dominate physical interactions among 0.1–10 km planetesimals. In this regime, collisions between equal mass objects produce the most debris. For a power-law cumulative size distribution with q ≳ 3, small planetesimals are more likely to collide with one another than with large oligarchs (Kenyon & Bromley 2004a; Goldreich et al. 2004). Thus, cratering leads to a gradual erosion of small planetesimals (e.g., Tanaka & Ida 1996; Kobayashi & Tanaka 2010).

During this erosive phase, collisions between equal mass objects produce a single object with mass comparable to the target and copious debris. For these objects, the slope of the size distribution then depends on the rate collisions convert planetesimals into debris and the rate of debris production over a range of sizes. To estimate the rate planetesimals are lost, we consider the collision rate, nσvf_g. During oligarchic growth, small planetesimals have similar orbital eccentricities and f_g ≈ 1. For collisions between equal mass objects, the rate small planetesimals are destroyed is dN/dt∝N²r², where N is the number of planetesimals in a size range and r is a typical radius. Setting N∝r^−q,

$$\begin{equation} d {\rm log} N / dt \approx r^{2-q}. \end{equation} \tag{ A1 }$$

For an ensemble of planetesimals, small objects are destroyed relatively more (less) often when q ≳ 2 (q ≲ 2). More destructive collisions among small objects reduce the slope; fewer destructive collisions raise the slope. Without any debris, this process leads to an equilibrium size distribution with q ≈ 2. In our simulations, debris has a power-law cumulative size distribution with q_debris = 2.5. Thus, cratering produces a power-law size distribution with q_cr ≈ 2.0–2.5, as observed in our calculations.

As the largest objects continue to grow, collisions among small objects become more and more catastrophic. For the fragmentation parameters we adopt, Q*_D monotonically increases with r for ≳ 0.1 km. In this range, small objects fragment catastrophically before large objects. Thus, catastrophic fragmentation reduces the equilibrium slope from the cratering slope of 2.0–2.5.

To provide a rough estimate of the equilibrium slope for catastrophic collisions, we consider a collision rate where the fraction of completely destructive collisions is a decreasing function of r, e.g., f_d∝r^−p. The loss rate of planetesimals is then

$$\begin{equation} d {\rm log} N / dt \approx f_d r^{2-q} \approx r^{2-p-q}. \end{equation} \tag{ A2 }$$

From the expression for Q*_D (Equation (6)), p ≈ 1.5. In the absence of debris production, catastrophic collisions lead to a power-law size distribution with q_cd ≈ 0.5, close to the q ≈ 0–1 derived in our simulations.