1 Introduction
The numerical resolution of the Nbody problem is one of the most popular problems in high performance computing. The intrinsic complexity of the problem has challenged many astronomers and computer scientists since the beginning of the computing era. Nbody algorithms have evolved side by side with the increase of computing power and availability of new resources, such as Graphic Processing Units (GPUs). Stateoftheart implementations include parallel Nbody tree codes (Morbidelli, 2002; Stadel, 2001; Richardson et al., 2002b), hybrid codes (Aarseth, 1999, 2001; Wang et al., 2015), adaptive algorithms of optimal orders (Pruett et al., 2003), systolic algorithms (Dorband et al., 2003), or more generally symplectic codes (Wisdom and Holman, 1991; Duncan et al., 1998; Chambers, 1999). In order to be able to solve granular dynamics problems, the capability of gravitational Nbody codes needs to be extended to include contact interaction between nonpointlike bodies. Two main classes of methods are commonly used to deal with contact interactions: hard and softbody methods. The first consider impulsive contacts between nondeformable particles (J. and E., 1959; Jean and Moreau, 1987). This method is numerically very stable and allows to simulate effectively the dynamics of hundreds of thousands of bodies. In more recent years, the alternative class of softbody methods have been developed (Cundall and Strack, 1979). Based on a forcedriven approach, softbody methods are more suitable than hardbody ones for the simulation of smooth dynamics. These could be used to complement some deficiencies of hardbody models. For example, due to the impulsive nature of contact exchanged, hardbody methods are not adequate to reproduce phenomena such as wave propagation in granular media (Gilardi and Sharf, 2002)
. However, the numerical solution of softbody methods is more unstable and often requires extremely small time steps to adequately reproduce elastic forces at contact. Both hard and softbody models have their own advantages and shortcomings and their use shall be carefully weighted upon the application scenario. As a general rule, hardbody methods qualify to simulate nonsmooth dynamics, and vice versa: softbody is suitable for smooth problems. However, whether a problem should be classified in either one category is a controversial field of debate: both models rely on sets of nondirectly measurable parameters, which often offer poor physical interpretation when compared to realworld scenarios
(Dubois et al., 2018). To date, both hard and softbody methods have been successfully used to simulate a wide variety of planetary science scenarios: Richardson et al. (2002a) used their hardbody model to simulate gravitational aggregation of rubblepile asteroids and planetary ring dynamics (Porco et al., 2008); Wada et al. (2006) used a softbody approach to study cratering impacts under constant gravity field and more recently, softsphere methods with mutual gravity between bodies have been used to simulate rubblepile scenarios (Sánchez and Scheeres, 2011; Tancredi et al., 2012; Schwartz et al., 2012). Another key aspect for a realistic evaluation of contact interactions concerns the shape of interacting bodies. All aforementioned codes are limited to sphericallyshaped bodies. The use of spherical shapes is now considered to be a very relevant simplification for problems where contact between bodies happens. The use of spheres allows to drastically reduce the computational burden, but can heavily affect the realism of the model when solving for the contact interaction between bodies (Michel et al., 2004). As such, the determination of the final shape of the aggregates and fragments could not be addressed (Richardson et al., 2009). Few attempts have been made using hardbody polyhedral grains together with video game engines to solve for contact interactions (Movshovitz et al., 2012). Although polyhedral bodies provide a better modeling, compared to spheres, of some physical phenomena (e.g. related to interlocking between bodies) the computational accuracy provided by game engines is very low (single precision) and not able to satisfy basic energy or angular momentum conservation requirements. Recently Ferrari et al. (2017) implemented the Nbody aggregation problem with nonspherical shapes and nonsmooth contacts, using modules of Chrono::Engine (C::E) (Tasora, 2019; Tasora et al., 2016), an open software optimized for granular and multibody engineering problems able to handle contacts and collisions of a large number of complexshaped bodies. The results obtained were very promising and satisfactorily accurate in terms of energy and angular momentum conservation. In their work, Ferrari et al. (2017) made use of direct Nbody simulations to solve for gravitational dynamics and only few thousands of bodies could be handled by the simulator.The work presented here builds upon the work by Ferrari et al. (2017) and discusses the implementation of a parallel CUDAGPU octree structure to extend the capabilities of code to a higher number of complexshaped bodies. As its parentcode, the implementation is based on a rework of the C::E multiphysics simulation engine to simulate collisions between nonspherical bodies and integrate the dynamics of the problem. The code features both nonsmooth (impulsebased) and smooth (forcebased) methods to handle contact interactions. Also, a compliant nonsmooth method is featured and proposed for the first time for planetary science application, in the attempt to joint the advantages of numerical stability provided by its impulsebased formulation, together with its ability of simulating compliance at contact level. Section 2 discusses the general architecture of the code and presents methods available to solve gravitational and collisional dynamics. The performance of the code are discussed in Section 3, where the results of validation/test problems are reported. The applicability to asteroid aggregation problems is eventually discussed through simulation examples. We address the problem of rubble pile asteroid aggregation processes in section 4. The physics of gravitational and collisional processes involved are reviewed and their numerical modeling is discussed. We derive a set of qualitative and quantitative requirements from dynamical and numerical constraints of the gravitationalcollisional problem.
2 Numerical implementation
The code is built using a modular approach, taking advantage of the objectoriented C++ architecture of C::E library. Modules represent interchangeable units of software, easily accessible by the user. Each module contains methods and routines to execute specific tasks. The main modules developed and/or directly retrieved by C::E libraries are:

gravity

contact

rigidbody dynamics

body creation

numerical solvers

data input/output

visual interface

post processing
Modules 2, 5, 7 are directly adapted from C::E libraries, modules 3, 4 are based on C::E libraries and further extended with additional methods and routines, modules 1, 6, 8 are fully developed by the authors. Modules 6 and 7 provide routines that allow the user to interface with the code at a higher level than software implementation level. The visual interface is based on Irrlicht library (Gebhardt, 2016) and can be set to highlight desired geometrical or simulation features (forces, contacts), as well as to display realtime simulation outputs. The data input/output interface is straightforward: every usertunable parameter of the model can be assigned via an input text file that is parsed and acquired by the code, while output data is typically saved into dedicated text files. This section presents a general overview of the numerical methods available in each module, focusing on those relevant to the gravitational aggregation problem.
2.1 Gravitational dynamics
Module 1 encloses methods and routines used to compute the gravity forces acting on each body. These are based on the evaluation of gravity field produced by pointmass or shapebased (e.g. polyhedron, ellipsoid) gravity sources. To simulate the problem of gravitational aggregation, the classical Nbody problem with pointmass sources is considered. This is a common and reasonable assumption: for a high number of commensurate interacting bodies, the use of shapebased mutual gravity would impact dramatically on the computational effort and would not provide relevant benefits to the accuracy of results, as discussed in section 4.1. As detailed below, module 1 includes both direct N integration and a parallelGPU octree method.
The Nbody problem
The equations of motion of the Nbody problem are typically written by using Newton’s law to compute the gravitational interactions between masses:
(1) 
where
represents the position vector of the
th body in an inertial frame, its mass, the universal gravitational constant and . A method implementing all NtoN interactions is commonly referred as direct N or particleparticle method (Hockney and Eastwood, 1988; Hut and Makino, 2010). Such method is very simple to implement and provides exact computation of the motion of the Nbodies without any physical approximation. The accuracy of the result depends solely on the numerical errors of the solver. The simple dynamics of direct N method make it very easy to implement and parallelize. However, they are extremely inefficient from a computational point of view: interactions must be computed for each body involved in the simulation and this leads to a time complexity of . A direct consequence is that such codes are limited in terms of maximum number of bodies to be handled. The bottlenecks are simulation time that becomes prohibitively long, and hardware capabilities with computer running out of memory, both effects occurring very rapidly as N increases. The direct N method, even with advanced parallelism on processing units, is not suitable for efficient asteroid aggregation simulations that can involve more than few thousands of bodies.Hierarchical treecode algorithm
Algorithms based on tree data structures rely on more dynamic and adaptive computations that allow for a significant reduction of time complexity up to . The BarnesHut (also referenced as BH in the followings) algorithm (Barnes and Hut, 1986) groups particles using a hierarchy of cube structures. A node in the algorithm corresponds to a cube in physical space. Because of the use of octrees, each node has eight child nodes obtained by a simple homogeneous spatial subdivision performed along the three principal axis of the system. The tree is therefore built by recursive subdivision until each node of the tree contains zero or one particle. The structure is adaptive, implying that the size of the tree is not fixed but comes as a result of the repartition of the particles in the 3D space. The data structure can grow naturally to more levels in regions where the particle density is high. The octree is rebuilt at each time step ^{1}^{1}1See section 4.1 for more details on the choice of simulation time step(s) of the dynamic simulation and consists of bodies stored at terminal nodes, called leaf nodes, and intermediary internal nodes that behave as clusters of particles. The algorithm uses the newly obtained adaptive octree data structure to compute the center of mass of the cubes involved in the simulation and the forces applied to each body: bodytobody and bodytonode interactions are considered.
Clustering
The algorithm relies on the idea that the force generated by a cluster of bodies can be approximated by treating the cluster as a single body. The accuracy of the approximations depends on the distance of the cluster from the body and the radius of the cluster of particles, as shown in Figure 1. We can therefore define the accuracy through , a usertunable parameter used to set the error of the method. The position of center of mass and total mass of the node is computed for each leaf and internal node of the BarnesHut tree. For a given internal node :
(2) 
where () are the coordinates of the position vector along the three principal axes of the system, while and (, , ) are respectively the mass and the coordinates of the position vector along the principal axes of the th child node of .
Force computation
After the octree is built and all nodes contain the required information, parameter is computed for each bodynode pair. (, ):
(3) 
This is compared with the value
chosen by the user to set the accuracy and performance of the force estimation algorithm. A more detailed discussion on the effects of the usertunable parameter
on both computational time and accuracy of the algorithm is provided in section 4.1. For each body , the tree is traversed from its root downwards along multiple branches and checks are performed at each encounter with a node , in order to compute force contributions:
If is a leaf node, the force contribution is evaluated as a classical bodytobody interaction:
(4) 
If is an internal node and , the traversal is interrupted at this node and the force contribution is evaluated:
(5) where r is the distance from the particle to the center of mass of the cube softened according to a parameter :
(6) is called softening length and, as its name suggests, ensures the “softening” of gravitational forces when , to avoid very high values of speed and acceleration that the integrator could not handle. Studies show that, without opting for a timedependent adaptive softening parameter, must be a factor of at least twice the minimum distance between a body and the attractive cluster involved (Rodionov and Sotnikova, 2005). In the case of finitesize body interaction, with contact and collisions, the softening parameter is not needed, since no overlap between center of mass of a body and a node could occur.

If is an internal node and , the traversal continues along all eight child nodes of .
CUDA implementation of GPUparallel BarnesHut
The implementation of the BarnesHut algorithm on a GPU using CUDA language is inspired by the work of Burtscher and Pingali (Burtscher and Pingali, 2011). The physical domain is divided into subdomains and the bodies are grouped following an octree structure. Similar to what done by Burtscher and Pingali (2011), the numerical tasks to follow the BarnesHut algorithm have been divided among five kernels, which are written to minimize memory accesses. The tasks are executed sequentially on the GPU as follows:

compute bounding box around all bodies (root of the octree)

build the octree

compute mass properties of each node

sort bodies by distance (to speed up force evaluation)

compute forces
Each body is assigned to a thread of the GPU, which traverses the octree to compute all force contributions acting on the body. Thanks to sorting, neighbor bodies require the traversal of a similar subpart of the octree, whereas distant bodies may require completely different traversals.
2.2 Rigidbody dynamics and contact interactions
Modules 3 and 4 deal with the creation of bodies and the definition of their inertial and surface properties. The latter depend on the contact method the user choose to use, among those available in 2. The following paragraphs recall briefly the main features of the code concerning body handling and contact interactions, focusing on the parameters relevant to each method. References are provided for a more complete discussion of each method.
Creation of bodies and their properties
The shape of each body can be selected by the user among common geometries (sphere, box, cone,…) or directly input by the user as a triangulated mesh. The shape can also be generated as the convex envelope of a randomly created cloud of points. When dealing with a large number of bodies, methods exist to select their spatial distribution or pattern (regular grid, randomly uniform, Poissondisk sampling,…) and bounding domain (common geometry or fitted in a triangulated mesh). Each body possesses 6 degrees of freedom, including both translational and rotational 3D motion. Random distribution routines can be used to set physical properties (size, density,…) and dynamical state (position, velocity, angular position, angular velocity) of bodies as well. Finally, the surface properties of each body can be defined. These are directly correlated to the contact method in use. Surface/contact parameters are discussed below both for smooth and for nonsmooth methods. More details can be found in the documentation of
Tasora (2019).Collision detection
Collision detection is performed into two steps: a broad and a narrow phase. During the broad phase, pairs of bodies whose geometrical boundaries are close enough are identified. The bounding volume of a body is estimated based on its velocity and integration time step. If the bounding volume of two bodies overlap, the narrow phase is performed and contact points are precisely found using a GJK algorithm (Tasora and Anitescu, 2010; Ferrari et al., 2017).
Contact method: smooth contacts
The Smooth Contact Method (SMC) implemented in C::E (Fleischmann et al., 2015)
is a forcebased softbody DEM method. The equations of motion are formulated as a system of Differential Algebraic Equations (DAE): namely a system of Ordinary Differential Equations (ODE) to reproduce the dynamics and Algebraic Equations (AE) for the kinematic constraints. Forces are exchanged between bodies at contact using a twoway normaltangent springdashpot system. The constitutive model of the springdashpot system can be selected between the simple Hooke and the more realistic Hertzian model. Relevant parameters of the SMC are coefficients of friction (static, dynamic, spinning), cohesion (as an attractive force at contact points), stiffness and damping. As any softbody DEM, it is suitable for smooth problems with no discontinuities but requires very small time steps for cases with high stiffness to avoid numerical instability.
Contact method: nonsmooth contacts
The NonSmooth Contact (NSC) method implemented in C::E (Tasora and Anitescu, 2010; Anitescu and Tasora, 2010; Tasora and Anitescu, 2011) is an impulsebased hardbody method. Unlike the SMC, the equations of motion are formulated as Differential Variational Inequalities (DVI) and require the solution of a Cone Complementarity Problem (CCP) at each time step. The NSC shares with the SMC parameters defining friction (static, dynamic, spinning) and adhesion. The handling of contact interaction is however much simpler and solely relies on the restitution coefficient, defined as the ratio between velocity after and before the collision. Due to its impulsemomentum formulation, the NSC is best suited for problems with discontinuities or with nearly rigid contacts (high stiffness).
Contact method: nonsmooth contacts with compliance
In the context of nonsmooth dynamics, a method is available to simulate contacts with compliance and damping (Tasora et al., 2013). As for the NSC method, the equations of motion are formulated as DVI and are based on a impulsemomentum formulation. However, compliance and damping are enforced at the constraint level and the method is suitable not only for nonsmooth problems but also to simulate the elastic and smooth behavior typical of softbody DEM models. Also, since it does not rely on the solution of a DAE, but its solution is found after solving a CCPbased DVI, this method does not require the very small time step required by softbody DEM. Its parameters are analogous to that of SMC: friction (static, dynamic, spinning), cohesion, stiffness and damping.
2.3 Numerical solvers
As discussed in the previous section, when using the SMC approach the dynamics are written as a system of DAE, and when using the NSC approach the dynamics are written as a DVI problem. From the numerical point of view, DAE require to numerically integrate ODEs, while DVI require the solution of an optimization problem at each time step. Despite the radical difference in approaching the solution, both approaches require a time stepper and a nonlinear solver. A large variety of time steppers and solvers is available in C::E (Tasora, 2017): the most interesting for gravitationalgranular dynamics applications include symplectic methods (semiimplicit Euler, leapfrog) and RungeKutta methods (RK45, explicit Euler, implicit Euler, trapezoidal, Heun), to be complemented with either iterative or direct solvers. In particular, symplectic methods are very suitable for gravitational dynamics problem and semiimplicit Euler method provide the best performance among them, when paired with an iterative solver (Anitescu and Tasora, 2010; Mangoni et al., 2018). This is suitable and effective to solve both DAE and DVI problems.
2.4 Postprocessing: aggregate identification
Module 8 contains routines for the postsimulation analysis of results. We provide here one example that applies to the case of rubblepile aggregation simulations. One of the tasks to be performed after such simulations is to identify the aggregate(s) of bodies and determine their properties such as shape, mass, inertia and void fraction (porosity).
Aggregate identification
The single or multiple aggregates obtained are identified through a method based on graph theory (Wilson, 1996). First, the graph is initialized to represent the whole system. Bodies are vertices of and the size of the vertex set is (number of bodies). Edges of represents physical contacts between objects: two bodies and , represented by vertices and of , are in contact if and only if there is an edge linking to . At the end of the simulation, a contact container is filled with pairs corresponding to pairs of bodies in contact with each other. The graph is filled with all existing contacts, drawn as edges in the graph. A schematic example is shown in Figure 2.
Aggregates obtained correspond to the connected subgraphs of . Single bodies, that do not belong to any aggregate (e.g because they have escaped from the reaggregation processes) are represented by isolated vertices. The number of aggregates is defined as the number of connected components of with a size greater than one, that is the number of the subset made of the connected subgraphs of with more than one vertex. A Depth First Search (DFS) algorithm (Wilson, 1996) is used to provide an effective way of traversing simple graphs. The output is a classification of the edges and a spanning tree that can be used for finding the connected components of a graph. Using such algorithm results in a complexity of (where and are respectively the number of vertices and edges of ) for the differentiation process.
Shape of the aggregate
The problem of defining the shape of a body made up of several discrete elements is not trivial and has a nonunique solution: infinite surfaces exist to envelope a given set of points. After identifying the aggregate and its children bodies, we address the problem of defining its overall shape by means of an alphashape algorithm (Edelsbrunner and Mücke, 1994). This computes the surface by enveloping a given set of points: in our case all vertices of children bodies. Unlike the simpler convex hull envelope, the alphashape surface is, in general, nonconvex. The working principle of the alphashape algorithm can be visually represented as follows. We start from a cloud of points (vertices of children bodies). Let us imagine to roll a sphere of radius over the cloud of points. The alphashape envelope is the surface created by the sphere as it rolls over the cloud of points. In the limiting case where , the surface enveloping the cloud would be a convex hull. For a finite value of , nonconvex envelopes can be obtained. Typically, should be on the same order of magnitude of the mean distance between vertices. As mentioned, the problem has a nonunique solution and the choice of makes the process arbitrary. This is affecting in a relevant manner the computation of global properties of the aggregate such as volume, void fraction (porosity) and shaperelated quantities (elongation, axis ratios, …). To solve for this arbitrariness, we define an admissible range of enveloping surfaces, within a range of limiting cases. In general, two limiting surfaces can be defined: (i) the minimum volume surface, which is found using the minimum value of such to have no surfaces delimiting cavities inside the aggregate, and (ii) the maximum volume surface, which is the convex hull of the aggregate, obtained with . Global properties are then computed considering the uncertainty on their values due to the arbitrary definition of the surface, and range between values computed using minimum and maximum volume surfaces.
3 Validation scenarios
Original modules of C::E dealing with the modeling of granular dynamics have been widely tested and validated against benchmark and experimental tests (see “Validation Studies section in Tasora (2019), Heyn (2013) and others). The results of further test simulations are shown here to validate the gravity modules and their integration into the overall implementation. In particular, the accuracy of direct N and octree gravity models are assessed and their performance is compared in terms of computational time. Conservation of energy, linear and angular momentum are checked for the purely gravitational case and for the coupled gravitycollision problem.
3.1 Direct Nbody integration
The first set of tests concerns the accuracy of force computation for the direct N method. Unlike the octree method, direct integration considers all gravitational interactions with no physical approximation. We performed several simulations of direct integration between mutually interacting bodies and in particular:

conic solutions of the twobody problem are reproduced as function of orbital energy;

orbits of planets in the Solar System are propagated for twothousand years;

Halo and Lyapunov orbits around collinear libration points of the SunEarth and EarthMoon circular restricted threebody problems (CR3BP) are reproduced.
When compared to the expected analytical (twobody problem) or numerical solution (Nbody problem and CR3BP), the code is able to reproduce exactly the result, where the term exactly (here and from now on) implies to within roundoff accuracy or numerical integration truncation errors. As mentioned in section 2.3, the semiimplicit Euler method is a very effective choice to solve gravitational dynamics, thanks to its symplectic nature that guarantees conservation of energy within the system. Energy, linear and angular momentum are conserved exactly for all cases, consistently with the nature of the direct N formulation that introduces no approximation errors in the model. The successful outcome of above tests guarantee that gravitational forces are computed correctly for each body. Figure 3 shows an example related to one of the aforementioned scenarios. The expected inplane wobbling of the Sun around the barycenter of the Solar System (data from JPL’s Horizons ephemerides, solid line) is compared and matches well with the outcome of a direct N integration of the motion of Sun and planets during a 100year timespan (cross markers).
3.2 GPUparallel octree
We evaluate here the performance of the GPUparallel implementation and compare it to direct N integration in terms of computational time and accuracy. Results presented here were obtained using a mediumrange laptop with

Intel(R) Core(TM) i76500U CPU @3.1Ghz

Nvidia Geforce 940M (384 CUDA cores).
Although certainly dependent on hardware architecture, the results here shown identify trends and behaviors that are not dependent on hardware, but only on the algorithms used. To assess the performance of algorithms, two different sets of simulations are performed. For each simulation, bodies are created at time zero in random positions within a cubic domain and with a uniform distribution. Collision detection and interaction is disabled during the entire duration of the simulation, such that bodies only experience mutual gravity. The results of validation scenarios are in agreement with previous assessments concerning tree codes
(Barnes and Hut, 1986; Hernquist, 1987).Computational time
A first set of simulations is arranged to compare the computational time required by each algorithm. Simulations are performed for several scenarios, including ranging from 10 up to 10, and using both direct N integration and BarnesHut GPU algorithm with different level of accuracy (). We compare the time required by the code to compute one integration time step. In principle, results can be compared after computing only one time step. However, to avoid any bias and for the sake of a more accurate profiling, we integrate the dynamics forward for several time steps and take the average time required to compute one of them. Also, we want the dynamics to be steady such to not interfere with the evaluation of performance. To this goal we use a very short time step (10 s) and a total duration of the simulation of 10 s. Results are shown in Figure 4 using a logarithmic scale. As expected, the cost of direct integration increases with N. As discussed in section 2.1, the BHGPU algorithm executes a sequence of kernels and the largest amount of computational time is shared between the creation of the octree and the evaluation of the gravitational bodybody or bodycluster interactions. Both octree and gravity evaluation depend on the number of bodies in the system, but to very different extents. Figure 4 shows clearly that for a low number of bodies, the computational time required is dominated by the octree creation phase, which is independent from the choice of the accuracy parameter and increases very slowly with N. Conversely, for high number of bodies the time expense is dominated by the evaluation of the gravitational interactions, which instead depends on . When comparing direct N with BHGPU, it clearly appears that for a low number of bodies (approximately N<700, the same for which is not affecting the time expense of BHGPU), the process of building the octree makes the BHGPU algorithm much more expensive than direct N integration. Conversely, for a high number of bodies, the BHGPU is much faster than direct N. For N=1,000 the BHGPU turns to be nearly two times faster and for N=10 it is faster of more than one order of magnitude. For a further increased number of bodies, direct N simulations become practically unfeasible and BHGPU becomes the only feasible option. As discussed, BHGPU algorithm can be tuned in accuracy: this affects the computational time as shown for different values used in Figure 4. In this highN regime, the computational time of BHGPU ranges from (with , i.e. no clustering of bodies, equivalent to direct integration) up to as increases.
Accuracy
A second set of simulations is used to assess the accuracy of the BHGPU algorithm. In this case, we wait for the bodies to complete a full orbit around the barycenter of the system and come back to their initial state at rest. Since collisions are disabled, all bodies are undisturbed in their gravitydriven path and the system exhibits a nearly perfect energyconserved pulsating contractingdilating behavior. This simulation time frame gives us a commensurate estimate of typical time duration of aggregation phenomena and provides us with every possible gravityrelated condition, i.e. both, when bodies are far away and when they are close to each other. Compared to the first set of simulations, we use a larger time step and we integrate forward the dynamics for a much longer time. Results are here shown for the case of N=1,000 and for . A simulation using direct N integration is also performed and used as a realworld result to assess the accuracy of the BHGPU algorithm. At each time step, the position of each body in the system is compared to its corresponding realworld one and the relative error is evaluated as percentage of the latter
(7) 
where and are the positions of body at time when using, respectively, BHGPU and N algorithms. Figure 5 shows the maximum position error among bodies in the system at each time step of the simulation. As expected, the error increases monotonically with , which is a direct mean and usertunable parameter to select the desired level of accuracy. Few more features are worth of discussion. The error is shown to be timedependent and exhibit a peak at half simulation time, i.e. at the end of contraction phase. This is expected, as the relative position error depends on the magnitude of accelerations acting on the bodies. The higher the acceleration, the higher the relative error when evaluating bodybody vs bodycluster interactions. Accordingly, the error is shown to be higher in the middle part of the simulation, where bodies are closer to each other and it increases during contraction and decreases as they get further to each other. It is interesting to look at the worst case depicted in Figure 5. At the end of contraction phase, the error between the two extremes and is of about two orders of magnitude. However, the maximum relative error obtained using is still quite small, in the order of 0.3%. Also, it is worth noting that this simulation provides us with an extremely conservative estimate of the accuracy. This is because realworld scenarios include the dissipative effects introduced by collisions and often entail only a part of the contraction phase simulated here. For example, if collisions were enabled in the scenario considered here, they would become very relevant starting at approximately 2.510 s (transient phase) and the bodies would never get as close to each other as they are in the [3.5, 6.5]10 s range due to their physical extent (aggregate phase). They will therefore experience lower gravitational acceleration actions and position errors due to BHGPU clustering: surely lower than the worstcase 0.3% value found here.
3.3 Collisions
A set of simulations is performed to test conservation laws during collisions and contact interactions between bodies. In this case, we simulated the settling dynamics of an already formed selfgravitating aggregate with 1,000 convexhull bodies. The dynamics are propagated forward for one time step and contacts are solved using both smooth and nonsmooth methods. The difference of physical and dynamical quantities are monitored. Several simulations are performed using the SMC method and different time steps, ranging between 10 and 10 seconds. The same is done using the NSC method which is able to handle higher time steps, between 10 to 10 seconds. A more detailed discussion on time step selection criteria is reported in Section 4.1. Unlike the purely gravitational case, energy, angular and linear momentum are not conserved exactly and collisions introduce approximation errors in the models. Since the semiimplicit Euler method is a first order integrator, the error is expected to depend linearly on the time step. This is confirmed by our simulations, as shown in Figure 6, which refer to the SMC case. It shows the angular momentum error in percentage, normalized to its initial value, for different time steps. As expected, a clear linear trend is identified for the error, as it scales with the time step.
4 The asteroid aggregation problem
This section discusses how rubblepile scenarios are simulated using the code described above (section 2). We discuss the appropriate choice of methods to use, how to tune their parameters and how to set up the numerical simulation.
4.1 Addressing the problem
To properly address the problem of simulating rubblepile dynamics, we study the fundamental interactions that occur between bodies during gravitational aggregation phenomena. These can be associated to phases of the gravitational aggregation process:

gravitational phase: the dynamics of the bodies is only driven by mutual gravity and they do not experience contact interactions.

transient phase: this phase begins as bodies start to collide with each other and form small aggregates of few bodies each. This phase is characterized by an extremely chaotic motion of the bodies that are simultaneously subjected to both gravity and random collisions. From a macroscopic point of view, the overall dynamics of the system is dominated by mutual gravity between bodies and newly forming aggregates. However, at particle level, the dynamics is dominated by collisions, which creates impulsive (nonsmooth contact) or quasiimpulsive (smooth) forces on a much faster scale than gravitational one. The combined effects at system and body level, and the dissipative nature of collisions leads eventually the system to a steady state condition.

aggregate phase: after transient phase, the system reaches a steady state equilibrium where bodies are either clustered in one or more aggregates, or dispersed.
The most challenging phase to be simulated is certainly the transient phase, where both gravitational and collision dynamics play a major role. Since gravity and collisions act at different scales, we explore the capability and suitability of the code by studying such dynamics separately. This leads to meaningful and general conclusions on how to numerically simulate such interactions, also applicable to the fully coupled problem.
Gravitational dynamics
In this early stage of the aggregation process, bodies are not in contact with each other and only interact gravitationally. As discussed in Section 3, the choice of the most suitable method depends on the number of bodies involved: direct Nbody for few bodies, octree for a large number of bodies with a hardwaredependent threshold typically between 1,000 and 10,000 bodies. In order to correctly simulate the Nbody dynamics, the integration time step must be consistent to the characteristic time of the problem and must sample the dynamics at least twice per . In particular, the following must be satisfied (Ferrari et al., 2017):
(8) 
where is the universal gravitational constant and is the material density of the bodies. For typical values of asteroid material density ranging between 1 and 4 g/cm (Richardson et al., 2002a), maximum time steps are in the order of thousands of seconds.
Collisional dynamics
As done for gravitational dynamics, we study here the fundamental features related to collisional interactions between bodies. We analyze the different phases and numerical tasks involved when solving for collisions and derive numerical constraints. From the numerical and algorithmic point of view, collisional dynamics requires the solution of two different problems: collision detection and collision output determination.
The collision detection algorithm operates at each time step to identify contacts between body pairs. To ensure proper collision detection, the time step must be adequate: if the time step is too high, a collision could be missed. To quantify the dynamical constraint on the time step, we consider here the case of a collision between two spheres of radius . In particular, we consider the limiting case when a collision is not detected. It is easy to assess that, for the case of a direct collision occurring along the line connecting the centers of the spheres, the distance traveled during the time step must be lower than 4, and then
(9) 
where is the relative velocity between bodies. Eq. (9) can be generalized for the case of grazing collisions (Sánchez and Scheeres, 2011):
(10) 
where is the maximum overlap allowed between bodies: the smaller the overlap allowed, the smaller the time step required and the more precise the detection algorithm. Note that the case of direct collision in (9) can be retrieved from (10) when (complete overlapping). When considering a system of bodies, is the maximum relative velocity in the system and is the sum of the radii of the two smallest particles. Also, although (10) refers to a simplified case with spheres, the same relation can be conservatively used with convex hulls of any shape, where represents their minimum characteristic size.
The second task is to solve for the actual collisional dynamics when bodies are in contact, in order to compute the collision output. This task is performed by the contact method selected, as discussed in section 2.2. Hard and softbody methods deal with a very different modeling of the physics at contact. Hardbody methods model instantaneous collisions and deal with nonsmooth dynamics based on impulsemomentum formulation. Since the collision is instantaneous, it does not make sense to speak of characteristic time involved in the dynamical process and then to derive a requirement on the time step from that. However, the time step still plays a major role and has to be selected properly to ensure the stability of the numerical solver. Hence, unlike gravity, collision detection and, as discussed below, softbody methods, the choice of the time steps for hardbody methods is not driven by dynamical constraints, but by the properties of the numerical solver. We consider here the case of the semiimplicit Euler method. Looking at its stability region (Niiranen, 1999), we can derive a constraint on the time step to be used, as function of the roots of the characteristic equation. However, our problem involves many bodies and requires the solution of a CCP at each time step: the roots of the characteristic equation of such problem are not of easy determination. Alternatively, an empirical but very effective method is to determine a proper time step by comparison against known benchmark problems or known behavior of the system. For example, when benchmark problems are not available, a very convenient way is to tune the time step according to a desired level of accuracy in terms of conservation of total angular momentum of the system . As discussed in Section 3, the error depends linearly on the time step. This linear dependency provides us with a very effective means to estimate and select properly the time step. Further examples on how the time step can be set are provided in C::E technical documentation and validation studies (Tasora, 2019), with reference to nonsmooth (NSC) and nonsmooth viscoelastic methods.
On the other hand, unlike hardbody ones, softbody methods directly rely on the physical modeling of collisional dynamics. The collision occurs in a finite time and foresee the viscoelastic behavior of the bodies at contact. Collisions are modeled by means of a springdashpot system at the contact point. We can therefore identify a characteristic time related to the fundamental frequency of the springdashpot system (Tsuji et al., 1993):
(11) 
When considering a system of bodies, is the highest frequency in the system, i.e. made of the combination of highest stiffness and lowest reduced mass ^{4}^{4}4 is the reduced mass of the two least massive bodies in the system. Accordingly, the time step must be at least two times smaller than the characteristic time
(12) 
Practical applications typically make use of smaller time steps, ranging from to (Sánchez and Scheeres, 2011), or even smaller than (Herrmann and Luding, 1998).
To summarize, two different constraints on the simulation time step can be derived after the analysis of the collisional process. The first one is expressed by Eq. (10). It concerns the capability of the code to properly detect collisions and does not depend on the contact method in use. The second one depends on the contact method in use. In case of impulsive contacts, it depends on the stability of the numerical solver and can be determined empirically. In case of smooth dynamics, it depends on the dynamics of contact as expressed in Eq. (12). Both requirements shall be satisfied to properly simulate collisional processes. This duality is easily resolved by enforcing the most stringent requirement, i.e. by choosing the lowest time step required. As a general rule, collision detection will be the bottleneck for high velocity impacts and vice versa, an accurate collision output will be the driving criteria in case of low velocity collisions. In the attempt to quantify in a more rigorous way this general criterion, we investigate possible links between Eq. (10) and (12). Indeed, while the two effects appear not to be directly coupled for hardbody methods, some considerations can be made for the case of softbody methods, where both (10) and (12) are built upon the dynamics of the problem. In particular, it is reasonable to assume that the maximum overlap allowed during collision detection must be consistently smaller than the maximum contraction of the springdashpot system. If this is not satisfied, the viscoelastic behavior at contact will not be accurately resolved. The maximum contraction of a springmass system can be computed based on energy conservation consideration, by equating the total energy right before contraction (only kinetic energy, no contraction) to the total energy at maximum contraction (only potential energy, zero velocity):
(13) 
from which is easy to find that
(14) 
In the case of a springdashpot system, the total energy is no longer conserved and a fraction of the initial kinetic energy is dissipated. The maximum contraction will be lower in this case, reduced by a factor <1
(15) 
Now we define the maximum overlap allowed during collision detection as a small fraction of :
(16) 
with , or equivalently
(17) 
with . We can now substitute (17) into (10) and get
(18) 
This relation correlates the two time steps involved in the numerical simulation of collision processes. In particular, we compare (18) with (12). The collision detection time step and the stiffness time step are equal if
(19) 
which can be rewritten as
(20) 
and finally, after trivial algebraic manipulation, as
(21) 
This provides us with the quantification of the general criterion stated above. In particular, when selecting the time step, collision detection will be the driving criterion when , i.e. when
(22) 
and vice versa. Eq. (22) gives us a useful condition to evaluate the tightest constraint through a simple comparison between maximum velocity and highest frequency in the system. Also, it correlates them with the accuracy of collision detection through . As expected, the problem of missing grazing collisions is an issue for high velocities and vice versa. More in detail, the break even point of the criterion is driven by the behavior , where is very small and is typically very high. Note that the dependency on shown in (22) might be misleading, since appears at the denominator of as well. In particular, the mass is proportional to and then the right hand side of Eq. (22) is proportional to .
The last point to discuss is about the choice of the contact method (SMC vs NSC vs NSC with compliance) to be used. As mentioned in section 2.2, SMC are naturally suited for smooth problems (low velocity contacts), NSC for nonsmooth problems (highvelocity impacts) while in principle, the NSC method with compliance can be used for both. A more detailed comparison between the aforementioned methods, including simulation of different scenarios, evaluation of performance and criteria for selection of parameters is reported in C::E technical documentation (Tasora, 2019). Studies to assess the accuracy of such methods are also reported. Unlike gravity, where the exact dynamics of the system are known and can be computed, contact methods do not model in a comprehensive fashion the physics involved in collisional processes. While the results of gravity algorithms can be validated using an analytical exact formulation, the only way to assess the accuracy of contact methods is through experimental test. As mentioned in Section 3, the documentation of C::E (Tasora, 2019) reports the results of several benchmark scenarios, including cone penetration test, standard triaxial test and direct shear test. Further analytical and experimental validation scenarios are reported in Heyn (2013).
The gravitycollision problem
The full gravitational aggregation problem features both gravitational and collisional dynamics and, in principle, both of them constraint the choice of the simulation time step. However, gravity and collisions act on very different time scales and the characteristic time of gravity is typically several orders of magnitude higher than that of collisions. Also, their dynamics are not coupled. For these reasons, it is not required to perform the simulation using only one time step (which would be the smallest one, i.e. the one related to collisions ). Instead, as done in Sánchez and Scheeres (2012), the aggregation scenario can be simulated using two time steps: the gravity time step , determined using (8) and the collision time step , determined using either (10) or (12). In particular, we advance the simulation of to properly solve collisional dynamics, but we evaluate gravity (through interactions or by building the octree) every . This allows for consistent savings in terms of computational time, without jeopardizing the accuracy of results.
In this context, a more accurate modeling of the local gravity interactions between fragments would not be beneficial towards a more realistic representation of the granular Nbody problem. Past astrodynamics studies have clearly shown that the mutual gravitational dynamics between objects of irregular shape, represented using shapebased (e.g. polyhedron) models, are very different from the case of pointmass sources interaction (e.g. Fahnestock and Scheeres, 2006; Ferrari and Lavagna, 2018). However, the granular Nbody problem is very different from that of astrodynamics applications, which consider purely gravitational motion with no contact and collision interactions. In our case, due to the presence of both gravity and contact interactions, acting at different time scales, a noticeable effect of deviations due to enhanced gravity modeling would appear after a time when many contact/collision interactions have already occurred. In this case, the chaotic nature of the contact/collision interactions represents a much higher source of uncertainty compared to the accuracy of local gravity models. Because of this, all authors identify the better modeling of contact/collision interactions as the way to enhance the realism of simulations (e.g. Michel et al., 2004; Richardson et al., 2009; Movshovitz et al., 2012, and others) and agree to using octrees, at the cost of a lower accuracy of local gravity computations. Unlike astrodynamics applications, the goal here is not to accurately reproduce the coupled and orbitalattitude dynamics of each fragments, but rather to reproduce the global behavior of the system. Relevant parameters concerning gravity are indeed global ones, such as bulk density and total mass of the system, and lowerfidelity local models of gravity (as in octrees) are widely used.
4.2 Features of the problem
As discussed, the study of the dynamics of selfgravitating aggregates is a very complex problem and deals with challenges arising from the coupled and chaotic interactions occurring within the granular media. The study of asteroid shapes and rubblepile dynamics has its fundamentals in the continuum theory, extended by Holsapple (2001, 2004, 2007, 2010) to the study of rubblepile objects. However, selfgravitating granular systems are not strengthless fluid bodies and, for example, can spin faster than a perfect fluid before shedding mass (Richardson et al., 2005). Even for cohesionless aggregates, many effects can contribute to providing strength to the aggregate. These include surface interaction phenomena, such as friction, and other effects related to the geometry of contact interactions between particles such as interlocking (Sánchez and Scheeres, 2012). To summarize all these effects, the strength of a rubble pile can be quantified based on its effective angle of friction, which results from particle shapes, size distributions, packing and surface friction (Sánchez and Scheeres, 2012). Holsapple (2010) showed that the angle of friction is very important to establish the domain of admissible equilibrium shapes and that spins at limiting conditions become larger for a higher angle of friction, implying a larger strength of the aggregate. Accordingly, Sánchez and Scheeres (2016) observe that the angle of friction inhibits deformation. In Zhang et al. (2017) geometric interlocking is the main source of shear strength for crystal structures and higher values can be achieved for random packing and increased size heterogeneity of fragments. Among geometrical effects, extreme relevance is given by the angularity of nonspherical particles. As mentioned, spheres behaves very differently from angular particles, for which the contact dynamic problem is very different. Due to their nonsmooth shape, angular particles can foster interlocking and between particles. In particular, spheres have only one contact point per pair, while nonspherical objects can have several. This has a great impact on all aspects related to contact dynamics and to the exchange of interactions/forces at contact points. For example, this affects how the value of parameters such as surface friction and restitution coefficients act on a global scale, since they both act at contact points. This means that their effect can be greatly amplified depending on relative geometry. In practice friction and restitution act on a larger scale and contact interactions results more dissipative with respect to the case of spheres. This is confirmed by experiments showing that restitution appears to be much lower for irregular objects (Hartmann, 1978) due to the loss of centerofmass kinetic energy to rotational energy postcollision (Korycansky and Asphaug, 2006). In practice, angular shapes can increase the strength (Richardson et al., 2005; Jiang et al., 2015; Zhang et al., 2018) of the aggregate and, through a lower restitution coefficient, increase the likelihood of the formation of a stable aggregate during a gravitational reaccumulation process (Walsh et al., 2008). The resolution of the model, in terms of number of bodies in the rubble pile aggregate, is also very relevant. Richardson et al. (2005) show that coarse configurations consisting of a small number of larger particles are more resistant to tidal disruption or reshaping than fine configurations with many smaller particles, observing that the strength of the aggregate depends on the number of particles as well. A major role in the granular dynamics processes is played by friction. Sánchez and Scheeres (2012) found that friction increases substantially the limiting spin rate. For the case of angular bodies, where the overall effect of friction is higher (many contact points) the expected limiting spin rate would be even higher. Sánchez and Scheeres (2012) also observe that the addition of surface friction makes the aggregate more stable. They conclude that friction might modify the square root dependency between critical spin of a granular aggregate and bulk density. The same is found by Zhang et al. (2017), who conclude that higher interparticle friction can keep a spinning rubble pile stable at a higher spin rate. Furthermore, Pravec and Harris (2007) show how maximum spin rate increases (substantially) with friction. To summarize, the chaotic dependency between the dynamical history, initial shape/configuration of the particles, contact/surface parameters and properties of the final aggregate (Zhang et al., 2017) gives a clear picture of the complexity of the problem and the unpredictability of its dynamics.
4.3 Simulation examples
A set of simulations is performed here to extend results presented by Ferrari et al. (2017) to an increased number of bodies and to verify the correlations between parameters of the simulation, initial conditions of the N bodies and aggregation outcome. To simulate realistic aggregation scenarios, it is important to carefully select the physical properties of the N bodies and their initial dynamics. Initial conditions play a crucial role to the formation of the aggregates and their properties. The initial dynamical state includes the position and velocity of the center of mass of each body, as well as the angular position and spin rate of each body. Typical asteroid aggregation scenarios consider these bodies as fragments created from a parent larger body after, for example, a collision event (Michel et al., 2001; Campo Bagatin et al., 2018). In such case, the fragments could share a common residual angular momentum due to spinning motion of the parent body. Parameter is introduced to simulate such effect, and it is defined as the initial residual angular velocity of the system of particles: when the fragments have a nonzero orbital velocity around the center of mass of the system. In order to compare and extend the results in Ferrari et al. (2017), similar scenarios are set up, using nonsmooth dynamics and a time step of 10 s. In our simulations, bodies are convex hulls whose shape is created as the geometrical envelope of randomly generated points. In particular, we use 16 points and corresponding convex hulls have, on average, 10 vertices. Surface friction is based on a simple Coulomb model with coefficient of 0.6. To foster the formation of aggregates, restitution is set to zero (perfectly inelastic collisions). Results presented here refer to simulations with 10 bodies, with average characteristic size of each fragment of 3 km and maximum initial distance between bodies in the order of 100 km. The physical properties of fragments are chosen among typical values of objects belonging to the main asteroid belt or Near Earth Asteroids population. The material density is set to 3 g/cm, a common choice for asteroid aggregation simulations (Sánchez and Scheeres, 2011), suitable to produce values of bulk density between 1.2 and 2.5 g/cm, typical of the asteroid population (Richardson et al., 2002a). At initial time, the angular and center of mass position of bodies are randomly generated with uniform spatial distribution, as well as the directions of their linear and angular velocity vectors.
An example of aggregation sequence is shown in Figure 7, for a simplified case with 10 bodies (to facilitate the visualization of images). While not in contact, the motion of the bodies is driven solely by their mutual attraction. When fragments start to interact at closer range, few small bodies are scattered away due to collision mechanisms. Depending on initial conditions, one, multiple or no stable aggregates are eventually formed. Aggregation is not observed when relative velocities between bodies are too high and in particular when parameters , and are above certain threshold values, above which no aggregation is possible. As for the case under study, these values are identified approximately with m/s, rad/s and rad/s. After reaching a steady state, aggregates are considered as single asteroids. Different kind of aggregates have resulted from the simulation performed.
Figure 8 shows the properties of the largest aggregate formed after the numerical simulation, as function of the initial dynamics of the system. In particular it shows (a) the inertial elongation
, defined as the ratio between maximum and minimum moment of inertia of the aggregate
(Ferrari et al., 2017), (b) the number of fragments in the aggregate and (c) its rotation period , as function of initial conditions (red diamonds), (blue asterisks) and (green stars). The values of are shown on the abscissa normalized to their maximum values, so that each parameter ranges from 0 to 1. Polynomial (for and ) and exponential (for ) fitting functions are used to highlight the trend of data distribution.The measure of the global properties and geometry of the final aggregate is affected by the uncertainty on the computation of its enveloping surface. As discussed in section 2.4, quantities such as volume, porosity, bulk density and size of the aggregate are better defined within a range of values computed between a minimum and a maximum volume surface. As a consequence, the inertial elongation and the angular velocity (or equivalently the rotation period ) of the aggregate are affected by the uncertainties of this measurement and should also be defined in within a range of values. In particular, the angular velocity of the aggregate is measured after computing its total angular momentum and inertia matrix , as , as in Richardson et al. (2005)
. The computation of the inertia tensor of the aggregate (and its principal inertia axes) rely on the definition of the enveloping surface. The uncertainty of this measurement depends on the size of the single particle with respect to the total aggregate. Hence, the accuracy depends on the number of particles in the aggregate and it is higher for smaller aggregates (those that spin faster in our simulations).
Richardson et al. (2005) estimate that the error in axis measurement could be as large as on particle radius along each axis, with errors of about 10% for the case of 1,000 bodies. When propagated to the computation of angular velocity, they estimate errors larger than 20%. Similarly, we estimated the error in the computed values of and , based on the number of bodies in the final aggregate, and show their uncertainties around the computed mean value in Figures 8 and 9. In our case the smaller aggregates are the fastest spinning ones: they are made of tens or few hundreds bodies, meaning potentially very large measurement errors. In other words, the resolution of our smaller aggregates is too coarse and does not allow to estimate accurately their spin rate, elongation, porosity and bulk density. Also, due to their low number of bodies, it is rather inaccurate to use the term rubble pile when referring to such small aggregates. To a closer inspection, the structure of these aggregates is closer to a monolith rather than a rubble pile: they are made of one (or few) big boulders and tens of smaller particles resting on their surface. This structure is by no means a typical selfgravitating aggregate and it’s rather closer to the “test particle on a rigid sphere” case reported in Weidenschilling (1981) and Richardson et al. (2000), where critical spin may be estimated using the simple balance between gravitational and centrifugal force (Hestroffer et al., 2019). This is indeed the result we find for them. The spin rate for aggregates with many more particles are comparable to what has been found by other authors within the error prescribed by the uncertainties mentioned above. Among these, our fastest spinning aggregates have spin rates in the order of 4.20.9 h.General trends of Figures 8 and 9 agree with expected behavior and are coherent with what discussed in section 4.2. As expected, higher elongations are obtained for higher values of , and . Maximum elongation is reached for values between 75%85% of their limiting value. For higher values of initial conditions (75%85%) the simulations result in the formation of multiple aggregates of smaller size, which are more typically regularly shaped than larger ones. It is indeed observed that maxima in Figure 8(a) represent thresholds after which the formation of more than one significant aggregate occurs. As confirmed in Figure 8(b), the higher the velocities, the lower number of bodies is found in the main aggregate, which is smaller for higher values of , and . Figure 8(c) confirms a wellknown result, with larger aggregates having a higher rotation period (lower spin rate) with respect to smaller ones. Figure 9 shows the rotation period of the final aggregate as function of (a) inertial elongation and (b) mass of the aggregate. As expected, fast rotators are more elongated and less massive, while slow rotators have a more regular rounded shape and are more massive.
5 Conclusion
The paper presents the numerical implementation of a code for coupled gravitygranular dynamics problems and reports numerical methods available to solve gravitational and collisional dynamics. We build upon the work by Ferrari et al. (2017) and extend the capabilities of the code by introducing a parallel CUDAGPU octree structure, to be able to evaluate mutual gravity for a higher number of bodies. Unlike existing Nbody codes, the implementation presented here has the capability to handle nonspherically shaped bodies to a full double precision accuracy. Compared to classical spherebased codes, this allows for a more realistic simulation of contacts between bodies and open to new opportunities and scenarios to be simulated. After presenting the results of test scenarios to validate newly implemented modules, we discuss the features and performance of numerical methods and derive requirements arising from the dynamics of asteroid aggregation scenarios. The problem of numerical simulation of gravitational aggregation is addressed: we provide guidelines and quantify criteria to properly set up such numerical simulations. We show examples of possible applications by performing numerical simulations of gravitational aggregation. The parameter space explored includes initial conditions of bodies in terms of relative velocity between them, spinning rate and residual angular velocity of the system. The outcome of the analysis is discussed showing the properties of the final selfgravitating aggregate after its stable formation in terms of inertial elongation, period of rotation and number of bodies in the aggregate.
Acknowledgements
This project has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie SkłodowskaCurie grant agreement No 800060. Part of the research work was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration. The authors would like to thank the anonymous reviewers for their comments and suggestions that helped to increase the quality of the paper.
References
 NBODY2: a direct nbody integration code. New Astronomy 6, pp. 277–291. Cited by: §1.
 From NBODY1 to NBODY6: the growth of an industry. 111 (765), pp. 1333–1346. External Links: Document, Link Cited by: §1.
 An iterative approach for cone complementarity problems for nonsmooth dynamics. 47 (2), pp. 207–235. Cited by: §2.2, §2.3.
 A hierarchical o(n log n) forcecalculation algorithm. Nature 324, pp. 446–449. External Links: Link Cited by: §2.1, §3.2.
 GPU computing gems. pp. 75 – 92. Cited by: §2.1.
 Internal structure of asteroid gravitational aggregates. 302, pp. 343 – 359. External Links: ISSN 00191035, Document, Link Cited by: §4.3.
 A hybrid symplectic integrator that permits close encounters between massive bodies. Monthly Notices of the Royal Astronomical Society 304, pp. 793–799. Cited by: §1.
 A discrete numerical model for granular assemblies. 29 (1), pp. 47–65. External Links: Document, https://doi.org/10.1680/geot.1979.29.1.47, Link Cited by: §1.
 Systolic and hypersystolic algorithms for the gravitational nbody problem, with an application to brownian motion. Journal of Computational Physics 185, pp. 484–511. Cited by: §1.
 The contact dynamics method: a nonsmooth story. Comptes Rendus Mécanique 346 (3), pp. 247 – 262. Note: The legacy of JeanJacques Moreau in mechanics / L’héritage de JeanJacques Moreau en mécanique External Links: ISSN 16310721, Document, Link Cited by: §1.
 A multiple time step symplectic algorithm for integrating close encounters. Astronomical Journal 116, pp. 2067–2077. Cited by: §1.
 Threedimensional alpha shapes. Transactions on Graphics 13 (1), pp. 43–72. Cited by: §2.4.
 Simulation of the full two rigid body problem using polyhedral mutual potential and potential derivatives approach. 96 (3), pp. 317–339. External Links: ISSN 15729478, Document, Link Cited by: §4.1.
 Nbody gravitational and contact dynamics for asteroid aggregation. Multibody System Dynamics 39 (1), pp. 3–20. External Links: ISSN 1573272X, Document, Link Cited by: §1, §1, §2.2, §4.1, §4.3, §4.3, §5.
 Ballistic landing design on binary asteroids: the aim case study. 62 (8), pp. 2245 – 2260. Note: Past, Present and Future of Small Body Science and Exploration External Links: ISSN 02731177, Document, Link Cited by: §4.1.
 On the importance of displacement history in softbody contact models. 11 (4), pp. 044502–044502–5. External Links: ISSN 15551415, Document, Link Cited by: §2.2.
 IrrLicht release 1.8.4. Note: http://irrlicht.sourceforge.netLast accessed 20161201 Cited by: §2.
 Literature survey of contact dynamics modelling. Mechanism and Machine Theory 37 (10), pp. 1213 – 1239. External Links: ISSN 0094114X, Document, Link Cited by: §1.
 Planet formation: mechanism of early growth. 33 (1), pp. 50 – 61. External Links: ISSN 00191035, Document, Link Cited by: §4.2.
 Performance characteristics of tree codes. 64, pp. 715–734. External Links: Document Cited by: §3.2.
 Modeling granular media on the computer. 10 (4), pp. 189–231. External Links: ISSN 14320959, Document, Link Cited by: §4.1.
 Small solar system bodies as granular media. 27 (1), pp. 6. External Links: ISSN 14320754, Document, Link Cited by: §4.3.
 On the modeling, simulation, and visualization of manybody dynamics problems with friction and contact. Ph.D. dissertation, University of WisconsinMadison. Cited by: §3, §4.1.
 Computer simulation using particles. CRC Press. Cited by: §2.1.
 Equilibrium configurations of solid cohesionless bodies. 154 (2), pp. 432 – 448. External Links: ISSN 00191035, Document, Link Cited by: §4.2.
 Equilibrium figures of spinning bodies with selfgravity. 172 (1), pp. 272 – 303. Note: Special Issue: CassiniHuygens at Jupiter External Links: ISSN 00191035, Document, Link Cited by: §4.2.
 Spin limits of solar system bodies: from the small fastrotators to 2003 el61. 187 (2), pp. 500 – 509. External Links: ISSN 00191035, Document, Link Cited by: §4.2.
 On yorpinduced spin deformations of asteroids. 205 (2), pp. 430 – 442. External Links: ISSN 00191035, Document, Link Cited by: §4.2.
 Moving stars around. Cited by: §2.1.
 Studies in molecular dynamics. i. general method. 31 (2), pp. 459–466. External Links: Document, https://doi.org/10.1063/1.1730376, Link Cited by: §1.
 Dynamics in the presence of unilateral contacts and dry friction: a numerical approach. In Unilateral Problems in Structural Analysis — 2, G. Del Piero and F. Maceri (Eds.), Vienna, pp. 151–196. External Links: ISBN 9783709129678 Cited by: §1.
 A novel threedimensional contact model for granulates incorporating rolling and twisting resistances. 65, pp. 147 – 163. External Links: ISSN 0266352X, Document, Link Cited by: §4.2.
 Lowspeed impacts between rubble piles modeled as collections of polyhedra. 181 (2), pp. 605 – 617. External Links: ISSN 00191035, Document, Link Cited by: §4.2.
 A primal–dual predictor–corrector interior point method for nonsmooth contact dynamics. 330, pp. 351 – 367. External Links: ISSN 00457825, Document, Link Cited by: §2.3.
 Catastrophic disruption of asteroids and family formation: a review of numerical simulations including both fragmentation and gravitational reaccumulations. Planetary and Space Science 52, pp. 1109–1117. Cited by: §1, §4.1.
 Collisions and gravitational reaccumulation: forming asteroid families and satellites. Science 294 (5547), pp. 1696–1700. Cited by: §4.3.
 Modern integrations of solar system dynamics. Annual Reviews of Earth and Planetary Sciences 30 (1), pp. 89–112. Cited by: §1.
 Numerical modeling of the disruption of comet d/1993 f2 shoemakerlevy 9 representing the progenitor by a gravitationally bound assemblage of randomly shaped polyhedra. 759 (2), pp. 93. External Links: Link Cited by: §1, §4.1.
 Cited by: §4.1.
 Simulations of the dynamical and lightscattering behavior of saturn’s rings and the derivation of ring particle and disk properties. 136 (5), pp. 2172. External Links: Link Cited by: §1.
 Binary asteroid population: 1. angular momentum content. 190 (1), pp. 250 – 259. External Links: ISSN 00191035, Document, Link Cited by: §4.2.
 An adaptive nbody algorithm of optimal order. Journal of Computational Physics 187, pp. 298–317. Cited by: §1.
 Numerical simulations of asteroids modelled as gravitational aggregates. Planetary and Space Science 57, pp. 183–192. Cited by: §1, §4.1.
 Gravitational aggregates: evidence and evolution. In W.F. Bottke Jr., A. Cellino, P. Paolicchi, R.P Binzel (Ed) Asteroids III, Tucson (Ed.), Vol. 160, pp. 501–515. Cited by: §1, §4.1, §4.3.
 Direct largescale Nbody simulations of planetesimal dynamics. Planetary and Space Science 30 (1), pp. 45 – 59. Cited by: §1.
 Numerical experiments with rubble piles: equilibrium shapes and spins. 173 (2), pp. 349 – 361. External Links: ISSN 00191035, Document, Link Cited by: §4.2, §4.3.
 Direct largescale nbody simulations of planetesimal dynamics. 143 (1), pp. 45 – 59. External Links: ISSN 00191035, Document, Link Cited by: §4.3.
 Optimal choice of the softening length and time step in nbody simulations. Astronomy Reports 49 (6), pp. 470–476. External Links: ISSN 15626881, Document, Link Cited by: 2nd item.
 DEM simulation of rotationinduced reshaping and disruption of rubblepile asteroids. 218 (2), pp. 876 – 894. External Links: ISSN 00191035, Document, Link Cited by: §4.1, §4.2.
 Simulating asteroid rubble piles with a selfgravitating softsphere distinct element method model. 727 (2), pp. 120. External Links: Link Cited by: §1, §4.1, §4.1, §4.3.
 Disruption patterns of rotating selfgravitating aggregates: a survey on angle of friction and tensile strength. 271, pp. 453 – 471. External Links: ISSN 00191035, Document, Link Cited by: §4.2.
 An implementation of the softsphere discrete element method in a highperformance parallel gravity treecode. 14 (3), pp. 363–380. External Links: ISSN 14347636, Document, Link Cited by: §1.
 Cosmological Nbody simulations and their analysis. Ph.D. Thesis, University of Washington, Seattle. Cited by: §1.
 Granular physics in lowgravity environments using discrete element method. Monthly Notices of the Royal Astronomical Society 420 (4), pp. 3368–3380. External Links: Document, /oup/backfile/content_public/journal/mnras/420/4/10.1111_j.13652966.2011.20259.x/2/mnras04203368.pdf, Link Cited by: §1.
 A matrixfree cone complementarity approach for solving largescale, nonsmooth, rigid body dynamics. 200, pp. 439–453. Cited by: §2.2.
 A compliant viscoplastic particle contact model based on differential variational inequalities. 53, pp. 2 – 12. Note: Multibody System Dynamics: A Selective Perspective External Links: ISSN 00207462, Document, Link Cited by: §2.2.
 Chrono::engine project. Note: http://projectchrono.org/ Cited by: §1, §2.2, §3, §4.1, §4.1.
 External Links: Link Cited by: §2.3.
 A convex complementarity approach for simulating large granular flows. 5 (3), pp. 031004–031004–10. External Links: ISSN 15551415, Document, Link Cited by: §2.2, §2.2.

Chrono: an open source multiphysics dynamics engine
. In High Performance Computing in Science and Engineering, T. Kozubek, R. Blaheta, J. Šístek, M. Rozložník, and M. Čermák (Eds.), Cham, pp. 19–49. External Links: ISBN 9783319403618 Cited by: §1.  Discrete particle simulation of twodimensional fluidized bed. 77 (1), pp. 79 – 87. External Links: ISSN 00325910, Document, Link Cited by: §4.1.
 Numerical simulation of impact cratering on granular material. Icarus 180 (2), pp. 528 – 545. External Links: ISSN 00191035, Document, Link Cited by: §1.
 Rotational breakup as the origin of small binary asteroids. 454, pp. 188. External Links: Link Cited by: §4.2.
 nbody6++gpu: ready for the gravitational millionbody problem. 450 (4), pp. 4070–4080. External Links: ISSN 00358711, Document, http://oup.prod.sis.lan/mnras/articlepdf/450/4/4070/5782448/stv817.pdf, Link Cited by: §1.
 How fast can an asteroid spin?. 46 (1), pp. 124 – 126. External Links: ISSN 00191035, Document, Link Cited by: §4.3.
 Introduction to graph theory. Prentiss Hall. Cited by: §2.4, §2.4.
 Symplectic maps for the nbody problem. Astronomical Journal 102, pp. 1528–1538. Cited by: §1.
 Creep stability of the proposed aida mission target 65803 didymos: i. discrete cohesionless granular physics model. 294, pp. 98 – 123. External Links: ISSN 00191035, Document, Link Cited by: §4.2.
 Rotational failure of rubblepile bodies: influences of shear and cohesive strengths. 857 (1), pp. 15. External Links: Document, Link Cited by: §4.2.
Comments
There are no comments yet.