## Abstract

There is theoretical and observational evidence that asteroids and comets are conglomerations of particles ranging in size from dust grains to boulders. It is well known that energy added to such systems is dissipated by friction, plasticity and fracture. In addition to these physical phenomena, we find that energy can be dissipated in the form of particle kinetic energy due to random velocity distributions. ‘Dissipation’ in this manner is measured by what is called a granular temperature owing to its similarities with kinetic gas theory. This work has implications on our understanding of the growth of asteroids and comets.

## 1. Introduction

Observation of asteroid spins, asteroid densities, the break-up of comets and the theory of planetesimal growth imply that a considerable number of asteroids and comets are not cohesive masses of ice, rock and/or metal but instead a conglomeration of particles that range in size from micrometre-sized dust to kilometre-sized boulders (Richardson *et al*. 2002).

The purpose of this paper is to examine how these objects dissipate energy, i.e. the energy dissipation of rotating, self-gravitating granular media. This matter is relevant but not limited to the discussion of collisional coagulation during growth in the early Solar System (Weidenschilling 1997), the evolution of asteroids and comets spinning about non-principal axes of inertia (Burns & Safronov 1973) and situations of accretion by the gravitational collapse of a swarm of planetesimals (Tanga *et al*. 2004). Specifically, the aim is to investigate the role of a granular temperature in the energy dissipation of such systems. This is achieved by placing a pack of rigid, frictionless spheres in a non-equilibrium energy state using a discrete element method (Munjiza 2004). In the limit of many particles (with the properties just described), the results of these experiments will converge to that expected of an incompressible, inviscid fluid. Theoretical considerations of such fluids under centripetal, tidal and gravitational forces are well developed—a class of which are the Maclaurin spheroids.

## 2. Energy states of Maclaurin spheroids

Consider a self-gravitating ellipsoid, depicted in figure 1, with semi-axis lengths *a*_{1}, *a*_{2} and *a*_{3} (*a*_{1}≥*a*_{2}≥*a*_{3}) that is rotating about its maximum principal axis of inertia at non-dimensional angular velocity −*Λ* relative to a Cartesian coordinate system that is rotating at non-dimensional angular velocity *Ω* with respect to an inertial reference frame (Chandrasekhar 1969, §53, eqns (235) and (247))(2.1)where *λ* and *ω* are angular velocities; *G* is the universal gravitation constant; and *ρ* is the ellipsoid density.

If the ellipsoid is an incompressible and inviscid fluid and its shape is restricted to that of a spheroid (*a*_{1}=*a*_{2}≥*a*_{3}), then for a given value of non-dimensional angular momentum(2.2)(where *L* is the angular momentum and *M* is the mass of the spheroid), there exists an eccentricity(2.3)such that the body is in equilibrium.

These equilibrium configurations are known as the Maclaurin spheroids and are stationary points of the specific potential energy *W* subject to the constraints of constant mass and self-adjointness in the sense of Dedekind's theorem (Chandrasekhar 1969, §53, eqn (239) and §53.a)(2.4)

For a given value of , any perturbation in *e* away from equilibrium is equivalent to an increase in potential energy. Indefinite oscillations of the form(2.5)will result because the fluid is inviscid (zero energy dissipation). Oscillations of this type are known as spheroidal oscillations, the stability of which depends on the magnitude of the perturbation and whether the equilibrium configuration is a minimum, maximum or a saddle point (with respect to energy). It has been shown by these methods that small (infinitesimal) perturbations of the form of equation (2.5) are stable for all values of and that large (finite) perturbations of the form of equation (2.5) are stable as long as *q*<1, where *q* is the ratio of the difference in energy of the perturbed and equilibrium state to the energy at equilibrium (Chandrasekhar 1969; §53 eqn 285).

Table 1 shows examples of nine different non-equilibrium configurations (characterized by eccentricity *e*_{prtb}) partnered with the respective equilibrium configuration (characterized by eccentricity *e*_{equil}). Oscillations are stable for all cases because the modes are spheroidal and will continue indefinitely because the fluid is inviscid. To parallel the theoretical results using the discrete element method, a sufficient number of discrete elements would be needed so that collectively the system behaves as a continuum. Instead, in this work, the nine different cases are modelled using a system of *N*∼1000 particles and observations are made how the particle pack evolves from a state of non-minimal potential energy.

### (a) Case 1

We start with case 1 because it is the least complicated (zero spin) and easiest to visualize (the equilibrium shape is a sphere). This scenario is modelled using a collection of rigid, smooth and spherical particles arranged in a simple cubic structure. The only forces experienced by the particles are normal contact forces and self-gravity. By self-gravity we mean that every particle exerts a gravitational force on every other particle. We choose the radius of the spherical particles to be 100 m with a density of 1000 kg m^{−3}, from which it follows that the ellipsoid with *e*_{prtb}=0.8 containing the closest to *N*=1000 particles whose centres are within or on its boundary has dimensions *a*_{1}=*a*_{2}=1480 m and *a*_{3}=888 m. The version of a discrete element method employed is that detailed in the text by Munjiza (2004).

The value of *q* for this case is 0.022 and (as shown in table 1), so that one may expect indefinite oscillations (in time) about a sphere (the state of equilibrium). This is because all the particles are symmetrically aligned and there is no energy dissipation (e.g. friction, plastic deformation, fracture, etc.). However, the results of numerical experiments performed in this work demonstrate that this is not what happens; instead, the object collapses directly to the state of minimum energy. The only explanation for this is that energy is dissipated.

In figure 2, the sum of the kinetic energies and the gravitational potential energies of the individual particles are shown. At first glance, it appears that the plots of figure 2 are in conflict with the virial theorem because the system's kinetic energy does not approach one-half the absolute value of its potential energy. But this is only to be expected for the special case of self-gravity in the absence of all other forces. A more general form of the virial theorem, applicable to these simulations owing to the presence of contact forces, only requires that(2.6)where *m* is the particle mass; *v* is the particle speed; *x*, *y* and *z* denote the location of each particle in a Cartesian coordinate system; *X*, *Y* and *Z* are the components of the force acting on each particle resolved along the *x*, *y* and *z* directions, the summations are over all particles; and the overbars represent the quantity averaged over an infinite amount of time (Clausius 1870).

The potential energy of figure 2 is not dissipated in the traditional sense because the total energy, which is the sum of the potential energy and the kinetic energy, is constant. But what does happen is that the gravitational potential energy is converted into kinetic energy. This is not different than what happens for an inviscid, incompressible fluid continuum, but what is different is the form of the kinetic energy. It is ‘chaotic’ as evident from the small amplitude of its oscillations. An alternative measure of this type of kinetic energy is(2.7)which is simply the kinetic energy in terms of the thermal speed *c*′ averaged over all particles. For granular media, this is referred to as the *granular temperature* (e.g. Campbell 1990) because it is analogous to the temperature definition for a monatomic ideal gas(2.8)where *k* is the Boltzmann constant (Bird 1994). In figure 3, the system's granular temperature and the eccentricities both as a function of time are shown. Eccentricities were derived from the principal moments of inertia. It is clear from figure 2 that the system of particles has reached some state of minimum energy, but from figure 3*b* it is not clear if this state is the global minimum energy because the shape is not perfectly spherical. The value of averaged over frames 200–5500 (assumed to be steady state) is 0.1730 and that of =0.2411; in absolute terms, these correspond to *a*_{2}/*a*_{1}=0.985 and *a*_{3}/*a*_{1}=0.971, which is remarkably close to spherical, given the fact that the system comprises only 1000 spheres.

### (b) Cases 2–9

The time-averaged, steady-state eccentricities and the system granular temperature for the remaining cases are shown in table 2.1 An eccentricity of 0.31 corresponds to a 5% deviation from a sphere (*a*_{2}/*a*_{1}=0.95) and that of 0.44 to a 10% deviation from a sphere (*a*_{2}/*a*_{1}=0.90). By these criteria, we conclude that the discrete results start to significantly deviate from the continuum theory for case 6 .

This observation agrees well with the concept of granular temperature. A result stemming from kinetic gas theory is that the distribution of particle speeds is Maxwellian(2.9)(Bird 1994) with a mean *μ* and a variance *σ*^{2} of(2.10)(Hoel 1984); hence, high temperatures imply high particle velocities. These arguments are not directly applicable to the experiments at hand because kinetic gas theory assumes that the potential energy of the particle is negligible (Garrod 1995). However, the distribution of particle speeds at a particular instant of time is shown in figure 4 for cases 2 and 6 after steady state has been reached. The distributions are Maxwellian and show the same trend with temperature as that expected in an ideal gas.A more general form of Maxwell's distribution is(2.11)(Weisstein 2006), which can be fitted to the discrete data by equating the mean of the velocities and area under the curve; for case 2, *a*=0.0572, *c*=0.0152 and for case 6, *a*=0.1033, *c*=0.0152. The histogram bin size is equal to 0.0159 and was calculated by applying the Freedman–Diaconis rule (Izenman 1991) to data from case 6.

The data presented in figure 4 suggest that as the granular temperature rises, so does the probability that a particle will near or exceed its escape velocity. As this probability increases, the ellipsoid's ‘goodness of fit’ decreases because the particle velocities are large enough to launch them into orbit about the main pack. From this, we conclude that the discrete results start to deviate from the continuum results starting at case 6 because the granular temperature has reached a level where an ellipsoid is no longer a good fit to the shape of the main pack of particles. This is demonstrated in figure 5, which shows the particle locations at a specific instant of time as viewed along the axis of rotation (an ‘overhead’ view). At this point, it is appropriate to contrast the results of this work with that of Richardson *et al*. (2005). In the latter, a discrete element method was used to perform a series of numerical experiments with inelastic behaviour2 for comparison with the equilibria of incompressible, inviscid fluids. It was found that the presence of energy dissipation allowed a range of equilibrium shapes for a given level of the angular momentum. Configurations not initially in equilibrium were found to collapse to a lower energy state; sometimes retaining 100% of the mass and other times losing mass (as in figure 5). The concept of granular temperature provides an explanation of these observations.

## 3. Conclusions

We have shown using a discrete element method that the energy of a pack of self-gravitating and rotating rigid, smooth and spherical particles may be ‘dissipated’ by a rise in granular temperature.

Of course, real particles of the Solar System are not so idealized. Friction, fracture and plasticity are all present, which dissipate energy as heat. A more realistic picture of the process by which a granular system attains an energy minimum is as follows: potential energy is converted to ‘random’ or ‘chaotic’ particle motions, the kinetic energy of which is measured by a granular temperature. The granular temperature decreases as a function of time because the number of particle collisions increases, and during each collision kinetic energy is dissipated as heat. Hence, energy stored by a rise in granular temperature is only an interim state before dissipation by heat. The time scale on which a granular temperature is important is usually short, except for the case when there is a continuous input of energy that is greater than or equal to the amount of energy dissipated as heat (e.g. Campbell 1990). This view of the Solar System could help clarify the interpretation of experimental and observational data and provide insights into currently unsolved problems.

An example of such exists when considering the growth process of the early Solar System. Approximately 4.6 billion years ago, all masses were in the form of dust and gas that orbited the Sun in a large cloud called the solar nebula. The growth of kilometre-sized planetesimals from submicrometre-sized dust grains occurred by the collisions of successively more massive particles. Growth in this manner is known as a process of accretion (Weidenschilling 2000). Accretion of 1 cm–1 km bodies is not completely understood because relative impact velocities are too large for particles to ‘stick’ by atomic forces (e.g. van der Waals, hydrogen bonding, etc.) and are also larger than particle escape velocities. Mechanical interlocking of particles, ice coatings and magnetic and coulombic forces are hypothesized to enable growth in this regime (Weidenschilling 1997). However, a simpler explanation is that the impact energy is dissipated by an increase in granular temperature—yielding sufficiently low post-collisional velocities to allow capture by the growing body's gravity field.

## Acknowledgments

The authors would like to thank Prof. Pete Washabaugh of the University of Michigan for his thoughts and insight.