The Soft-Sphere Discrete Element Method in pkdgrav3 for Small-Body Simulations

Originally developed for darkmatter-only, cosmological N-body simulations, pkdgrav is a massively-parallel N -body code that uses a k-D binary tree for domain decomposition and for approximating the gravity of distant particles [1]. A treatment for collisions was added which enabled applications in planetary science [2], such as the first full simulations of asteroid family formation by asteroid disruption and gravitational reaccumulation [3]. In the planetary science community, pkdgrav is also known for its implementation of the soft-sphere discrete element method (SSDEM) which is a technique used in granular mechanics for modeling the frictional forces between contacting particles [4,5]. In this model, interpenetrating particles feel a restoring normal force based on a Hooke’s law spring, as well as rolling, twisting, and sliding frictional forces based on user-defined coefficients that can be set to represent realistic granular media. Over the last∼25 years, pkdgrav has been been a cutting-edge toolin the planetary sciences, and used to make significant advancements in formation of planetesimals and planets [2,6,7], asteroid families [3], and asteroid satellites [8-10], as well as the study of Saturn’s rings [11] and asteroid surface processes and seismology [12,13]. It has also been used to support small body space missions including the OSIRIS-REx, New Horizons, DART, and Hera missions. Owing to its k−D tree, pkdgrav’s runtime scales as O(N logN), where N is the number of particles in the simulation, which is a significant improvement from a simple brute-force N -body code which has O(N2) scaling. However, the original pkdgrav architecture is now more than two decades old and increasingly limited relative to modern computational approaches. We present a new implementation of SSDEM in pkdgrav3. pkdgrav3 is a massively parallel, GPU-accelerated N-body and hydrodynamics code designed for a range of astrophysical problems [14]. Building on the legacy of the original pkdgrav framework, pkdgrav3 combines a hierarchical tree algorithm with a Fast Multipole Method (FMM) gravity solver that reduces the runtime scaling from O(Nlog⁡N) to O(N), enabling simulations with billions to trillions of self-gravitating particles. The code is optimized for modern heterogeneous supercomputing architectures through hybrid MPI/pthreads parallelization, asynchronous communication, SIMD vectorization, and GPU acceleration, allowing efficient scaling across thousands of CPU cores and GPUs. Its scalability, computational efficiency, and low memory footprint have enabled some of the largest cosmological simulations ever performed, including simulations with four trillion self-gravitating particles for the Euclid mission’s flagship mock galaxy catalog [15]. Smoothed Particle Hydrodynamics (SPH) has recently been implemented within the code, enabling shock physics simulations on planetary and asteroid size scales [16]. We will describe the numerical algorithms underlying the pkdgrav3 SSDEM implementation and present benchmarking results demonstrating the code’s performance. We will report the current status of the code’s development and discuss features planned for future releases. We will also discuss some planned applications for this code, including high-resolution simulations of asteroid family formation, planetesimal formation, and asteroid tidal encounters in preparation for Apophis’ close Earth encounter. References:[1] J. Stadel, “Cosmological N-body simulations and their analysis,” University of Washington, 2001.[2] D. Richardson, “Direct Large-Scale N-Body Simulations of Planetesimal Dynamics,” Icarus, vol. 143, no. 1, pp. 45–59, Jan. 2000, doi: 10.1006/icar.1999.6243.[3] P. Michel, W. Benz, P. Tanga, and D. C. Richardson, “Collisions and Gravitational Reaccumulation: Forming Asteroid Families and Satellites,” Science, vol. 294, no. 5547, pp. 1696–1700, Nov. 2001, doi: 10.1126/science.1065189.[4] S. R. Schwartz, D. C. Richardson, and P. Michel, “An implementation of the soft-sphere discrete element method in a high-performance parallel gravity tree-code,” Granular Matter, vol. 14, no. 3, pp. 363–380, Mar. 2012, doi: 10.1007/s10035-012-0346-z.[5] Y. Zhang et al., “Creep stability of the proposed AIDA mission target 65803 Didymos: I. Discrete cohesionless granular physics model,” Icarus, vol. 294, pp. 98–123, Sep. 2017, doi: 10.1016/j.icarus.2017.04.027.[6] D. Nesvorný, A. N. Youdin, and D. C. Richardson, “FORMATION OF KUIPER BELT BINARIES BY GRAVITATIONAL COLLAPSE,” The Astronomical Journal, vol. 140, no. 3, pp. 785–793, Aug. 2010, doi: 10.1088/0004-6256/140/3/785.[7] J. C. Marohnic et al., “Constraining the final merger of contact binary (486958) Arrokoth with soft-sphere discrete element simulations,” Icarus, vol. 356, p. 113824, Mar. 2021, doi: 10.1016/j.icarus.2020.113824.[8] K. J. Walsh, D. C. Richardson, and P. Michel, “Rotational breakup as the origin of small binary asteroids,” Nature, vol. 454, no. 7201, pp. 188–191, Jul. 2008, doi: 10.1038/nature07078.[9] K. J. Walsh, R.-L. Ballouz, H. F. Agrusa, J. Hanus̆, M. Jutzi, and P. Michel, “Satellite Formation around the Largest Asteroids,” The Astrophysical Journal Letters, vol. 986, no. 1, p. L12, Jun. 2025, doi: 10.3847/2041-8213/adc562.[10] H. F. Agrusa et al., “Direct N-body Simulations of Satellite Formation around Small Asteroids: Insights from DART’s Encounter with the Didymos System,” The Planetary Science Journal, vol. 5, no. 2, p. 54, Feb. 2024, doi: 10.3847/psj/ad206b.[11] R.-L. Ballouz, D. C. Richardson, and R. Morishima, “Numerical Simulations of Saturn’s B Ring: Granular Friction as a Mediator between Self-gravity Wakes and Viscous Overstability,” The Astronomical Journal, vol. 153, no. 4, p. 146, Mar. 2017, doi: 10.3847/1538-3881/aa60be.[12] Y. Kim, J. V. DeMartini, D. C. Richardson, and M. Hirabayashi, “Tidal resurfacing model for (99942) Apophis during the 2029 close approach with Earth,” Monthly Notices of the Royal Astronomical Society, vol. 520, no. 3, pp. 3405–3415, Feb. 2023, doi: 10.1093/mnras/stad351.[13] J. V. DeMartini et al., “Using a discrete element method to investigate seismic response and spin change of 99942 Apophis during its 2029 tidal encounter with Earth,” Icarus, vol. 328, pp. 93–103, Aug. 2019, doi: 10.1016/j.icarus.2019.03.015.[14] D. Potter, J. Stadel, and R. Teyssier, “PKDGRAV3: beyond trillion particle cosmological simulations for the next era of galaxy surveys,” Computational Astrophysics and Cosmology, vol. 4, no. 1, May 2017, doi: 10.1186/s40668-017-0021-1.[15] F. J. Castander et al., “Euclid,” Astronomy & Astrophysics, vol. 697, p. A5, Apr. 2025, doi: 10.1051/0004-6361/202450853.[16] T. Meier, D. Potter, C. Reinhardt, and J. Stadel, “Smoothed Particle Hydrodynamics in pkdgrav3 for Shock Physics Simulations. I. Hydrodynamics,” The Astrophysical Journal, vol. 1000, no. 2, p. 266, Mar. 2026, doi: 10.3847/1538-4357/ae4e29.