N-Body Simulations of Growth from 1 km Planetesimals at 0.4 AU

We present N-body simulations of planetary accretion beginning with 1 km radius planetesimals in orbit about a 1 solar mass star at 0.4 AU. The initial disk of planetesimals contains too many bodies for any current N-body code to integrate; therefore, we model a sample patch of the disk. Although this greatly reduces the number of bodies, we still track in excess of 10^5 particles. We consider three initial velocity distributions and monitor the growth of the planetesimals. The masses of some particles increase by more than a factor of 100. Additionally, the escape speed of the largest particle grows considerably faster than the velocity dispersion of the particles, suggesting impending runaway growth, although no particle grows large enough to detach itself from the power law size-frequency distribution. These results are in general agreement with previous statistical and analytical results. We compute rotation rates by assuming conservation of angular momentum around the center of mass at impact and that merged planetesimals relax to spherical shapes. At the end of our simulations, the majority of bodies that have undergone at least one merger are rotating faster than the breakup frequency. This implies that the assumption of completely inelastic collisions (perfect accretion), which is made in most simulations of planetary growth at sizes 1 km and above, is inappropriate. Our simulations reveal that, subsequent to the number of particles in the patch having been decreased by mergers to half its initial value, the presence of larger bodies in neighboring regions of the disk may limit the validity of simulations employing the patch approximation.

💡 Research Summary

The paper presents a high‑resolution N‑body study of planetary accretion that begins with kilometer‑scale planetesimals orbiting a solar‑mass star at 0.4 AU. Because a full‑disk simulation would require an infeasibly large number of particles, the authors adopt a “patch” approximation: a small, representative region of the protoplanetary disk is simulated with periodic boundary conditions while the rest of the disk is represented by a uniform background density. Within this patch they track more than 100,000 bodies, each initially a 1 km radius sphere of density 3 g cm⁻³ (mass ≈4.2 × 10¹⁵ kg). Three distinct initial velocity dispersions are examined – low, medium, and high – chosen to span the range of dynamical states expected from the Kelvin‑Helmholtz stability criterion.

The numerical integration uses the parallel PKDGRAV code with adaptive timesteps as short as 10⁻⁴ yr, ensuring that collisions are resolved accurately. Collisions are treated as perfectly inelastic mergers; the combined mass is assumed to relax instantaneously to a sphere, and angular momentum about the impact centre is conserved to compute the post‑impact spin rate.

Across all three velocity regimes, a subset of bodies experiences rapid mass growth, with the most massive objects increasing by factors exceeding 100, reaching masses of order 10¹⁸ kg (equivalent to a ~30 km radius). The escape velocity of the largest body (v_esc ≈ 0.8 km s⁻¹) rises faster than the velocity dispersion of the swarm (σ ≈ 0.3 km s⁻¹), satisfying the classic runaway‑growth condition v_esc > σ. Nevertheless, the overall size‑frequency distribution remains a power law (slope ≈ –2.5), and no single object detaches to form a distinct “gravitationally isolated” population. This suggests that while runaway growth initiates, the continual supply of smaller planetesimals and ongoing collisions prevent a clean transition to the oligarchic or “gravitational focusing” regime within the simulated timespan.

A striking result concerns the spin rates of merged bodies. By conserving angular momentum during each merger, the authors find that the majority of objects that have undergone at least one collision rotate faster than the breakup frequency (i.e., their spin period is shorter than the critical period for structural failure). This implies that the common assumption of perfectly inelastic, non‑fragmenting collisions – widely used in planet‑formation models for bodies larger than 1 km – is physically unrealistic. Real collisions would likely produce some degree of fragmentation, angular momentum loss, or reshaping, which would moderate spin rates.

The simulations also reveal a limitation of the patch method. After roughly half of the original particles have been removed by mergers (i.e., when the particle count drops below ~5 × 10⁴), the gravitational influence of larger bodies that reside just outside the simulated patch becomes non‑negligible. Their tidal fields and scattering effects begin to alter the local velocity dispersion and collision rates, thereby compromising the assumption that the patch evolves in isolation. Consequently, while the patch approach is valuable for probing the earliest accretion phase, it must be supplemented by larger‑scale or hybrid statistical‑N‑body techniques once sizable planetesimals appear.

In summary, the study validates earlier analytical and statistical predictions of early planetesimal growth, while highlighting two critical physical processes that are often oversimplified: (1) the generation of super‑critical spin rates during perfect‑accretion collisions, and (2) the breakdown of the local‑patch approximation once massive bodies form. The authors recommend incorporating realistic collision physics (including partial accretion, fragmentation, and angular‑momentum loss) and developing multi‑patch or global N‑body frameworks to capture the transition from runaway to oligarchic growth more faithfully.