Inelastic collapse of perfectly inelastic particles

One of the most intensively discussed subjects in the dynamics of dissipative hard sphere systems is the effect of inelastic collapse, where the entire kinetic energy of the relative motion of a set of particles is dissipated in finite time due to an infinite sequence of collisions. The known collapse scenarios imply two preconditions: inertia of the particles and at least some degree of elasticity. For completely inelastic particles, collapse scenarios degenerate to a single sticky contact. By considering the overdamped motion of a frictional particle along the steepest descent in a rigid landscape, here we show that there exist collapse scenarios of novel type even if neither of these preconditions hold true. By means of numerical simulations we show that such collapses are no rare events due to particular particle shape and/or initial conditions and, thus, may be considered as an alternative scenario of granular cluster formation. Inelastic collapse, whereby particles can collide an infinite number of times in finite time, has been extensively studied, and occurs under the condition of particle inertia and elasticity. The authors propose a scenario for inelastic collapse in the physics of granular materials by numerically studying the deposition of totally inelastic particles in an overdamped dynamics showing that the two conditions above are not necessary to have inelastic collapse.

T he problem of inelastic collapse of identical, ideally hard, purely repulsive, and dissipatively colliding spheres reads as follows: given N particles in force-free space, are there initial conditions for the positions and velocities such that the entire energy of their relative motion will be dissipated? If dissipation is quantified by the coefficient of restitution, ε, relating the pre-collisional and post-collisional normal velocity, the question implies that there occurs an infinite sequence of collisions resulting in a cluster where each particle is in permanent contact with at least one neighbor. The answer to the question is "yes", provided the coefficient of restitution does not exceed a certain limit, ε ¼ 7 À 4 ffiffi ffi 3 p 1-3 . This result, obtained first for the case of three aligned particles attracted much attention and subsequently many related phenomena were studied, for instance the one-dimensional collapse of N particles 4,5 , frictional particles 6 , non-perfectly aligned particles 7 , collapses in two and threedimensional systems 8 , in the presence of solid walls [9][10][11] , under the influence of noise 5,[12][13][14][15] , and others. While collapse phenomena seem to be omnipresent for simple particles, the effect may be suppressed for more complex particles, for instance when the coefficient of restitution depends on the impact rate 16 or for particles having internal degrees of freedom 17 .
While the mathematical existence of inelastic collapse scenarios is free of doubt, the implications of a collapse for the physics of dissipative gases of hard spheres, frequently called granular gases, is still under debate, for example, see refs. [18][19][20] . An extended discussion on the collapse and its relation to granular clustering 2,21 , can be found in ref. 22 .
The infinite sequence of collisions leading to a collapse occurs in finite time 1,3 which lead to severe problems for event-driven molecular dynamics where time progresses in steps of collisions. There are several numerical tricks, which, however, fail in certain situations 23,24 or change the dynamics in an undesired way 25,26 . An event driven algorithm following the dynamics consistently must come to rest (that is, cannot propagate in time) due to the inelastic collapse 27 .
A close relative to the inelastic collapse revealing all its interesting features 20 is the problem of an ideally hard particle dropped from height h 0 to a plane where dissipation is again described by the coefficient of restitution, ε. The particle performs an infinite number of jumps of height h n+1 = ε 2 h n , (n = 0…∞) in finite time ffiffiffiffiffiffiffiffiffiffiffi ffi 2h 0 =g p ð1 þ εÞ=ð1 À εÞ, where g is gravity. Common to all collapse scenarios are two preconditions: (a) a certain degree of elasticity is needed and (b) the dynamics must be inertial. If violated, the sequence of collisions terminates after the first contact and the collapse scenario renders trivial.
In the present work, we show that for complex shaped, frictional particles moving in a rigid landscape with overdamped dynamics, there exist collapse scenarios, although both preconditions, (a) and (b), are not fulfilled. In the course of such a collapse, the particle experiences an infinite number of collisions occurring in finite time, in that, these scenarios resemble the collapse scenarios discussed above.

Results
Collapse in numerical simulations. To demonstrate the effect, particles as shown in Fig. 1 each consisting of four spheres are dropped from a point source to pile up to a heap due to the Visscher-Bolsterli algorithm 28,29 . Each dropped particle moves in the landscape formed by previously sedimented particles (see Fig. 1) until it finds its stable position in a local minimum. The dynamics is assumed overdamped, e.g., under the influence of an ambient viscous fluid 30 , corresponding to vanishing coefficient of restitution, ε = 0, irrespective of the particle material properties since the particles move always at terminal velocity, in the limit of large damping. Note that for the limit of large damping, longrange forces and lubrication forces vanish. No slip is assumed for the contact between the particle and the landscape. Consequently, due to overdamped dynamics where inertial forces are unimportant, the motion of the particle is solely determined by the instantaneous force, that is, by the gradient of the potential. Once a contact is established, the particle rolls such that its center of mass follows the steepest descent. Figure 1 shows the evolution of the vertical position of the particle's center of mass following the steepest descent. Time progresses in units of contact events, that is, a data point is drawn whenever a previously open contact between particles closes. At a certain time, t = 33, the particle enters a cycling motion, at three contact points in periodic sequence, and eventually experiences a collapse (dashed line) such that the simulation comes to rest. When entering the collapse scenario, the decrements, |Δz cm |, change qualitatively ( Fig. 1): For t ≥ 33, the decrement becomes regular; here the periodically active contact points are indicated by red, green and blue color. For t < 33 other contacts are active (black dots) and the decrement is on the scale of the particle size |Δz cm |~1.
The situation resembles the collapse scenario for granular gases where an infinite number of collisions (also occurring in cyclic sequence) takes place in finite time 1-3 , such that event-driven simulations cannot progress beyond this time. For event-driven Molecular Dynamics, the inelastic collapse is a severe problem, therefore, methods were developed [23][24][25][26] to abandon the infinite cyclic sequence by a small perturbation of the system which helps to overcome the collapse without noticeably changing the macroscopic system dynamics. Similarly, our system cannot progress beyond a collapse as one can see in Fig. 1 (dashed line). Therefore, in simulations we need a way to identify and overcome the collapse. A natural way to identify the contact of two objects is by the vanishing of the distance between their surfaces. In a collapse scenario, however, this condition is approached through an infinite sequence. Therefore, we need another condition to decide when particles touch one another. For the simulation in Fig. 1 (full line), we defined the contacts closed when the distance the center of mass traveled in one cycle drops below 10 −10 particle diameters. This condition is fulfilled at t = 66, where we assume all contacts of the The upper inset shows part of the heap's surface with the currently moving particle drawn in red. Blue and red lines: Vertical component of the particle's center of mass, z cm , over time in units of contact events. The dashed line shows the collapse scenario, the full line shows z cm with the regularization described in the text. Circles: decrement |Δz cm | of z cm in logscale. The color of the data points indicate the active contact. We see that for t ≥ 33 the motion is regular and cyclic particle closed and assess whether the particle found a stable position. This is not the case for the situation shown in Fig. 1, t = 66, such that the particle continues its motion along the steepest descent, that is, we overcame the collapse. Remarkably, a collapsed position does not necessarily correspond to the minimum of potential energy. Once all three contacts are closed, due to collapse, the particle can continue rolling on two contacts, following the steepest descent. The described collapse scenario could be observed for many particle shapes, except for particles whose surface is convex everywhere.
We wish to point out that the described model was successfully applied to a wide class of sedimentation phenomena with applications, ranging from geophysical systems 31,32 to nanostructured materials 29,33 . The numerical method was first described for spheres 28 and later generalized to more complex particles 29,34 .
Analytical characterization. We assume a particle shape as sketched in Fig. 2 which may look rather artificial but allows us to introduce the notation in a simple way. We will show that its motion along the steepest descent has the same properties as the motion of the particle in the sedimentation process presented above.
Consider the steepest descent of the particle in a landscape indicated by three patches, shown in Fig. 2. Figure 2a shows the initial state, with contacts C 1 and C 3 closed and C 2 open. The position is unstable, thus, the particle rolls in C 1 by opening C 3 until C 2 closes, Fig. 2b. Now the particle rolls in C 2 by opening C 1 until C 3 closes, Fig. 2c. To close the cycle, the particle rolls in C 3 by opening C 2 until C 1 closes. Now the particle is again in the situation sketched in Fig. 2a and the sequence repeats cyclically. Note that each step of this motion is related to lowering the center of mass, thus, each of the situations, Fig. 2a-c is unstable. The cycling motion terminates when simultaneously all three contacts are closed.
For quantitative analysis of the cyclic motion, we consider the evolution of the gap sizes g i for one cycle. Here, g i denotes the distance of the tip of the particle at contact i in the state when the contact is fully open, see Fig. 3. Consider one cycle, starting at the position shown in Fig. 2a. The ratio of initial gap size, g 2 , and the corresponding gap size in the successive cycle, g′ 2 reads In order to obtain g 3 /g 2 we look to the motion starting from the initial position, Fig. 2a, and the related rotation of the particle by the angle Δφ, see Fig. 3.
Assuming the gap size much smaller than the particle size, g 2 ≪ r 12 , we obtain in a linear approximation where r ij is the distance between the current point of rotation, C i , and the tip of the particle at the open contact, C j . α ij is the angle between the surface C j and the line connecting the rotation point, C i , with the point where the particle will contact the surface at C j . Figure 3 sketches the introduced variables. At the same time the tip of the particle at position C 3 is also rotated by Δφ around C 1 , therefore, g 3 r 13 cos α 13 j j ¼ Δφ ¼ g 2 r 12 cos α 12 j j ð3Þ The other quantities needed for Eq. (1), g 1 /g 3 and g′ 2 =g 1 , can be obtained in the same way by cyclic change of indexes. With r ij = r ji we obtain from Eq. (1):  Fig. 2 Contact patches of a complex particle and cyclic motion. A particle (blue) in overdamped motion in a rigid landscape specified by three patches (gray). Panel a shows the initial state, with contacts C 1 and C 3 closed and C 2 open. The position is unstable, thus, the particle rolls in C 1 by opening C 3 until C 2 closes, as sketched in panel b. Panel b: Now the particle rolls in C 2 by opening C 1 until C 3 closes, as sketched in panel c. Panel c: To close the cycle, the particle rolls in C 3 by opening C 2 until C 1 closes. Now the particle is again in the situation sketched in panel a. Red symbols indicate the particle's center of mass, its direction of motion and the current point of rotation. Black arrows show the direction of motion of the particle's edges contacting the landscape in cyclic sequence (a) → (b) → (c) → (a) → …  Fig. 3 Relation between g 2 and Δφ. Relation between g 2 (distance of the tip of the particle and contact point C 2 in the state when the contact is fully open) and Δφ (related rotation angle of the particle). α 12 is the angle between the surface C 2 and the line connecting the rotation point, C 1 , with the point where particle will contact the surface at C 2 . θ 1 is the angle between the vertical (dashed line) and the line connecting the center of mass with C 1 the gaps will reduce in geometric progression, thus, the series g → (g′ = δg) → (g″ = δ 2 g) → etc., converges to zero and the motion of the particle comes to rest after an infinite number of cycles when all contacts are closed, (g 1 , g 2 , g 3 ) → 0. It may be shown that δ given by Eq. (5) is always smaller than unity, consequently, if the particle performs one full cycle, it will necessarily enter a collapse scenario.
To compute the total distance traveled by the center of mass during the collapse scenario, again we consider the situation sketched in Fig. 2a. Since g 2 ≪ r 12 (and likewise for the other gaps), r cm moves along straight lines. The distance traveled during the first cycle, is Δr ð1Þ cm jΔr ð1Þ where subscripts denote the contact in which the particle rolls. From geometry we obtain and likewise Δ 2 r ð1Þ cm and Δ 3 r ð1Þ cm with cyclic permutation of the indexes where r icm is the distance between contact C i and the center of mass.
With Eq. (5), we obtain a recursion relation for the motion of the center of mass during the kth cycle, Δr ðkÞ cm ¼ δΔr cm . Summing up, we obtain the total path traveled by the center of mass: The distance traveled in vertical direction during the first cycle is where angles θ i are defined in Fig. 3. The total descent is, therefore, The analysis presented here applies to a wider class of systems: Indeed, the characteristic that the systems from Figs. 1 and 2 share that leads them to collapse is the identical relation for opening/closing contacts (up to permutation of (C 1 , C 2 , C 3 )): rolling on C 1 implies ð _ g 1 ¼ 0; _ g 2 <0; _ g 3 >0Þ, rolling on C 2 implies ð _ g 1 >0; _ g 2 ¼ 0; _ g 3 <0Þ, rolling on C 3 implies ð _ g 1 <0; _ g 2 >0; _ g 3 ¼ 0Þ, representing, thus, together with δ < 1 sufficient conditions for the collapse. For collapse scenarios involving three contacts, these conditions are also necessary -if any of them is violated, the sequence ceases. There may, however, exist more complex scenarios involving more than three contacts obeying different rules.
We wish to mention that there is a difference between the behavior of two and three-dimensional systems. Indeed, in three dimensions, a particle can undergo a collapse scenario and end up in a situation where all contacts are closed but the particle's position is not stable. This is, what we see in Fig. 1 at time t ≈ 66. In this case, after a collapse, the particle continues its descent in the landscape created by previously sedimented particles until it either finds its stable position or again undergoes a collapse after which it may be in a stable or unstable position. In contrast, in a two-dimensional system, a collapse scenario leads always to a stable position. Nevertheless, while the post-collapse behavior may be different in two-dimensional and three-dimensional systems, the mechanism of the collapse scenario described in this section is the same.
Numerical test. We assume the system sketched in Fig. 2a with the initial positions of the contact points and the center of mass C 1 = (0;10), C 2 = (5;5), C 3 = (0;0), and r cm = (5/3;5). The surface patches are inclined by 60°with respect to the horizontal. The initial position of the tip of the particle close to C 2 is (5.1;5.1). Starting from this initial configuration, in agreement with (8), the vertical component of the center of mass decays and saturates at a finite value while the decrements −Δz cm in each cycle decay exponentially, Fig. 4. Thus, after an infinite number of cycles, the particle assumes a position with all three contacts closed.
For the specified situation, from Eq. (5) we obtain an analytic expression for the recursion parameter, δ ≈ 0.066796, in good agreement with the fit shown in Fig. 4.
The described collapse scenario is not an exotic phenomenon occurring for a small set of particular initial conditions. To estimate the frequency of occurrence, we generated 10 9 sets of system configurations specified by the positions of C 1 , C 2 , C 3 , and r cm all chosen randomly from the interval [0, 1], and slopes of the surface patches from [−π/2, π/2]. Out of this set, approximately 3.6 × 10 6 systems led to a collapse. Given that a granular system may consist of many millions of complex shaped particles 29 , therefore, collapse scenarios are no rare events but will occur almost certainly in medium and large scale systems. Similar as for the known inelastic collapse in granular gases, a numerical simulation fails if the system enters a collapse scenario, since due to the infinite sequence of cycles, the algorithm would not progress in (real) time.

Discussion
Inelastic collapse scenarios of dissipative and purely repulsive hard-sphere systems, as extensively discussed in the literature so far, may be described as cyclic sequences of pairwise particle contacts. When a collapse occurs, an infinite sequence of collisions takes place in finite time such that eventually all energy of the relative velocity ceases. Perfectly inelastic particles interacting with coefficient of restitution, ε = 0, cannot undergo such an infinite sequence since the relative velocity of the particles would vanish already after the first contact. For the same reason, noninertial (overdamped) systems cannot cyclically collapse.
In the present work, we describe a novel collapse scenario occurring for a particle moving along the path of steepest descent of potential energy in a rigid landscape. The particle is assumed perfectly inelastic, ε = 0, and the dynamics is assumed overdamped, such that both preconditions of collapsing mentioned above do not hold true. For a model particle we could rigorously show that the collapse occurs via an infinite cyclic sequence of contacts where the decrement of the vertical component of the center of mass decreases geometrically such that the infinite sum of the decrements converges to a finite value. Our model assumes a rigid particle moving in a rigid landscape. For finite stiffness, in the first stage of the scenario, the gap size decays exponentially, since the mechanism driving the cycling motion is not affected by finite stiffness. Later, however, the sequence ceases, when the gap size approaches the typical deformation of particles in contact. Therefore, for finite stiffness, up to this criterion fulfilled, the arguments of our analysis remain invariant.
By extensive numerical simulations we showed that the collapse is not an exotic phenomenon but occurs for a non-negligible fraction (ca. 0.4%) of randomly generated initial conditions for our model particle. It was also shown that the phenomenon occurs for rather common systems, that is, for a particle consisting of five spheres toppling down the surface of a heap built from identical particles.
By now, the theoretically predicted phenomenon was not reported in the experimental literature. We believe that good candidates for experimental confirmation would be the analysis of the sound emission of a particle during sedimentation or the analysis of time series from electrical impedance measurements.

Data availability
All relevant data are available from the authors.