Size distribution of particles in Saturn’s rings from aggregation and fragmentation

Brilliantov, Nikolai; Krapivsky, P. L.; Bodrova, Anna; Spahn, Frank; Hayakawa, Hisao; Stadnichuk, Vladimir; Schmidt, Jürgen

doi:10.1073/pnas.1503957112

Earth, Atmospheric, and Planetary Sciences

Size distribution of particles in Saturn’s rings from aggregation and fragmentation

Nikolai Brilliantov, P. L. Krapivsky, Anna Bodrova, +3 , Frank Spahn, Hisao Hayakawa, Vladimir Stadnichuk, and Jürgen Schmidt [email protected]Authors Info & Affiliations

Edited by Neta A. Bahcall, Princeton University, Princeton, NJ, and approved June 12, 2015 (received for review February 26, 2015)

July 16, 2015

112 (31) 9536-9541

https://doi.org/10.1073/pnas.1503957112

PDF/EPUB

Significance

Although it is well accepted that the particle size distribution in Saturn’s rings is not primordial, it remains unclear whether the observed distribution is unique or universal, that is, whether it is determined by the history of the rings and details of the particle interaction or whether the distribution is generic for all planetary rings. We show that a power-law size distribution with large-size cutoff, as observed in Saturn’s rings, is universal for systems where a balance between aggregation and disruptive collisions is steadily sustained. Hence, the same size distribution is expected for any ring system where collisions play a role, like the Uranian rings, the recently discovered rings of Chariklo and Chiron, and possibly rings around extrasolar objects.

Abstract

Saturn’s rings consist of a huge number of water ice particles, with a tiny addition of rocky material. They form a flat disk, as the result of an interplay of angular momentum conservation and the steady loss of energy in dissipative interparticle collisions. For particles in the size range from a few centimeters to a few meters, a power-law distribution of radii,

\sim r^{- q}

with

q \approx 3

, has been inferred; for larger sizes, the distribution has a steep cutoff. It has been suggested that this size distribution may arise from a balance between aggregation and fragmentation of ring particles, yet neither the power-law dependence nor the upper size cutoff have been established on theoretical grounds. Here we propose a model for the particle size distribution that quantitatively explains the observations. In accordance with data, our model predicts the exponent q to be constrained to the interval

2.75 \leq q \leq 3.5

. Also an exponential cutoff for larger particle sizes establishes naturally with the cutoff radius being set by the relative frequency of aggregating and disruptive collisions. This cutoff is much smaller than the typical scale of microstructures seen in Saturn’s rings.

Bombardment of Saturn’s rings by interplanetary meteoroids (1–3) and the observation of rapid processes in the ring system (4) indicate that the shape of the particle size distribution is likely not primordial or a direct result of the ring-creating event. Rather, ring particles are believed to be involved in active accretion–destruction dynamics (5–13) and their sizes vary over a few orders of magnitude as a power law (14–17), with a sharp cutoff for large sizes (18–21). Moreover, tidal forces fail to explain the abrupt decay of the size distribution for house-sized particles (22). One wishes to understand the following: (i) Can the interplay between aggregation and fragmentation lead to the observed size distribution? And (ii) is this distribution peculiar for Saturn’s rings, or is it universal for planetary rings in general? To answer these questions quantitatively, one needs to elaborate a detailed model of the kinetic processes in which the ring particles are involved. Here we develop a theory that quantitatively explains the observed properties of the particle size distribution and show that these properties are generic for a steady state, when a balance between aggregation and fragmentation holds. Our model is based on the hypothesis that collisions are binary and that they may be classified as aggregative, restitutive, or disruptive (including collisions with erosion); which type of collision is realized depends on the relative speed of colliding particles and their masses. We apply the kinetic theory of granular gases (23, 24) for the multicomponent system of ring particles to quantify the collision rate and the type of collision.

Results and Discussion

Model.

Ring particles may be treated as aggregates built up of primary grains (9) of a certain size

r_{1}

and mass

m_{1}

. [Observations indicate that particles below a certain radius are absent in dense rings (16).] Denote by

m_{k} = k m_{1}

the mass of ring particles of “size” k containing k primary grains and by

n_{k}

their number density. For the purpose of a kinetic description we assume that all particles are spheres; then the radius of an agglomerate containing k monomers is

r_{k} = r_{1} k^{1 / 3}

. (In principle, aggregates can be fractal objects, so that

r_{k} \sim k^{1 / D}

, where

D \leq 3

is the fractal dimension of aggregates. For dense planetary rings it is reasonable to assume that aggregates are compact, so

D = 3

.) Systems composed of spherical particles may be described in the framework of the Enskog–Boltzmann theory (25–27). In this case the rate of binary collisions depends on particle dimension and relative velocity. The cross-section for the collision of particles of size i and j can be written as

σ_{i j}^{2} = {(r_{i} + r_{j})}^{2} = r_{1}^{2} {(i^{1 / 3} + j^{1 / 3})}^{2}

. The relative speed [on the order of

0.01 - 0.1 cm / s

(16)] is determined by the velocity dispersions

〈 v_{i}^{2} 〉

and

〈 v_{j}^{2} 〉

for particles of size i and j. The velocity dispersion quantifies the root-mean-square deviation of particle velocities from the orbital speed (

\sim 20 km / s

). These deviations follow a certain distribution, implying a range of interparticle impact speeds and, thus, different collisional outcomes. The detailed analysis of an impact shows that for collisions at small relative velocities, when the relative kinetic energy is smaller than a certain threshold energy,

E_{agg}

, particles stick together to form a joint aggregate (11, 28, 29). This occurs because adhesive forces acting between ice particles’ surfaces are strong enough to keep them together. For larger velocities, particles rebound with a partial loss of their kinetic energy. For still larger impact speeds, the relative kinetic energy exceeds the threshold energy for fragmentation,

E_{frag}

, and particles break into pieces (29).

Using kinetic theory of granular gases one can find the collision frequency for all kinds of collisions and the respective rate coefficients:

C_{i j}

for collisions leading to merging and

A_{i j}

for disruptive collisions. The coefficients

C_{i j}

give the number of aggregates of size

(i + j)

forming per unit time in a unit volume as a result of aggregative collisions involving particles of size i and j. Similarly,

A_{i j}

quantify disruptive collisions, when particles of size i and j collide and break into smaller pieces. These rate coefficients depend on masses of particles, velocity dispersions, and threshold energies,

E_{agg}

and

E_{frag}

:

C_{i j} = ν_{i j} (1 - (1 + B_{i j} E_{agg}) \exp (- B_{i j} E_{agg})) A_{i j} = ν_{i j} \exp (- B_{i j} E_{frag}) ν_{i j} = 4 σ_{i j}^{2} \sqrt{\frac{π (〈 v_{i}^{2} 〉 + 〈 v_{j}^{2} 〉)}{6}} B_{i j} = 3 \frac{m_{i}^{- 1} + m_{j}^{- 1}}{〈 v_{i}^{2} 〉 + 〈 v_{j}^{2} 〉} .

[1]

These results follow from the Boltzmann equation, which describes evolution of a system in terms of the joint size–velocity distribution function (section below and SI Text). The governing rate equations for the concentrations of particles of size k read

\frac{d n_{k}}{d t} = \frac{1}{2} \sum_{i + j = k} C_{i j} n_{i} n_{j} - \sum_{i = 1}^{\infty} C_{k i} n_{i} n_{k} - \sum_{i = 1}^{\infty} A_{k i} n_{i} n_{k} (1 - δ_{k 1}) + \sum_{i = 1}^{k} n_{i} \sum_{j = k + 1}^{\infty} A_{i j} n_{j} x_{k} (j) + \frac{1}{2} \sum_{i, j \geq k + 1} A_{i j} n_{i} n_{j} [x_{k} (i) + x_{k} (j)] .

[2]

The first term on the right-hand side of Eq. 2 describes the rate at which aggregates of size k are formed in aggregative collisions of particles i and j (the factor

\frac{1}{2}

avoids double counting). The second and third terms give the rates at which the particles of size k disappear in collisions with other particles of any size i, due to aggregation and fragmentation, respectively. The fourth and fifth terms account for production of particles of size k due to disruption of larger particles. Here

x_{k} (i)

is the total number of debris of size k, produced in the disruption of a projectile of size i. We have analyzed two models for the distribution of debris

x_{k} (i)

. One is the complete fragmentation model,

x_{k} (i) = i δ_{1 k}

, when both colliding particles disintegrate into monomers; another is a power-law fragmentation model, when the distribution of debris sizes obeys a power law,

x_{k} (i) \sim B (i) k^{- α}

, in agreement with experimental observations (30, 31); the impact of collisions with erosion is also analyzed.

Decomposition into Monomers.

In the case of complete fragmentation,

x_{k} (i) = i δ_{1 k}

, the general kinetic equations [2] become

\frac{d n_{k}}{d t} = \frac{1}{2} \sum_{i + j = k} C_{i j} n_{i} n_{j} - \sum_{i \geq 1} (C_{i k} + A_{i k}) n_{i} n_{k}

[3]

\begin{matrix} \frac{d n_{1}}{d t} = - n_{1} \sum_{j \geq 1} C_{1 j} n_{j} + n_{1} \sum_{j \geq 2} j A_{1 j} n_{j} \\ + \frac{1}{2} \sum_{i, j \geq 2} A_{i j} (i + j) n_{i} n_{j} . \end{matrix}

[4]

Mathematically similar equations modeling a physically different setting (e.g., fragmentation was assumed to be spontaneous and collisional) have been analyzed in the context of rain drop formation (32).

Constant rate coefficients.

The case of constant

C_{i j} = C_{0}

and

A_{i j} = λ C_{0}

can be treated analytically, providing useful insight into the general structure of solutions of Eqs. 3 and 4, explicitly showing the emergence of the steady state. The constant λ here characterizes the relative frequency of disruptive and aggregative collisions. Without loss of generality we set

C_{0} = 1

. Solving the governing equations for monodisperse initial conditions,

n_{k} (t = 0) = δ_{k, 1}

, one finds

n_{1} (t) = λ_{1} [1 + λ^{- 1} {(λ_{2}^{- 1} e^{λ t} - \frac{λ^{- 1}}{2})}^{- λ_{2} / λ_{1}}],

[5]

where

λ_{1} = λ / (1 + λ)

and

λ_{2} = 2 λ / (1 + 2 λ)

. Using the recursive nature of Eq. 3, one can determine

n_{k} (t)

for

k > 1

. The system demonstrates a relaxation behavior: After a relaxation time that scales as

λ^{- 1}

, the system approaches a steady state with

n_{1} = λ_{1}

, the other concentrations satisfying

0 = \frac{1}{2} \sum_{i + j = k} n_{i} n_{j} - (1 + λ) n_{k} N .

[6]

Here

N = λ_{2}

is the steady-state value of the total number density of aggregates,

N = \sum_{k \geq 1} n_{k}

. We solve [6] using the generation function technique to yield

n_{k} = \frac{N}{\sqrt{4 π}} (1 + λ) {[\frac{2 n_{1}}{(1 + λ) N}]}^{k} \frac{Γ (k - 1 / 2)}{Γ (k + 1)} .

[7]

Now we assume that disruptive collisions in rings are considerably less frequent than aggregative ones, so that

λ ≪ 1

(this assumption leads to results that are consistent with observations); moreover,

k ≫ 1

for most of the ring particles. Using the steady-state values

n_{1} = λ_{1}

and

N = λ_{2}

, one can rewrite Eq. 7 for

k ≫ 1

as

n_{k} = \frac{λ}{\sqrt{π}} e^{- λ^{2} k} k^{- 3 / 2} .

[8]

Thus, for

k < λ^{- 2}

, the mass distribution exhibits power-law behavior,

n_{k} \sim k^{- 3 / 2}

, with an exponential cutoff for larger k.

Size-dependent rate coefficients.

For a more realistic description, one must take into account the dependence of the rate coefficients on the aggregate size (Eq. 1). Here we present the results for two basic limiting cases that reflect the most prominent features of the system:

i)

The first case corresponds to energy equipartition,

〈 E_{k} 〉 = (1 / 2) m_{k} 〈 v_{k}^{2} 〉 = const

, which implies that the energy of random motion is equally distributed among all species, like in molecular gases. In systems of dissipatively colliding particles, like planetary rings, this is usually not fulfilled, the smaller particles being colder than suggested by equipartition (33, 34). We also assume that the threshold energies of aggregation and fragmentation are constant:

E_{agg} = const

and

E_{frag} = const

; the latter quantities may be regarded as effective average values for all collisions. Then, as follows from Eq. 1, we have

λ = A_{i j} / C_{i j} = const

, and the kinetic coefficients read

C_{i j} = C_{0} {(i^{1 / 3} + j^{1 / 3})}^{2} {(i^{- 1} + j^{- 1})}^{1 / 2},

[9]

where

C_{0} = const

, so that the

C_{i j}

are homogeneous functions of the masses of colliding particles

C_{a i, a j} = a^{ϰ} C_{i j} .

[10]

The specific form [9] implies that the homogeneity degree is

ϰ = 1 / 6

.

ii)

The second limiting case corresponds to equal velocity dispersion for all species,

〈 v_{i}^{2} 〉 = v_{0}^{2} = const

. In planetary rings the smaller particles do have larger velocity dispersions than the larger ones but they are by far not as hot as equipartition would imply (33). Thus, this limiting case of equal velocity dispersions is closer to the situation in the rings. For the dependence of the fragmentation threshold energy

E_{frag}

on the masses of colliding aggregates we use the symmetric function

E_{frag} = E_{0} (i j / (i + j))

, which implies that

E_{frag}

is proportional to the reduced mass of the colliding pair,

μ_{i j} = m_{1} (i j / (i + j))

. This yields

B_{i j} E_{frag} = (3 E_{0} / 2 m_{1} v_{0}^{2}) = const

and allows a simplified analysis. We assume that the aggregation threshold energy

E_{agg}

for all colliding pairs is large compared with the average kinetic energy of the relative motion of colliding pairs,

(1 / 2) μ_{i j} v_{0}^{2}

(our detailed analysis confirms this assumption) (SI Text). Then

\exp (- B_{i j} E_{agg}) ≪ 1

and Eq. 1 yields

C_{i j} = ν_{i j}

. Therefore, the ratio

A_{i j} / C_{i j} = \exp (- B_{i j} E_{frag})

is again constant. Thus, the relative frequency of disruptive and aggregative collisions is also characterized by the constant

λ = A_{i j} / C_{i j}

. The kinetic coefficients attain now the form

C_{i j} = {\tilde{C}}_{0} {(i^{1 / 3} + j^{1 / 3})}^{2},

[11]

which is again a homogeneous function of i and j but with different homogeneity degree

ϰ = 2 / 3

.

An important property of the kinetic equations, where the rate coefficients

C_{i j}

and

A_{i j} = λ C_{i j}

are homogeneous functions of i and j, is that these equations possess a scaling solution for

i, j ≫ 1

. The latter is determined by the homogeneity degree ϰ and is practically insensitive to the detailed form of these coefficients (35, 36). We use this property and replace the original rate coefficients [9] and [11] by the generalized product kernel

C_{i j} = {\hat{C}}_{0} {(i j)}^{μ}, {\hat{C}}_{0} = const .

[12]

For this kernel, the homogeneity degree is

ϰ = 2 μ

. To match it with the homogeneity degree of [9] and [11] we choose

μ = 1 / 12

for the first limiting case and

μ = 1 / 3

for the second. The advantage of the product kernel [12] is the existence of an analytic solution for the steady-state distribution. Indeed, with the homogeneous coefficients [12] the steady-state version of Eq. 3 reads

0 = \frac{1}{2} \sum_{i + j = k} l_{i} l_{j} - (1 + λ) l_{k} L,

[13]

where we have used the shorthand notations

l_{k} = k^{μ} n_{k}, L = \sum_{k \geq 1} l_{k} .

[14]

With the substitute,

n_{k} \to l_{k}

and

N \to L

, the system of equations [13] is mathematically identical to the system of equations with a constant kernel [6], so that the steady-state solution reads

n_{k} = \frac{L}{2 \sqrt{π}} e^{- λ^{2} k} k^{- 3 / 2 - μ},

[15]

again a power-law dependence with exponential cutoff.

Our analytical findings are confirmed by simulations. In Fig. 1, the results of a direct numerical solution of the system of rate Eqs. 3 and 4 are shown for both limiting kernels [9] and [11], together with their simplified counterparts [12]. The stationary distributions for the systems with the complete kinetic coefficients [9] and [10] have exactly the same slope as the systems with the simplified kernel [12] of the same degree of homogeneity and hence quantitatively agree with the theoretical prediction [15]. Moreover, the numerical solutions demonstrate an exponential cutoff for large k, in agreement with the theoretical predictions.

Fig. 1.

Particle size distribution in the case of complete decomposition into monomers. (A) The limiting case of energy equipartition $〈 E_{k} 〉 = const$ for all species. The solid black line corresponds to the system with complete kinetic coefficients [9], and the solid green line corresponds to the simplified coefficients [12] with the same degree of homogeneity, $2 μ = 1 / 6$ . The dashed line has slope $n_{k} \propto k^{- 19 / 12}$ , predicted by the theory. (B) The limiting case of equal velocity dispersion for all species, $〈 v_{k}^{2} 〉 = const$ . The solid black line refers to the system with the complete kinetic coefficients [11], and the solid green line corresponds to the simplified coefficients [12] with the same degree of homogeneity, $2 μ = 2 / 3$ . The dashed line shows the theoretical dependence, $n_{k} \propto k^{- 11 / 6}$ . In all cases the power-law distribution turns into an exponential decay for large sizes.

Kernels [9] and [11] with homogeneity degree

ϰ = 1 / 6

and

ϰ = 2 / 3

correspond to two limiting cases of the size dependence of the average kinetic energy

〈 E_{k} 〉 = (1 / 2) m_{k} 〈 v_{k}^{2} 〉 \sim k^{β}

. Namely,

β = 0

corresponds to

ϰ = 1 / 6

and

β = 1

to

ϰ = 2 / 3

. Physically, we expect that β is constrained within the interval

0 \leq β \leq 1

. Indeed, negative β would imply vanishing velocity dispersion for very large particles, which is possible only for the unrealistic condition of the collision-free motion. The condition

β > 1

is unrealistic as well. We conclude that β must be limited within the interval

[0,1]

, and therefore

μ = ϰ / 2

varies in the interval

1 / 12 \leq μ \leq 1 / 3

.

Power-Law Collisional Decomposition and Erosion.

Next we consider more realistic models, where collisions with some size distribution of debris and erosion take place.

Power-law decomposition.

Experiments (30, 31) show that the number of debris particles of size k produced in the fragmentation of a particle of size i scales as

x_{k} (i) \sim B (i) / k^{α}

. If the distribution

x_{k} (i)

of the debris size is steep enough, the emerging steady-state particle distribution should be close to that for complete fragmentation into monomers. A scaling analysis, outlined below, confirms this expectation, provided that

α > 2

; moreover, in this case

B (i) = i

(SI Text).

Substituting the debris size distribution

x_{k} (i) \sim i / k^{α}

into the basic kinetic equations [2], we note that the equation for the monomer production rate coincides with Eq. 4, up to a factor in the coefficients

A_{i j}

. At the same time, the general equations [2] for

n_{k}

have the same terms as Eq. 3 for complete decomposition, but with two extra terms—the fourth and fifth terms in Eq. 2. These terms describe an additional gain of particles of size k due to decomposition of larger aggregates. Assuming that the steady-state distribution has the same form as for monomer decomposition,

n_{k} \sim k^{- γ} e^{- ν k}

, one can estimate (up to a factor) these extra terms for the homogeneous kinetic coefficients,

A_{i j} = λ C_{i j} \sim {(i j)}^{μ}

(Eq. 12). One gets

\sum_{i = 1}^{k} \sum_{j = k + 1}^{\infty} A_{i j} n_{i} n_{j} B (j) k^{- α} \sim k^{- α}

[16]

\sum_{i, j \geq k} A_{i j} n_{i} n_{j} [B (i) + B (j)] k^{- α} \sim k^{μ - γ + 1 - α} .

[17]

Here we also require that

ν k ≪ 1

, which is the region where the size distribution exhibits a power-law behavior. The above terms are compared with the other three terms in Eq. 2 or Eq. 3, which are the same for monomer and power-law decomposition:

\sum_{i + j = k} C_{i j} n_{i} n_{j} \sim \sum_{i \geq 1} C_{i k} n_{i} n_{k} \sim \sum_{i \geq 1} A_{i k} n_{i} n_{k} \sim k^{μ - γ} .

[18]

If the additional terms [16] and [17] were negligible, compared with the terms [18] that arise for both models, the emergent steady-state size distributions would be the same. For

k ≫ 1

, one can neglect [16] and [17] compared with [18] if

α > γ - μ

and

α > 1

. Taking into account that the equations for the monomers for the two models coincide when

α > 2

, we arrive at the following criterion for universality of the steady-state distribution:

α > max {γ - μ, 2}

. In the case of complete decomposition into monomers we have

γ = μ + 3 / 2

([15]). Hence the above criterion becomes

α > 2

. In other words, if

α > 2

, the model of complete decomposition into monomers yields the same steady-state size distribution as the model with any power-law distribution of debris.

Collisions with erosion.

In collisions with erosion only a small fraction of a particle mass is chipped off (31, 37, 38). Here we consider a simplified model of such collisions: It takes place when the relative kinetic energy exceeds the threshold energy

E_{eros}

, which is smaller than the fragmentation energy

E_{frag}

. Also, we assume that the chipped-off piece always contains l monomers. Following the same steps as before one can derive rate equations that describe both disruptive and erosive collision. For instance, for complete decomposition into monomers the equation for

n_{k}

with

k \geq l + 2

acquires two additional terms,

λ_{e} \sum_{i = 1}^{\infty} C_{i k + l} n_{i} n_{k + l} - λ_{e} \sum_{i = 1}^{\infty} C_{i k} n_{i} n_{k},

with similar additional terms for

l + 2 > k > 1

and for the monomer equation. Here

λ_{e}

gives the ratio of the frequencies of aggregative and erosive collisions, which may be expressed in terms of

E_{eros}

(SI Text). We assume that

λ_{e}

is small and is of the same order of magnitude as λ. We also assume that

λ_{e}

is constant and that

l λ_{e} ≪ 1

. Then we can show that for

k ≫ 1

the size distribution of aggregates

n_{k}

has exactly the same form, Eq. 15, as for the case of purely disruptive collisions (SI Text).

Universality of the steady-state distribution.

The steady-state size distribution of aggregates [15] is generally universal: It is the same for all size distributions of debris, with a strong dominance of small fragments, independently of its functional form. Moreover, it may be shown analytically (SI Text) that the form [15] of the distribution persists when collisions with erosion are involved. We checked this conclusion numerically, solving the kinetic equations [2] with a power-law, exponential size distribution of debris and for collisions with an erosion (Fig. 2). We find that the particle size distribution [15] is indeed universal for steep distributions of debris size. Fig. 2 also confirms the condition of universality of the distribution [15], if

α > 2

for power-law debris size distributions.

Fig. 2.

Particle size distribution in the case of power-law decomposition. The black solid line depicts the particle size distribution for the case of complete decomposition. The other solid lines show the emergent steady-state distributions for the following size distributions of debris: the power-law distribution, $x_{k} (i) \sim k^{- α}$ , with $α = 1$ (gray), $α = 1.5$ (red), $α = 2$ (blue), $α = 3$ (green), exponential distribution $x_{k} (i) \sim exp (- k)$ (yellow), and collisions with an erosion, $λ_{e} = 0.05$ , $l = 3$ (violet). The dashed lines indicate the corresponding fit $n_{k} \sim k^{- γ}$ . Note that for steep size distributions of debris (power law with $α > 2$ , exponential distribution) all slopes coincide with the one for the case of complete decomposition into monomers. All curves correspond to the case of constant kinetic coefficients with $λ = C_{i j} / A_{i j} = 0.01$ .

A steep distribution of debris size, with strong domination of small fragments, appears plausible since the aggregates are relatively loose objects, with a low average coordination number.

Size Distribution of the Ring Particles.

The distribution of the ring particles’ radii,

F (R)

, is constrained by space- and earth-bound observations (16). To extract

F (R)

we use the relation

R^{3} = r_{k}^{3} = k r_{1}^{3}

(for spherical particles) in conjunction with

n_{k} d k = F (R) d R

. We find that

n_{k} \sim k^{- 3 / 2 - μ} \exp (- λ^{2} k)

implies

F (R) \sim R^{- q} e^{- {(R / R_{c})}^{3}}, q = \frac{5}{2} + 3 μ

[19]

R_{c} = \frac{r_{1}}{λ^{2 / 3}} .

[20]

Thus, for

R ≪ R_{c}

the distribution is algebraic with exponent

q = 5 / 2 + 3 μ

, and the crossover to exponential behavior occurs at

R \sim R_{c}

.

We have shown that the exponent μ can vary within the interval

1 / 12 \leq μ \leq 1 / 3

, and hence the exponent q for the size distribution varies in the range

2.75 \leq q \leq 3.5

. This is in excellent agreement with observations, where values for q in the range from 2.70 to 3.11 were reported (15, 17). Fitting the theory to the particle size distribution of Saturn’s A ring inferred from data obtained by the Voyager Radio Science Subsystem (RSS) during a radio occultation of the spacecraft by Saturn’s A ring (15), we find

R_{c} = 5.5

m (Fig. 3). For

r_{1}

in the plausible range from 1 cm to 10 cm (16) we get (Eq. 20) λ on the order of

10^{- 4}

–

10^{- 3}

, which is the ratio of the frequencies of disruptive and coagulating collisions. It is also possible to estimate characteristic energies and the strength of the aggregates. Using the plausible range for random velocity,

v_{0} = 0.01 - 0.1 cm / s

(16), we obtain values that agree with the laboratory measurements (SI Text).

Fig. 3.

Particle size distribution for Saturn’s A ring. The dashed line represents the particle size distribution of Saturn’s A ring inferred from data obtained by the Voyager Radio Science Subsystem (RSS) during a radio occultation of the spacecraft by the rings (15). A fit of the theoretical curve, Eq. 19, is shown as a solid line, giving a cutoff radius of $R_{c} = 5.5 m$ .

SI Text

Collision Integrals

The collision integral for aggregative collisions has the following form:

\begin{array}{l} I_{k}^{agg} ({\vec{v}}_{k}) = \frac{1}{2} \sum_{i + j = k} σ_{i j}^{2} \int d {\vec{v}}_{i} \int d {\vec{v}}_{j} \int d \vec{e} Θ (- {\vec{v}}_{i j} \cdot \vec{e}) | {\vec{v}}_{i j} \cdot \vec{e} | \\ \times f_{i} ({\vec{v}}_{i}) f_{j} ({\vec{v}}_{j}) Θ (E_{agg} - E_{i j}) δ (m_{k} {\vec{v}}_{k} - m_{i} {\vec{v}}_{i} - m_{j} {\vec{v}}_{j}) \\ - \sum_{j} σ_{k j}^{2} \int d {\vec{v}}_{j} \int d \vec{e} Θ (- {\vec{v}}_{k j} \cdot \vec{e}) | {\vec{v}}_{k j} \cdot \vec{e} | \\ \times f_{k} ({\vec{v}}_{k}) f_{j} ({\vec{v}}_{j}) Θ (E_{agg} - E_{k j}) = I_{k}^{agg, 1} - I_{k}^{agg, 2} . \end{array}

[S1]

The first sum on the right-hand side of Eq. S1 refers to collisions where an aggregate of mass k is formed from smaller aggregates of masses i and j, whereas the second sum describes the collisions of k aggregates with all other particles. In the first sum

m_{k} = m_{i} + m_{j}

and

m_{k} {\vec{v}}_{k} = m_{i} {\vec{v}}_{i} + m_{j} {\vec{v}}_{j}

due to mass and momentum conservation. The rest of the notation is standard (e.g., ref. 24):

σ_{i j}^{2} = r_{1} {(i^{1 / 3} + j^{1 / 3})}^{2}

quantifies the collision cross-section and

| {\vec{v}}_{i j} \cdot \vec{e} |

is the length of the collision cylinder, where the unit vector

\vec{e}

specifies the direction of the intercenter vector at the collision instant;

Θ (- {\vec{v}}_{i j} \cdot \vec{e})

selects only approaching particles. The factor

Θ (E_{agg} - E_{i j})

in the integrands guarantees that the relative kinetic energy does not exceed

E_{agg}

to cause the aggregation. The kinetic energy of the relative motion is

E_{i j} = \frac{1}{2} μ_{i j} {\vec{v}}_{i j}^{2}

, with the relative velocity

{\vec{v}}_{i j} = {\vec{v}}_{i} - {\vec{v}}_{j}

and reduced mass

μ_{i j} = m_{i} m_{j} / (m_{i} + m_{j})

. The notations in the second sum on the right-hand side of Eq. S1 have a similar meaning.

For the collisions leading to fragmentation we have

\begin{array}{l} I_{k}^{frag} ({\vec{v}}_{k}) = \frac{1}{2} \sum_{i, j \geq k + 1} σ_{i j}^{2} \int d {\vec{v}}_{j} \int d {\vec{v}}_{i} \int d \vec{e} Θ (- {\vec{v}}_{i j} \cdot \vec{e}) \\ \times | {\vec{v}}_{i j} \cdot \vec{e} | f_{j} ({\vec{v}}_{j}) f_{i} ({\vec{v}}_{i}) Θ (E_{i j}^{n} - E_{frag}) \\ \times (q_{k i} ({\vec{v}}_{k}, {\vec{v}}_{i}, {\vec{v}}_{j}) + q_{k j} ({\vec{v}}_{k}, {\vec{v}}_{i}, {\vec{v}}_{j})) \\ + \sum_{i = 1}^{k} \sum_{j \geq k + 1} σ_{i j}^{2} \int d {\vec{v}}_{j} \int d {\vec{v}}_{i} \int d \vec{e} Θ (- {\vec{v}}_{i j} \cdot \vec{e}) \\ \times | {\vec{v}}_{i j} \cdot \vec{e} | f_{j} ({\vec{v}}_{j}) f_{i} ({\vec{v}}_{i}) Θ (E_{i j}^{n} - E_{frag}) q_{k j} ({\vec{v}}_{k}, {\vec{v}}_{i}, {\vec{v}}_{j}) \\ - \sum_{i} (1 - δ_{k, 1}) σ_{k i}^{2} \int d {\vec{v}}_{i} \int d \vec{e} Θ (- {\vec{v}}_{k i} \cdot \vec{e}) | {\vec{v}}_{k i} \cdot \vec{e} | \\ \times f_{k} ({\vec{v}}_{k}) f_{i} ({\vec{v}}_{i}) Θ (E_{k i}^{n} - E_{frag}) = I_{k}^{frag, 1} + I_{k}^{frag, 2} - I_{k}^{frag, 3}, \end{array}

[S2]

where we define the kinetic energy of the relative normal motion,

E_{i j}^{n} = (1 / 2) μ_{i j} {({\vec{v}}_{i j} \cdot \vec{e})}^{2}

, with

({\vec{v}}_{i j} \cdot \vec{e})

being the normal relative velocity. In contrast to the case of aggregation where both normal and tangential components must be small, so that the total energy of the relative motion

E_{i j}

matters, for fragmentation only the relative normal motion is important: Only normal motion causes a compression and the subsequent breakage of particles’ material. Hence the kinetic energy of the relative normal motion,

E_{i j}^{n}

must exceed some threshold. The first sum in Eq. S2 describes the collision of particles

i > k

and

j > k

with the relative kinetic energy of the normal motion above the fragmentation threshold

E_{frag}

; both particles give rise to fragments of size k. Further,

q_{k i} ({\vec{v}}_{k}, {\vec{v}}_{i}, {\vec{v}}_{j})

indicates the number of debris of mass

m_{k} = m_{1} k

with the velocity

{\vec{v}}_{k}

, when a particle of mass

m_{i} = m_{1} i

disintegrates in a collision with a particle of mass

m_{j} = m_{1} j

, provided that the precollision velocities are

{\vec{v}}_{i}

and

{\vec{v}}_{j}

. Obviously,

q_{k i} = 0

if

k \geq i

. The function

q_{k i} ({\vec{v}}_{k}, {\vec{v}}_{i}, {\vec{v}}_{j})

depends on a particular collision model. The second sum describes the process, when only one particle with

j > k

(but not with

i < k

) gives rise to debris of size k. Finally, the third term accounts for the breakage of particles of size

k > 1

in collisions with all other particles.

In the present study we focus on the evolution of particle concentrations

n_{k}

. Therefore, (i) the particular forms of two other collision integrals—for bouncing collisions,

I_{k}^{b}

, and for the one describing viscous heating,

I_{k}^{heat}

—are not important, because these terms do not change densities of the species and (ii) it is sufficient to use a more simple distribution,

x_{k} (i)

, defined as

\int q_{k i} ({\vec{v}}_{k}, {\vec{v}}_{i}, {\vec{v}}_{j}) d {\vec{v}}_{k} = x_{k} (i),

[S3]

where

x_{k} (i)

gives the total number of fragments of size k in all possible disruptive collisions of a particles of size

i > k

. In Eq. S3 we exploit the “mean-field” approximation; that is, we assume that the averaged distribution

x_{k} (i)

depends neither on the size of the colliding partner j nor on the velocities of

{\vec{v}}_{i}

and

{\vec{v}}_{j}

, provided a fragmentation occurs.

Note, that whereas

I_{k}^{b}

describes the loss of energy in dissipative collisions,

I_{k}^{heat}

characterizes the energy input, so that the average kinetic energy of all species

〈 E_{k} 〉

,

k = 1,2, \dots

, is kept in a steady state.

Maxwell Approximation for the Velocity Distribution Functions

The interplay between aggregation and fragmentation results in a dynamically sustained mixture of particles of different mass. Mixtures of dissipative particles, generally, have different velocity dispersion or mean kinetic energy (“granular temperature”) of each species. The partial number density (concentration)

n_{i}

and the respective mean kinetic energy

〈 E_{i} 〉

of the species read (e.g., ref. 54)

n_{i} = \int d {\vec{v}}_{i} f_{i} ({\vec{v}}_{i}), n_{i} 〈 E_{i} 〉 = \int \frac{m_{i} {\vec{v}}_{i}^{2}}{2} f_{i} ({\vec{v}}_{i}) d {\vec{v}}_{i} .

[S4]

We assume that the distribution function

f_{i} ({\vec{v}}_{i}, t)

may be written as (11, 54, 55)

f_{i} ({\vec{v}}_{i}, t) = \frac{n_{i} (t)}{v_{0, i}^{3} (t)} ϕ_{i} ({\vec{c}}_{i}), {\vec{c}}_{i} \equiv \frac{{\vec{v}}_{i}}{v_{0, i}},

[S5]

where

v_{0, i}^{2} (t) = 2 〈 E_{i} 〉 (t) / m_{i}

is the thermal velocity and

ϕ (c_{i})

is the reduced distribution function. For force-free granular mixtures (56) and interacting particles (which suffer ballistic annihilation) (57) the reduced distribution function is represented in the form of the Sonine polynomial expansion,

ϕ_{i} (\vec{c}) = ϕ_{M} (\vec{c}) [1 + \sum_{k = 1}^{\infty} a_{k}^{(i)} S_{k} (c^{2})] .

Here

ϕ_{M} (\vec{c})

is the Maxwell distribution function,

ϕ_{M} (\vec{c}) = π^{- 3 / 2} \exp (- c^{2}),

[S6]

and

S_{k} (c^{2})

are the Sonine polynomials. The coefficients

a_{k}^{(i)}

have been computed in a few systems. In all examples they were rather small, e.g., in the case of dissipative collisions (56) and in the case of reacting particles (57). Therefore, we use the Maxwell distribution function [S6] in all further calculations. Integration of Eq. 21 from the main text over

{\vec{v}}_{k}

with the use of the Maxwell distribution function [S6] is rather straightforward, because all arising integrals are Gaussian. This integration, discussed in detail in the next section, yields Eqs. 1 and 2 of the main text.

Derivation of the Rate Coefficients

To derive the rate equations [2] of the main text, we integrate the Boltzmann Eq. 21 of the main text over

{\vec{v}}_{k}

. Because

n_{k} = \int d {\vec{v}}_{k} f_{k} ({\vec{v}}_{k}, t)

, the left-hand side of the Boltzmann equation turns then into

d n_{k} / d t

and gives the rate of change of the concentrations

n_{k}

. The right-hand side gives the contributions to

d n_{k} / d t

from different parts of the collision integral. Because bouncing collisions and the heating term do not change the number of particles, we easily obtain (24)

\int d {\vec{v}}_{k} I_{k}^{b} = \int d {\vec{v}}_{k} I_{k}^{heat} = 0.

We use the Maxwellian distribution for the distribution function

f_{k}

,

f_{k} ({\vec{v}}_{k}, t) = \frac{n_{k}}{π^{3 / 2} v_{0, k}^{3}} e^{- v^{2} / v_{0, k}^{2}},

where

v_{0, k}^{2}

is the thermal velocity of aggregates composed of k monomers. Then the integral over

{\vec{v}}_{k}

of the second part of the aggregation integral,

I_{k}^{agg, 2}

(Eq. S1), may be written as

\begin{matrix} \int d {\vec{v}}_{k} I_{k}^{agg, 2} = \sum_{j} σ_{k j}^{2} \int d {\vec{v}}_{k} \int d {\vec{v}}_{j} \int d \vec{e} Θ (- {\vec{v}}_{k j} \cdot \vec{e}) \\ \times | {\vec{v}}_{k j} \cdot \vec{e} | f_{k} ({\vec{v}}_{k}) f_{j} ({\vec{v}}_{j}) Θ (E_{agg} - E_{k j}) \\ = \sum_{j} \frac{σ_{k j}^{2} n_{k} n_{j}}{π^{3} v_{0, k}^{3} v_{0, j}^{3}} \int d {\vec{v}}_{k} d {\vec{v}}_{j} d \vec{e} Θ (- {\vec{v}}_{k j} \cdot \vec{e}) | {\vec{v}}_{k j} \cdot \vec{e} | \\ \times e^{- v_{k}^{2} / v_{0, k}^{2} - v_{j}^{2} / v_{0, j}^{2}} Θ (E_{agg} - \frac{1}{2} μ_{k j} v_{k j}^{2}), \end{matrix}

[S7]

where

μ_{k j} = m_{k} m_{j} / (m_{k} + m_{j})

,

{\vec{v}}_{k j} = {\vec{v}}_{k} - {\vec{v}}_{j}

, and

E_{k j} = μ_{k j} v_{k j}^{2} / 2

. The integrals in the above equation are Gaussian and hence may be straightforwardly calculated. We perform this calculation for a particular pair k and j. With the substitute

{\vec{v}}_{k} = \vec{u} + \vec{w} (\frac{μ_{k j}}{m_{k}} - p_{k j})

{\vec{v}}_{j} = \vec{u} - \vec{w} (\frac{μ_{k j}}{m_{j}} + p_{k j}),

where

p_{k j} = μ_{k j} [{(m_{k} v_{0, k}^{2})}^{- 1} - {(m_{j} v_{0, j}^{2})}^{- 1}] / [v_{0, k}^{- 2} + v_{0, j}^{- 2}]

, the above integral with respect to

{\vec{v}}_{k}

,

{\vec{v}}_{j}

, and

\vec{e}

may be written as

\int d \vec{u} d \vec{w} d \vec{e} Θ (- \vec{w} \cdot \vec{e}) | \vec{w} \cdot \vec{e} | e^{- a u^{2} - b w^{2}} Θ (E_{agg} - \frac{1}{2} μ_{k j} w^{2}),

where

a \equiv a_{k j} = v_{0, k}^{- 2} + v_{0, j}^{- 2}

and

b \equiv b_{k j} = 1 / (v_{0, k}^{2} + v_{0, j}^{2})

and we take into account that the Jacobian of transformation from

({\vec{v}}_{k}, {\vec{v}}_{j})

to

({\vec{u}}_{k}, \vec{w})

is equal to unity. Integration over

\vec{u}

gives

{(π / a)}^{3 / 2}

. Integration over the unit vector

\vec{e}

gives

4 π

and we are left with the integral over

\vec{w}

. Integration over directions of the vector

\vec{w}

gives π, so finally we need to calculate the remaining integral:

h_{k j} = \int_{0}^{\sqrt{2 E_{agg} / μ_{k j}}} d w w^{3} e^{- b w^{2}} = \frac{1}{2 b^{2}} [1 - e^{2 b E_{agg} / μ_{k j}} (1 + \frac{2 b E_{agg}}{μ_{k j}})] .

As the result we obtain

\int d {\vec{v}}_{k} I_{k}^{agg, 2} = \sum_{j} \frac{σ_{k j}^{2} n_{k} n_{j}}{π^{3} v_{0, k}^{3} v_{0, j}^{3}} 4 π^{2} {(\frac{π}{a_{k j}})}^{3 / 2} h_{k j} .

[S8]

When we integrate the first part of

I_{k}^{agg, 1}

in Eq. S1 over

{\vec{v}}_{k}

, we observe that

\int d {\vec{v}}_{k} δ (m_{k} {\vec{v}}_{k} - m_{i} {\vec{v}}_{i} - m_{j} {\vec{v}}_{j}) = 1

, because the other part of the integrand does not depend on

{\vec{v}}_{k}

. Then the remaining integration is exactly the same as for

I_{k}^{agg, 2}

, which has been already performed; therefore we find

\int d {\vec{v}}_{k} I_{k}^{agg, 1} = \frac{1}{2} \sum_{i + j = k} \frac{σ_{i j}^{2} n_{i} n_{j}}{π^{3} v_{0, i}^{3} v_{0, j}^{3}} 4 π^{2} {(\frac{π}{a_{i j}})}^{3 / 2} h_{i j} .

[S9]

Turn now to the calculation of the integral over

{\vec{v}}_{k}

of the fragmentation integral

I_{k}^{frag}

. We note that according to Eq. S3, the integration over

{\vec{v}}_{k}

of

q_{k i} ({\vec{v}}_{k}, {\vec{v}}_{i}, {\vec{v}}_{j})

yields

x_{k} (i)

. Therefore, for all three parts of

\int d {\vec{v}}_{k} I_{k}^{frag}

we need to compute the integrals

\begin{array}{l} \int d {\vec{v}}_{j} \int d {\vec{v}}_{i} \int d \vec{e} Θ (- {\vec{v}}_{i j} \cdot \vec{e}) \\ \times | {\vec{v}}_{i j} \cdot \vec{e} | f_{i} ({\vec{v}}_{i}) f_{j} ({\vec{v}}_{j}) Θ (E_{i j}^{n} - E_{frag}), \end{array}

[S10]

which have the same structure as the integrals in [S7]. The only difference is that instead of the factor

Θ (E_{agg} - E_{i j})

in [S7] we have now

Θ (E_{i j}^{n} - E_{frag})

. Therefore, calculations of the integrals [S10] may be performed as done previously. We apply the same transformation of variables as before and arrive at the integral

\int d \vec{u} d \vec{w} d \vec{e} Θ (- \vec{w} \cdot \vec{e}) | \vec{w} \cdot \vec{e} | e^{- a u^{2} - b w^{2}} Θ (E_{frag} - \frac{μ_{i j} w_{n}^{2}}{2}),

where

w_{n} = (\vec{w} \cdot \vec{e})

, and the same notations as previously are used. Integration over

\vec{u}

and

\vec{e}

gives again

4 π {(π / a)}^{3 / 2}

. If we choose the direction of the vector

\vec{e}

along the z axis (the direction of the z axis is arbitrary), integration over

w_{x}

and

w_{y}

yields

{(π / b)}^{1 / 2} \cdot {(π / b)}^{1 / 2}

and integration over

w_{z}

leads to the integral

g_{i j} = \int_{- \infty}^{- \sqrt{2 E_{frag} / μ_{i j}}} d w_{z} | w_{z} | e^{- b w_{z}^{2}} = \frac{1}{2 b} e^{- (2 b E_{frag} / μ_{i j})} .

Hence we obtain

\begin{matrix} \int d {\vec{v}}_{k} I_{k}^{frag} = \frac{1}{2} \sum_{i, j \geq k + 1} (x_{k i} + x_{k j}) \frac{σ_{i j}^{2} n_{i} n_{j}}{π^{3} v_{0, i}^{3} v_{0, j}^{3}} 4 π {(\frac{π}{a_{i j}})}^{3 / 2} \frac{π}{b_{i j}} g_{i j} \\ + \sum_{i = 1}^{k} \sum_{j \geq k + 1} x_{k j} \frac{σ_{i j}^{2} n_{i} n_{j}}{π^{3} v_{0, i}^{3} v_{0, j}^{3}} 4 π {(\frac{π}{a_{i j}})}^{3 / 2} \frac{π}{b_{i j}} g_{i j} \\ - \sum_{i} (1 - δ_{k, 1}) \frac{σ_{k i}^{2} n_{k} n_{i}}{π^{3} v_{0, k}^{3} v_{0, i}^{3}} 4 π {(\frac{π}{a_{k i}})}^{3 / 2} \frac{π}{b_{k i}} g_{k i} . \end{matrix}

[S11]

Combining Eqs. S8, S9, and S11 and using the above definitions of the quantities

a_{i j}

,

b_{i j}

,

μ_{i j}

,

h_{i j}

, and

g_{i j}

, we arrive at the rate equations [2] with the according rate coefficients [1] of the main text. To write the rate coefficients

C_{i j}

and

A_{i j}

in the form of Eq. 2 of the main text we take into account that the thermal velocity

v_{0 i}

is related to the mean square velocity

〈 v_{0}^{2} 〉

, termed here as the velocity dispersion, as

〈 v_{0}^{2} 〉 = \frac{3}{2} v_{0, i}^{2},

which is a direct consequence of the Maxwellian distribution.

Self-Gravity Wakes and Averaged Kinetic Equations

Saturn’s rings are not uniform but exhibit a large variety of structures (58). One example is the self-gravity wakes (59–62), arising from self-gravitational instability, forming a transient and fluctuating pattern in the surface mass density of the rings. These are canted relative to the azimuthal direction, with a typical length scale of about

L_{w} \sim 10^{2}

m, one Toomre critical wavelength (63). Thus, to describe adequately particle kinetics one needs, in principle, to take into account effects of dense packing and nonhomogeneity.

Two important comments are needed in this respect. First, due to the low-velocity dispersion of particles in the dense parts of the wakes, the collision duration is still significantly smaller than the time between particle collisions. This implies the validity of the assumption of binary collisions, as the dominant mechanism of particles’ kinetics. Therefore, a kinetic description in terms of the Enskog–Boltzmann equation is possible (25). Although this Markovian equation ignores memory effects in particle kinetics, it may be still applicable, when the mean free path is comparable to or even smaller than the particle size (64). Second, the characteristic length scale of the density wakes,

L_{w}

, and the upper cutoff radius (less then

10 m

) are well separated. This allows one to neglect variations of density and distribution functions on the latter length scale and use a local approximation for the distribution function of two particles at contact:

f_{2} ({\vec{v}}_{k}, \vec{r} - \vec{e} r_{k}, {\vec{v}}_{l}, \vec{r} + \vec{e} r_{l}, t) ≃ g_{2} (σ_{l k}) f_{k} ({\vec{v}}_{k}, \vec{r}, t) f_{l} ({\vec{v}}_{l}, \vec{r}, t) .

[S12]

Here

f_{2}

is the two-particle distribution function corresponding to particles of radii

r_{k}

and

r_{l}

, which have a contact at point

\vec{r}

. The unit vector

\vec{e}

joins the centers of particles and

g_{2} (σ_{l k})

, with

σ_{l k} = r_{l} + r_{k}

, is the contact value of the pair distribution function. It may be well approximated by the corresponding equilibrium value for the hard sphere fluid; explicit expressions for

g_{2} (σ_{l k})

can be found in ref. 24. In the local approximation, Eq. S12, the collision integrals depend only on local values at a particular space point

\vec{r}

. This significantly simplifies the kinetic description of a high-density gas, because the density effects are taken into account in this approach by the multiplicative Enskog factor,

g_{2} (σ_{l k})

only, leaving the structure of the collision integrals unchanged. Therefore, the Enskog–Boltzmann equation valid for the case of wakes reads

\begin{array}{l} \frac{\partial}{\partial t} f_{k} ({\vec{v}}_{k}, \vec{r}, t) + {\vec{v}}_{k} \cdot \vec{\nabla} f_{k} ({\vec{v}}_{k}, \vec{r}, t) + {\vec{F}}_{k} (\vec{r}) \cdot \frac{\partial}{\partial {\vec{v}}_{k}} f_{k} ({\vec{v}}_{k}, \vec{r}, t) \\ = I_{k}^{agg} (\vec{r}) + I_{k}^{b} (\vec{r}) + I_{k}^{frag} (\vec{r}) . \end{array}

[S13]

Here

f_{k} ({\vec{v}}_{k}, \vec{r}, t)

is the velocity distribution function of particles of size k, which depends on the space coordinate

\vec{r}

. Further,

{\vec{F}}_{k} (\vec{r})

is the total gravitational force, acting on the particle of size k, which includes the gravitational force from the central planet as well as self-gravitation of the ring particles. In what follows we do not need an explicit expression for this term. The collision integrals on the right-hand side of Eq. S13 have the same form as in the previous case of a uniform system, with the only difference being that they depend on local parameters taken at a point

\vec{r}

and that the collision cross-sections are renormalized according to the rule

σ_{i j}^{2} \to σ_{i j}^{2} g_{2} (σ_{i j}),

[S14]

which accounts for the high-density effects in the local approximation (24). We do not need here the term

I_{k}^{heat} (\vec{r})

that mimics the heating for the model of a uniform gas, because Eq. S13 implicitly contains the spacial gradients and fluxes responsible for the heating.

Integrating the kinetic Eq. S13 over

{\vec{v}}_{k}

, we find rate equations for the concentrations

n_{k} (\vec{r})

,

\begin{array}{l} \frac{\partial}{\partial t} n_{k} (\vec{r}) + \vec{\nabla} \cdot {\vec{j}}_{k} (\vec{r}) \\ = \frac{1}{2} \sum_{i + j = k} C_{i j} (\vec{r}) n_{i} (\vec{r}) n_{j} (\vec{r}) \\ - n_{k} (\vec{r}) \sum_{i \geq 1} C_{k i} (\vec{r}) n_{i} (\vec{r}) - \sum_{i \geq 1} A_{k i} (\vec{r}) n_{k} (\vec{r}) n_{i} (\vec{r}) (1 - δ_{k 1}) \\ + \sum_{i = 1}^{k} n_{i} (\vec{r}) \sum_{j \geq k + 1} A_{i j} (\vec{r}) n_{j} (\vec{r}) x_{k} (j) \\ + \frac{1}{2} \sum_{i, j \geq k + 1} A_{i j} (\vec{r}) n_{i} (\vec{r}) n_{j} (\vec{r}) (x_{k} (i) + x_{k} (j)), \end{array}

[S15]

where

{\vec{j}}_{k} (\vec{r}) = \int {\vec{v}}_{k} f_{k} ({\vec{v}}_{k}, \vec{r}, t) d {\vec{v}}_{k}

is the macroscopic (hydrodynamic) flux associated with particles of size k. The important feature of Eq. S15 is the spatial dependence of the kinetic coefficients

A_{i j}

and

C_{i j}

. Although the structure of these coefficients coincides with that of kinetic coefficients in the uniform system (apart from the trivial renormalization, Eq. S14), all quantities here are local. Naturally, the local velocity dispersion

〈 v_{i}^{2} 〉 (\vec{r})

in the dense parts of the wakes significantly differs from that of the dilute regions in between.

Now we average Eq. S15 over a suitable control volume V, which contains a large number of wakes. Applying Green’s theorem,

\frac{1}{V} \int_{V} \vec{\nabla} \cdot {\vec{j}}_{k} d \vec{r} = \frac{1}{V} \int_{S} {\vec{j}}_{k} \cdot d \vec{s} \sim \frac{S}{V},

we note that the contribution of the term containing the flux

{\vec{j}}_{k}

vanishes as

S / V \to 0

for large enough volume. As a result we arrive at the set of equations with the space-averaged quantities

\begin{array}{l} \frac{d}{d t} {\bar{n}}_{k} = \frac{1}{2} \sum_{i + j = k} {\bar{C}}_{i j} {\bar{n}}_{i} {\bar{n}}_{j} - {\bar{n}}_{k} \sum_{i} {\bar{C}}_{k i} {\bar{n}}_{i} \\ - \sum_{i} {\bar{A}}_{k i} {\bar{n}}_{k} {\bar{n}}_{i} (1 - δ_{k 1}) + \sum_{i = 1}^{k} {\bar{n}}_{i} \sum_{j = k + 1}^{\infty} {\bar{A}}_{i j} {\bar{n}}_{j} x_{k} (j) \\ + \frac{1}{2} \sum_{i, j \geq k + 1} {\bar{A}}_{i j} {\bar{n}}_{i} {\bar{n}}_{j} (x_{k} (i) + x_{k} (j)) . \end{array}

[S16]

Here, by definition,

{\bar{n}}_{k} = \frac{1}{V} \int n_{k} (\vec{r}) d \vec{r}

and

{\bar{C}}_{i j} = \frac{1}{{\bar{n}}_{i} {\bar{n}}_{j}} \frac{1}{V} \int n_{i} (\vec{r}) n_{j} (\vec{r}) C_{i j} (\vec{r}) d \vec{r},

where

C_{i j} (\vec{r})

are defined by Eq. 1 of the main text, with the local velocity dispersions

〈 v_{i}^{2} 〉 (\vec{r})

. A similar expression holds true for the coefficients

{\bar{A}}_{i j}

.

It is important to note that the coefficients

{\bar{C}}_{i j}

and

{\bar{A}}_{i j}

are density-weighted quantities. Therefore, the contribution to the average value is proportional to the local density. This in turn implies that the values of these coefficients practically coincide with these for the dense part of the wakes,

{\bar{C}}_{i j} = C_{i j}^{(densepart)}, {\bar{A}}_{i j} = A_{i j}^{(densepart)} .

Hence we conclude that the kinetic equations for the average concentrations of particles

{\vec{n}}_{k}

coincide with the previously derived equations for

n_{k}

for the case of a uniform system. In what follows we use

n_{k}

,

C_{i j}

, and

A_{i j}

for notation brevity, keeping in mind, however, that they correspond to the average values that are almost equal to these values in the dense part of the wakes.

Estimates of the Characteristic Energies and the Aggregates Strength

Using the data for λ reported in the main text, we perform here some estimates.

Estimate of the Fragmentation Energy and of the Aggregates Strength.

First we estimate the effective value of

E_{frag}

, assuming that

〈 B_{i j} E_{agg} (R) 〉 ≫ 1

. According to Eq. 1 of the main text

C_{i j} ≃ ν_{i j}

and

A_{i j} ν_{i j} \exp (- B_{i j} E_{frag})

, which yields

λ = \frac{A_{i j}}{C_{i j}} ≃ \exp (- B_{i j} E_{frag}) .

[S17]

We estimate the average value of

B_{i j}

. For two particles of equal mass

m_{i} = m_{j}

with characteristic square velocity

v_{0}^{2}

this quantity reads, according to Eq. 1 of the main text,

B_{i j} = \frac{3}{2} \frac{(m_{i}^{- 1} + m_{j}^{- 1})}{v_{0}^{2}} = \frac{9}{4 π ρ ϕ v_{0}^{2}} R^{- 3},

[S18]

where

m_{i} = (4 π / 3) ρ ϕ R^{3}

,

ρ = 900 kg / m^{3}

is the material density of ice, and

ϕ = 0.3

is the approximate packing fraction of aggregates. Let us estimate the average value of

B_{i j}

, which we define as

〈 B_{i j} (R) 〉 = \frac{\int_{r_{1}}^{R_{c}} B_{i j} (R) F (R) d R}{\int_{r_{1}}^{R_{c}} F (R) d R} .

Here

F (R) ≃ const . R^{- q}

is the radii distribution function. It behaves as a power law with

q ≃ 3

for

R < R_{c}

, with

R_{c} ≫ r_{1}

being the cutoff radius for the distribution. The averaging for

q = 3

yields

〈 B_{i j} 〉 = \frac{9}{10 π ρ ϕ v_{0}^{2}} r_{1}^{- 3}

and respectively the average fragmentation energy,

〈 E_{frag} 〉 = - \log λ {〈 B_{i j} 〉}^{- 1} = \frac{5}{3} π ρ ϕ v_{0}^{2} r_{1}^{3} log (\frac{R_{c}}{r_{1}}) .

[S19]

Here we use Eq. 20 of the main text,

R_{c} = r_{1} / λ^{2 / 3}

. As may be seen from the above equation, the estimate of the effective fragmentation energy sensitively depends on the monomer size

r_{1}

and the characteristic square velocity

v_{0}^{2}

. The plausible range for these values is

1 \leq r_{1} \leq 10 cm

and

0.01 \leq v_{0} \leq 0.1 cm / s

(16).

To be consistent with the laboratory measurements we choose the particular values

r_{1} = 7 cm

and

v_{0} = 0.07 cm / s

from the above intervals for

r_{1}

and

v_{0}

(other combinations of these parameters are also consistent with the laboratory data) to obtain

〈 E_{frag} 〉 = 1.04 \cdot 10^{- 6} J .

This fragmentation energy is equal to the product of the fragmentation energy of a single contact between aggregates

E_{b}

times the average number of contacts between monomers

〈 N_{c} 〉

in the aggregates. It may be estimated as follows. If the radius of an aggregate composed of monomers of radius

r_{1}

is R and the packing fraction is ϕ, the number of contacts reads

N_{c} (R) = {(R / r_{1})}^{3} ϕ z_{c}

, where

z_{c}

is the average number of contacts with neighbors. For a random packing of spheres

z_{c} = 4.7

(65) and the averaged contact number is

〈 N_{c} 〉 = \frac{\int_{r_{1}}^{R_{c}} N_{c} (R) F (R) d R}{\int_{r_{1}}^{R_{c}} F (R) d R} = 2 ϕ z_{c} (\frac{R_{c}}{r_{1}}) .

Here we again use

F (R) ≃ const . R^{- q}

with

q = 3

. Hence we obtain the average contact energy for a monomer–monomer bond:

E_{b} = \frac{〈 E_{frag} 〉}{〈 N_{c} 〉} = 4.68 \cdot 10^{- 9} J .

Now we assume that the adhesive contacts between the monomers occur in accordance with the overlapping frost layer model, as follows from the laboratory measurements of refs. 66 and 67. This model of cohesion has been also used in the numerical simulation of the Saturn rings, where aggregation and fragmentation processes have been taken into account (51, 52). The typical thickness of the frost layer is about

d = 20 μ

, which yields the estimate of the adhesion force

f_{b}

,

f_{b} = \frac{E_{b}}{d} = 23.4 \cdot 10^{- 5} N = 23.4 dyne,

in good agreement with the laboratory data (51, 52), where the force of the order

30 - 50 dyne

has been reported.

The contact area

S_{b}

between the monomers, composing an aggregate, is equal to (51)

S_{b} = π r_{1}^{2} β (1 - \frac{β}{4}),

where the parameter

β = d / r_{1}

has been introduced in ref. 51 and characterizes the ratio of the frost layer thickness and the particle radius. From the laboratory experiments it follows that

β = 10^{- 3}

(51). This gives the estimate of the strength of icy aggregates:

P_{b} = \frac{f_{b}}{S_{b}} \sim 1.5 \cdot 10^{6} Pa .

Indeed, this stress exists in the aggregates and keeps the constituents together. Therefore, it is reasonable to assume that any external stress smaller than

P_{b}

would not destroy it. However, the applied external stress exceeding

P_{b}

would most probably cause fragmentation.

Estimate of the Aggregation Energy.

Here we estimate the value of

E_{agg}

for the overlapping frost layer model for particles contact. First we estimate the apparent adhesive coefficient γ, which is equal to twice the surface free energy per unit area of a solid in vacuum. We use the laboratory data of refs. 66 and 67 performed on the ice particles of radius

R_{0} = 2.5 cm

. The characteristic force due to the frost layer of the thickness

d = 10 - 30 μ

was

f_{0} = 30 - 50 dynes

; therefore the according energy,

f_{0} \cdot d

, is about

40 \cdot 10^{- 5} N \times 20 \cdot 10^{- 6} m = 8 \cdot 10^{- 9} J

. Therefore, the apparent coefficient of adhesion reads

γ_{eff} = \frac{f_{0} d}{π R_{0}^{2} β} = 0.0041 J / m .

Now we wish to estimate the effective Young modulus and Poisson ratio of aggregates composed of random packing of monomers, which have cohesive bonds as follows from the overlapping frost layer model.

Generally it may be shown that the effective Young modulus

Y_{eff}

and Poisson ratio

ν_{eff}

of a random packing of spheres of radius R that interact with a force

f (r)

read

Y_{eff} = {\frac{3 π^{2}}{10} R^{- 1 / 3} ϕ z_{c} \frac{d f (r)}{d r} |}_{r = r_{eq}}, ν_{eff} = \frac{1}{4},

where ϕ is the packing fraction and

r_{eq}

is the equilibrium distance between the spheres’ centers. Note that the effective Poison ratio does not depend on the microscopic detail and is determined by the geometry of random packing only. To derive the above relations, one needs to consider two different types of deformation: uniform compression (without shear) and uniform shear (without change of a volume). Then it is straightforward to compute the change of the system’s elastic energy in terms of the deformations. This is done for the case of randomly packed spheres, with a given interparticle force, and an average over configurations must be performed. The linear coefficients that relate the variation of the elastic energy to the respective deformation (for both types) yield the elastic coefficients, that is, the Young modulus and the Poisson ratio, as given in the above equation.

Using the experimental value of

d f (r) / d r = 5.5 \cdot 10^{4} dyne / cm

from ref. 66 and

R = r_{1} = 7 cm

, we obtain the effective Young modulus,

Y_{eff} = 557 Pa

for the aggregate material, treated as an elastic continuum medium.

With the obtained effective values for the Young modulus, the Poisson ratio, and the adhesive coefficient we can apply the effective Johnson, Kendall, and Roberts model, treating the colliding aggregates as continuum bodies with the effective material parameters (that is, ignoring the discrete structure of the aggregates). Then the threshold energy of aggregation for two particles of radii

R_{i}

and

R_{j}

has the following form (68)

E_{agg} = q_{0} {(π^{5} γ^{5} R_{eff}^{4} D^{2})}^{1 / 3} .

[S20]

Here

q_{0} = 1.457

is a constant; γ is the adhesion coefficient;

D = (3 / 2) (1 - ν^{2}) / Y

, where Y and ν are, respectively, the Young modulus and Poisson ratio; and

R_{eff} = R_{i} R_{j} / (R_{i} + R_{j})

. In the case of interest we need to use the effective values for all of the parameters. To estimate

〈 E_{agg} 〉

we apply the above expression for particles of equal radius to obtain

〈 E_{agg} 〉 = \frac{\int_{r_{1}}^{R_{c}} E_{agg} (R) F (R) d R}{\int_{r_{1}}^{R_{c}} F (R) d R} = 3 q_{0} {(π^{5} γ_{eff}^{5} {(\frac{r_{1}}{2})}^{2} D_{eff}^{2})}^{1 / 3},

where

D_{eff} = (3 / 2) (1 - ν_{eff}^{2}) / Y_{eff}

. This gives

〈 E_{agg} 〉 = 6.21 \cdot 10^{- 7} J,

that is,

〈 E_{frag} 〉 = 1.68 〈 E_{agg} 〉,

which implies that the aggregation and fragmentation energies are rather close.

Finally we compute

〈 B_{i j} E_{agg} (R) 〉

. Because we have already computed

〈 B_{i j} (R) 〉

and

〈 E_{agg} (R) 〉

, we just apply the approximation

〈 B_{i j} E_{agg} (R) 〉 \approx 〈 B_{i j} (R) 〉 〈 E_{agg} (R) 〉,

which yields

exp (- 〈 B_{i j} E_{agg} (R) 〉) \approx 0.02

and justifies the approximation

A_{i j} / C_{i j} = \exp (- B_{i j} E_{frag})

used in the main text.

Universality of Particles Size Distribution for Steep Distribution of Debris

As has been already mentioned in the main text, distribution of debris in a collision obeys in its main part a power law. That is, if an aggregate of size i suffers a disruption in an impact, plenty of fragments of size

k < i

appear. Let

x_{k} (i)

denote the number of fragments of size k; the power-law fragment distribution implies that

x_{k} (i) \sim k^{- α}

in the main part of the distribution. This allows us to quantify the prefactor of the distribution,

x_{k} (i) = B (i) k^{- α}

from the normalization condition, that is, from the condition that the total mass of all debris is equal to the mass of the parent body. Although we have a discrete mass spectrum of debris,

m_{k} = m_{1} k

,

k = 1,2, \dots

, for

i ≫ 1

one can approximate summation by integration to obtain

i ≃ \int_{1}^{i - 1} B k^{- α} k d k = {\begin{array}{l} β \frac{B}{(2 - α)} [{(i - 1)}^{2 - α} - 1] & if α \neq 2 \\ β_{1} B \log (i - 1) & if α = 2, \end{array}

[S21]

where the factors β and

β_{1}

stand for an approximate correction when the summation is approximated by the integration (more details in next section). This yields for

i ≫ 1

B (i) \sim {\begin{array}{l} i^{α - 1} & if α < 2 \\ i {(\log i)}^{- 1} & if α = 2 \\ i & if α > 2 . \end{array}

[S22]

Now we perform analysis of the general system of Eq. 2 in the main text to show that under certain conditions the solution to Eq. 2 (fragmentation with a particular debris-size distribution) coincides with the solution to Eqs. 3 and 4 of the main text (complete fragmentation into monomers). First, we note that if

x_{k} (i) \sim i / k^{α}

, which holds true for

α > 2

, the equations for monomers are identical (up to a factor at the coefficients

A_{i j}

) for both models. Next, we write Eq. 2 as

\frac{d n_{k}}{d t} = K_{1} - K_{2} - K_{3} + K_{4} + K_{5},

[S23]

where

K_{1} = \frac{1}{2} \sum_{i + j = k} C_{i j} n_{i} n_{j}

K_{2} = \sum_{i = 1}^{\infty} C_{i k} n_{i} n_{k}

K_{3} = \sum_{i = 1}^{\infty} A_{k i} n_{i} n_{k} (1 - δ_{k 1})

K_{4} = \sum_{i = 1}^{k} \sum_{j = k + 1}^{\infty} A_{i j} n_{i} n_{j} x_{k} (j)

K_{5} = \frac{1}{2} \sum_{i, j \geq k + 1} A_{i j} n_{i} n_{j} [x_{k} (i) + x_{k} (j)] .

In these notations Eq. 3 of the main text for the case of decomposition into monomers read

\frac{d n_{k}}{d t} = K_{1} - K_{2} - K_{3};

[S24]

that is, the two models differ by the two terms

K_{4}

and

K_{5}

only. Now we estimate the relative magnitude of the terms

(K_{1} - K_{2} - K_{3})

and

K_{4}

and

K_{5}

, using the scaling approach. We assume that under certain conditions, the distribution of aggregate concentrations in a steady state (when

{\dot{n}}_{k} = 0

) has the form

n_{k} \sim k^{- γ} e^{- a k},

the same as for the case of decomposition into monomers. We perform the analysis for the generalized product kernels

C_{i j} = {(i j)}^{μ}

and

A_{i j} = λ C_{i j}

. We are interested in the scaling regime,

k ≫ 1

, and we focus on the power-law domain where

k a ≪ 1

. Additionally, we assume that

1 < γ - μ < 2

; we will check all of the assumptions a posteriori. Approximating again the summation by the integration we obtain for the first term

K_{1}

,

\begin{matrix} K_{1} \sim \sum_{i + j = k} {(i j)}^{μ - γ} e^{- a k} \sim \int_{1}^{k - 1} i^{μ - γ} {(k - i)}^{μ - γ} d i \\ \sim k^{μ - γ} \int_{1}^{k / 2} i^{μ - γ} {(1 - \frac{i}{k})}^{μ - γ} d i \\ \sim k^{μ - γ} {[i^{μ - γ + 1} - b \frac{i^{μ - γ + 2}}{k} \dots]}_{1}^{k / 2} \sim k^{μ - γ}, \end{matrix}

[S25]

where

b = (μ - γ) (μ - γ + 1) {(μ - γ + 2)}^{- 1}

. Here we take into account the symmetry of the integrand around

k / 2

, make an expansion of the factor

{(1 - i / k)}^{μ - γ}

, and keep only the leading term in the obtained series. Now we evaluate the second and the third terms:

K_{2} + K_{3} \sim (1 + λ) k^{μ - γ} e^{- a k} \int_{1}^{\infty} i^{μ - γ} e^{- a i} d i \sim k^{μ - γ} .

[S26]

Similarly, we find for the fourth term,

\begin{matrix} K_{4} \sim λ \sum_{i = 1}^{k} i^{- γ} e^{- a i} \sum_{j = k + 1}^{\infty} {(i j)}^{μ} j^{- γ} e^{- a j} B (j) k^{- α} \\ \sim k^{- α} \int_{1}^{k} i^{μ - γ} e^{- a i} d i \int_{k}^{\infty} j^{μ - γ} e^{- a j} B (j) d j . \end{matrix}

Using

B (j)

from Eq. S22, it is straightforward to show that the forth term scales as

K_{4} \sim {\begin{array}{l} k^{- α} & if α > γ - μ, α \neq 2 \\ k^{- α} {(\log k)}^{- 1} & if α = 2 \\ k^{- α} | \log k a | & if α = γ - μ \\ k^{μ - γ} & if α < γ - μ . \end{array}

[S27]

Analogously, the fifth term is estimated to give

\begin{matrix} K_{5} \sim \sum_{i, j \geq k} {(i j)}^{μ} {(i j)}^{- γ} e^{- a (i + j)} [B (i) + B (j)] k^{- α} \\ \sim k^{- α} \int_{k}^{\infty} d i i^{μ - γ} e^{- a i} \int_{k}^{\infty} j^{μ - γ} e^{- a j} B (j) d j . \end{matrix}

Finally, we get

K_{5} \sim {\begin{array}{l} k^{μ + 1 - α - γ} & if α > γ - μ, α \neq 2 \\ k^{μ + 1 - α - γ} {(\log k)}^{- 1} & if α = 2 \\ k^{μ + 1 - α - γ} | \log k a | & if α = γ - μ \\ k^{1 - 2 (γ - μ)} & if α < γ - μ . \end{array}

[S28]

Comparing Eqs. S26 and S25 with Eqs. S27 and S28, we conclude that the fourth and fifth terms of the basic Eq. 2 of the main text are negligibly small for

k ≫ 1

compared with the first, second, and third terms of this equation, provided

α > (γ - μ)

under the condition

1 < γ - μ < 2

. If we additionally take into account that for

α > 2

the equations for the monomer concentration coincide for the two models, we conclude that if

α > max {γ - μ, 2}

, the steady-state size distribution of aggregates for the case of complete decomposition into monomers and for the power-law decomposition would be the same in the domain

k ≫ 1

and

k a ≪ 1

. Because for the case of monomer decomposition

γ = 3 / 2 + μ

and

a = λ^{2}

, the condition

1 < γ - μ < 2

holds true, and

k a ≪ 1

is fulfilled even for large k if

λ ≪ 1

. Hence it is expected that for a steep size distribution of debris with

α > α_{0} = 2

the steady-state distribution

n_{k} \sim k^{3 / 2 + μ}

[S29]

is universal; that is, it does not depend on the particular value of α.

Moreover, the same conclusion of the universality of the distribution [S29] holds true for any functional form of a steep distribution of debris size. If we write it as

x_{k} (i) = B (i) ϕ (k)

, where the function

ϕ (k)

is steep enough, so that

\int_{1}^{i - 1} ϕ (k) k d k ≃ \int_{1}^{\infty} ϕ (k) k d k = C^{- 1},

the prefactor

B (i)

reads

B (i) = C i

. Then for any function

ϕ (k)

satisfying the condition

ϕ (k) ≪ k^{- 3 / 2}

for

k ≫ 1

the resultant distribution of aggregates will have the universal form [S29]. This has been confirmed numerically for the exponential distribution of debris (Fig. 2 in the main text).

Particle Size Distribution in the Presence of Collisions with Erosion

Here we consider a more general case when in addition to the disruptive collisions that completely destroy aggregates, collisions with an erosion exist. In the erosive collisions a small fraction of the colliding particle mass is chipped off (31, 37). To understand what the impact is of the erosive collisions on the particle size distribution, we consider a simplified model: A disruptive collisions occurs, if the kinetic energy of the relative normal motion

E_{i j}^{n}

exceeds

E_{frag}

, whereas the erosive collision takes place if

E_{i j}^{n}

exceeds smaller energy,

E_{eros}

. That is, the condition for the erosive collision reads

E_{eros} \leq E_{i j}^{n} < E_{frag}

. We assume that in an erosive collision a piece of a fixed size is chipped off from one of the colliding partners. Let the chipped-off piece always contains l monomers. This piece may be further decomposed into smaller fragments. Below we consider two limiting cases, when the chipped-off piece remains intact and when it breaks into l monomers. Performing the same derivation steps as for the fragmentation without erosion, we arrive at the rate equations that may be written for the both cases uniformly. For simplicity we consider here the case of complete decomposition into monomers. Then the equation for monomers reads

\begin{matrix} \frac{d n_{1}}{d t} = - n_{1} \sum_{i \geq 1} C_{1 i} n_{i} + \frac{λ}{2} \sum_{i, j \geq 2} (i + j) C_{i j} n_{i} n_{j} \\ + λ n_{1} \sum_{i \geq 2} C_{1 i} i n_{i} + ϵ λ_{e} l \sum_{i \geq 1} \sum_{j \geq l + 2} C_{i j} n_{i} n_{j}, \end{matrix}

[S30]

where

ϵ = 1

if the chipped-off piece disintegrates into monomers and

ϵ = 0

if it remains intact. For

k \geq l + 2

we obtain

\begin{matrix} \frac{d n_{k}}{d t} = \frac{1}{2} \sum_{i + j = k} C_{i j} n_{i} n_{j} - (1 + λ) \sum_{i \geq 1} C_{i k} n_{i} n_{k} \\ + λ_{e} \sum_{i \geq 1} C_{i k + l} n_{i} n_{k + l} - λ_{e} \sum_{i \geq 1} C_{i k} n_{i} n_{k}, \end{matrix}

[S31]

and, correspondingly for

l + 1 \geq k \geq 2

,

\begin{matrix} \frac{d n_{k}}{d t} = \frac{1}{2} \sum_{i + j = k} C_{i j} n_{i} n_{j} - (1 + λ) \sum_{i \geq 1} C_{i k} n_{i} n_{k} \\ + λ_{e} \sum_{i \geq 1} C_{i k + l} n_{i} n_{k + l} + (1 - ϵ) δ_{k, l} λ_{e} \sum_{i \geq 1} \sum_{j \geq l + 2} C_{i j} n_{i} n_{j} . \end{matrix}

[S32]

The new coefficient

λ_{e}

characterizes the relative frequency of the erosive and aggregative collisions,

λ_{e} = \frac{e^{- B_{i j} E_{eros}} (1 - e^{- B_{i j} (E_{frag} - E_{eros})})}{1 - (1 + B_{i j} E_{agg}) e^{- B_{i j} E_{agg}}},

[S33]

where the coefficients

B_{i j}

are defined in Eq. 1 of the main text. (Note that in the above model of an erosive collision with the chip-off of an intact piece, the monomers do not appear. They are produced in the disruptive collisions.) We illustrate the derivation of size distribution for the case of complete disintegration of the chipped-off piece, which corresponds to

ϵ = 1

in the above equations; the case of

ϵ = 0

is more simple and yields qualitatively the same results. Similar to the analysis of the main text, we assume that the conditions that keep λ and

λ_{e}

constant are fulfilled. Moreover, we also assume that these constants are small and are of the same order of magnitude; that is,

λ_{e} = α λ

, where α is of the order of unity,

α \sim 1

. Because in erosive collisions the value of l is not very large, and

λ_{e}

is small, we further assume that

l λ_{e} ≪ 1

.

For notation simplicity we perform the analysis for the generic case of constant rate coefficients,

C_{i j} = 1

. Similar to that in the main text, the analysis may be undertaken for the more general case of

C_{i j} = {(i j)}^{μ}

. We are looking for the steady-state distribution, when

d n_{k} / d t = 0

. With

C_{i j} = 1

and

ϵ = 1

the above equations read for a steady state,

\begin{matrix} 0 = - n_{1} \sum_{i \geq 1} n_{i} + \frac{λ}{2} \sum_{i, j \geq 2} (i + j) n_{i} n_{j} \\ + λ n_{1} \sum_{i \geq 2} i n_{i} + λ_{e} l \sum_{i \geq 1} \sum_{j \geq l + 2} n_{i} n_{j}, \end{matrix}

[S34]

\begin{matrix} 0 = \frac{1}{2} \sum_{i + j = k} n_{i} n_{j} - (1 + λ) \sum_{i \geq 1} n_{i} n_{k} \\ + λ_{e} \sum_{i \geq 1} n_{i} n_{k + l} - θ (k - l - 2) λ_{e} \sum_{i \geq 1} n_{i} n_{k}, \end{matrix}

[S35]

where the Heaviside step function

θ (k) = 1

for

k \geq 0

and

θ (k) = 0

for

k < 0

.

Now we apply the generation function technique, with a slightly different definition of this function,

N (z) = \sum_{k \geq 1} n_{k} z^{k + l}

, so that the concentrations

n_{k}

are the coefficients of

N (z)

at

z^{k + l}

. It is straightforward to show that the generation function satisfies the quadratic equation

N {(z)}^{2} - 2 A (z) N (z) + B (z) = 0

with

A (z) = (1 + (1 + α) λ) N z^{l} - α λ N

[S36]

B (z) = 2 N ((1 + λ) n_{1} z^{2 l + 1} - α λ (1 - z^{l}) G (z)),

[S37]

where

G (z) = \sum_{i = 1}^{l + 1} n_{i} z^{i}

.

Because λ is small, we analyze the expansion of

N (z)

in terms of λ. We use this expansion for the total number of aggregates,

N = \sum_{i \geq 1} n_{i}

, and for the concentrations

N = N^{(0)} + λ N^{(1)} + λ^{2} N^{(2)} + λ^{3} N^{(3)} + \dots

[S38]

n_{i} = n_{i}^{(0)} + λ n_{i}^{(1)} + λ^{2} n_{i}^{(2)} + λ^{3} n_{i}^{(3)} + \dots .

[S39]

Substituting the above expansion for

n_{i}

into Eqs. S34 and S35 and summing these equations up, we find

N^{(0)} = n_{1}^{(0)} = 0

and

N^{(1)} = 2, N^{(2)} = - 4 + 2 α K, N^{(3)} = 8 - 4 α K

[S40]

n_{1}^{(1)} = 1, n_{1}^{(2)} = - 1 + α K, n_{1}^{(3)} = 1 - α K

[S41]

with

K = l \sum_{j \geq l + 2} n_{j}^{(1)}

.

Now we can analyze

N (z) = A (z) - \sqrt{A^{2} (z) - B (z)}

[S42]

in different orders of λ. Keeping only terms of the order of λ we find

N (z) = λ z^{l} (N^{(1)} - N^{(1)} \sqrt{1 - \frac{2 z n_{1}^{(1)}}{N^{(1)}}}) .

With

n_{1}^{(1)} = 1

and

N^{(1)} = 2

we obtain

N (z) = λ \sum_{k = 1}^{\infty} \frac{Γ (k - 1 / 2)}{\sqrt{π} Γ (k + 1)} z^{k + l},

[S43]

which implies, that in the first order in λ the concentrations of the aggregates

n_{k}

behave for

k ≫ 1

as

n_{k} ≃ \frac{λ}{\sqrt{π}} k^{- 3 / 2};

that is, we obtain exactly the same power-law behavior as previously, for the case of solely disruptive collision without erosion. To find the exponential cutoff for this power-law dependence, we need to consider next-order terms in λ.

To do this, we focus on the part of

N (z)

given by Eq. S42, corresponding to the terms

\sim z^{m}

with

m ≫ l, 1

. Because

A (z)

contains only terms up to

z^{l + 1}

, we analyze the behavior of

\sqrt{A^{2} (z) - B (z)} = \sqrt{f (z)}

.

Obviously,

f (z)

is a polynomial of the degree

2 l + 1

and may be written as

f (z) = (z_{1} - z) (z_{2} - z) \dots (z_{2 l + 1} - z)

, where

z_{1} \leq z_{2} \leq \dots z_{2 l + 1}

are the roots of

f (z)

. We are looking for the expansion of

\sqrt{f (z)}

in term of z in the vicinity of

z = 0

. Moreover, we are interested in the expansion coefficients for

z^{m}

with

m ≫ 1

. For

m ≫ 1

the expansion coefficients at

z^{m}

of

\sqrt{f (z)}

coincide with the expansion coefficients of the function

\sqrt{| z_{1} f^{'} (z_{1}) |} \sqrt{1 - z / z_{1}}

, which may be written as (69)

\sqrt{| z_{1} f^{'} (z_{1}) |} \sqrt{1 - \frac{z}{z_{1}}} = \sqrt{\frac{| z_{1} f^{'} (z_{1}) |}{4 π}} \sum_{k = 1}^{\infty} \frac{Γ (k - 1 / 2)}{Γ (k + 1)} {(\frac{z}{z_{1}})}^{k} .

[S44]

As follows from Eq. S42, the closest to the

z = 0

root of

f (z)

is equal to 1 (

z_{1} = 1

) in the first order in λ. Therefore, we need to find the next-order corrections to

z_{1}

in powers of λ:

z_{1} = 1 + λ ω_{1} + λ^{2} ω_{2} .

Substituting into

f (z) = A^{2} (z) - B (z) = 0

the expansions for N and

n_{i}

, given by Eqs. S40 and S41, and the above expansion for

z_{1}

, one can find the coefficients

ω_{1}

and

ω_{2}

,

ω_{1} = 0, ω_{2} = 1,

which means that

z_{1} = 1 + λ^{2} + \dots

. Therefore, the high-order terms of

N (z)

that contain

z^{k}

with

k ≫ l, 1

behave as

\frac{Γ (k - 1 / 2)}{Γ (k + 1)} {(\frac{z}{z_{1}})}^{k} \sim k^{- 3 / 2} e^{- k \log z_{1}} \sim k^{- 3 / 2} e^{- λ^{2} k} .

Such dependence of the concentrations

n_{k}

on k coincides with the one for the model of disruptive collisions discussed in the main text. This is confirmed by the numerical solution of the respective rate equations (Fig. S1).

Fig. S1.

Steady-state particle size distribution for completely disruptive collisions in the presence of collisions with an erosion. The solid lines depict the size distribution for the systems with the following parameters: $λ = 0.05$ , $λ_{e} = 0.05$ , $l = 9$ (black); $λ = 0.05$ , $λ_{e} = 0.08$ , $l = 3$ (red); and $λ = 0.01$ , $λ_{e} = 0.05$ , $l = 3$ (blue). The dashed lines indicate the respective size distribution for purely disruptive collisions for $λ = 0.05$ (black) and $λ = 0.01$ (blue). All curves correspond to the case of constant kinetic coefficients. As can be seen, the presence of the erosive collisions has practically negligible impact on the steady-state size distribution.

The case of

C_{i j} = {(i j)}^{μ}

may be analyzed analogously; in this case one obtains the exponent

- (μ + 3 / 2)

instead of

- 3 / 2

, that is, again the same behavior as for purely disruptive collisions. Hence we conclude that the presence of the collisions with an erosion does not change qualitatively the size distribution of particles in planetary rings.

Efficient Numerical Solution of the Kinetic Equations

We use Euler’s method for the numerical solution of the kinetic equations. This method is rather suitable for the case of interest, because we search stationary, continuous, and smooth solutions. The first problem in the numerical analysis of the infinite number of rate equations is the conservation of mass. Indeed, in any real simulation one can handle only a finite number of equations, say N, which describe evolution of particles of size

1, 2, \dots, N

(a particle of size k has mass

m_{k} = m_{1} k

). These equations have both aggregation and fragmentation terms. In particular, they have a term that describes aggregation of particles of size

i < N

and

j < N

, resulting in an aggregate of size

i + j > N

. Because the system of N equations does not account for particles larger than N, such processes would effectively lead to the leak of particles’ mass and violation of the mass conservation. To preserve the mass conservation we assume that all collisions of particles of mass i and j are fragmentative if

i + j \geq N

. We have checked that this assumption does not lead to any noticeable distortion of the numerical solution of the rate equations

n_{k}

, if k is smaller than some fraction of N.

Another problem is to handle efficiently a large number of equations, say up to

\sim 10^{6}

. One possible way is an application of the coarse graining, that is, grouping concentrations

n_{k}

–

n_{k + l}

into coarse-grained variables

{\tilde{n}}_{K}

with increasing l as k grows. In the case of interest, however, we have a drastic variation of the functional dependence of

n_{k} (k)

, which changes from a power-law to the exponential decay. This hinders an effective application of the coarse graining and we need to keep explicitly all individual concentrations. Hence we have to work with a large number of equations, which is computationally costly. To speed up the computations we have developed a recursive procedure.

In the case of fragmentation into monomers the system of kinetic equations has the following form:

\begin{array}{l} \frac{d n_{1}}{d t} = - n_{1} \sum_{j = 1}^{N} n_{j} + λ (1 - n_{1}) \sum_{j = 1}^{N} n_{j} \dots \\ \frac{d n_{k}}{d t} = \frac{1}{2} \sum_{i + j = k} C_{i, j} n_{i} n_{j} - (1 + λ) \sum_{i = 1}^{N} C_{i, k} n_{i} n_{k} \\ \frac{d n_{k + 1}}{d t} = \frac{1}{2} \sum_{i + j = k + 1} C_{i, j} n_{i} n_{j} - (1 + λ) \sum_{i = 1}^{N} C_{i, k + 1} n_{i} n_{k + 1} \dots \end{array}

[S45]

Taking into account that we search for a stationary solution,

d n_{k + 1} / d t = 0

, we obtain for the number density

n_{k + 1}

,

\frac{1}{2 (1 + λ)} \sum_{i + j = k + 1} C_{i, j} n_{i} n_{j} - \sum_{i = 1}^{N} C_{i, k + 1} n_{i} n_{k + 1} = 0 .

[S46]

The first sum in Eq. S46 contains only

n_{i}

with

i \leq k

, whereas we write the second sum as

\begin{matrix} \sum_{i = 1}^{N} C_{i, k + 1} n_{i} n_{k + 1} = n_{k + 1} \sum_{i = 1}^{k} C_{i, k + 1} n_{i} + C_{k + 1, k + 1} n_{k + 1}^{2} \\ + n_{k + 1} \sum_{i = k + 1}^{N} C_{i, k + 1} n_{i} . \end{matrix}

[S47]

Now we use the properties of the kinetic kernel

C_{i j}

and the steady-state distribution

n_{k} = n_{k} (m_{k})

, which we assume to be a decreasing function of k. Namely, we assume that the coefficients

C_{i j}

increase with i and j at a smaller rate than the rate at which

n_{k}

decreases with k. That is, we assume that for

k ≫ 1

the following condition holds true:

\sum_{i = 1}^{k} C_{i, k + 1} n_{i} ≫ \sum_{i = k + 1}^{N} C_{i, k + 1} n_{i} .

[S48]

This allows us to neglect the last sum in Eq. S47 and obtain the quadratic equation for

n_{k + 1}

:

C_{k + 1, k + 1} n_{k + 1}^{2} + n_{k + 1} \sum_{i = 1}^{k} C_{i, k + 1} n_{i} - \sum_{i + j = k + 1} \frac{C_{i, j} n_{i} n_{j}}{2 (1 + λ)} = 0 .

[S49]

Solving the above equation and choosing the positive root, we arrive at the recurrent relation for the concentrations

n_{k}

:

\begin{matrix} n_{k + 1} = \frac{\sqrt{D} - \sum_{i = 1}^{k} C_{i, k + 1} n_{i}}{2 C_{k + 1, k + 1}} \\ D = \frac{2 C_{k + 1, k + 1}}{(1 + λ)} \sum_{i + j = k + 1} C_{i, j} n_{i} n_{j} + {(\sum_{i = 1}^{k} C_{i, k + 1} n_{i})}^{2} . \end{matrix}

[S50]

Using the recurrent relation [S50], one can significantly accelerate computations. This may be done as follows: First, one solves explicitly the system of rate equations [S46] for

k ≪ N

, choosing the value of k to fulfill the condition [S48]. Then the concentrations

n_{i}

with

k < i \leq N

may be straightforwardly obtained from the recurrence [S50]. Performing a numerical solution of the rate equations with different kernels directly and with the use of the recurrence [S50], we proved the efficiency and accuracy of the above accelerating approach.

Numerical Calculation of the Distribution of Fragments

In the numerical solution we calculate

x_{k} (i)

, using the mass conservation and taking into account the discreteness of the distribution of particles:

k x_{k} (i) = B (i) \int_{k - 1 / 2}^{k + 1 / 2} k_{1} k_{1}^{- α} d k_{1} .

[S51]

Computing the integral, we find,

if

α \neq 2

,

x_{k} (i) = \frac{B (i)}{k} \frac{1}{2 - α} [{(k + \frac{1}{2})}^{2 - α} - {(k - \frac{1}{2})}^{2 - α}];

[S52]

if

α = 2

,

x_{k} (i) = \frac{B (i)}{k} [\ln (k + \frac{1}{2}) - \ln (k - \frac{1}{2})] .

[S53]

Note that

x_{k} (i) \to B (i) k^{- α}

for

k ≫ 1

, when the discreteness of the system becomes insignificant.

Here

B (i)

represents a normalization constant, which can be computed from

B (i) \int_{1 / 2}^{i - 1 / 2} k k^{- α} d k = i .

[S54]

Thus, we obtain,

if

α \neq 2

,

B (i) = \frac{i (2 - α)}{{(i - 1 / 2)}^{2 - α} - {(1 / 2)}^{2 - α}},

[S55]

so that

B (i) = 2^{2 - α} (α - 2) i

for

i ≫ 1

;

if

α = 2

,

B (i) = \frac{i}{ln (i - 1 / 2) - ln (1 / 2)} .

[S56]

From the last equations we see that the correction factor β, introduced above, reads, e.g., for the case of

α \neq 2

,

β = 2^{2 - α}

. In the case of exponential distribution we get analogously

\begin{matrix} k x_{k} (i) = B (i) \int_{k - 1 / 2}^{k + 1 / 2} k_{1} \exp (- k_{1}) d k_{1} \\ = B (i) [(k + \frac{1}{2}) e^{1 / 2 - k} - (k + \frac{3}{2}) e^{- 1 / 2 - k}] . \end{matrix}

[S57]

Here

B (i) = i / I_{0}

, with

I_{0} (i) = \int_{1 / 2}^{i - 1 / 2} k e^{- k} d k = (\frac{3}{2}) e^{- 1 / 2} - (i + \frac{1}{2}) e^{1 / 2 - i} .

Conclusion and Outlook

We have developed a kinetic model for the particle size distribution in a dense planetary ring and showed that the steady-state distribution emerges from the dynamic balance between aggregation and fragmentation processes. The model quantitatively explains properties of the particle size distribution of Saturn’s rings inferred from observations. It naturally leads to a power-law size distribution with an exponential cutoff (Eq. 19). Interestingly, the exponent

q = 2.5 + 3 μ

is universal, for a specific class of models we have investigated in detail. That is, q does not depend on details of the collisional fragmentation mechanism, provided the size distribution of debris, emerging in an impact, is steep enough; collisions with erosion do not alter q as well. The exponent q is a sum of two parts: The main part, 2.5, corresponds to the “basic” case when the collision frequency does not depend on particle size (

μ = 0

); such slope is generic for a steady size distribution, stemming from the aggregation–fragmentation balance in binary collisions. The additional part,

3 μ

, varying in the interval

0.25 \leq 3 μ \leq 1

, characterizes size dependence of the collision frequency. The latter is determined by the particles’ diameters and the mean square velocities

〈 v_{k}^{2} 〉

of their random motion. We have obtained analytical solutions for the limiting cases of energy equipartition,

(1 / 2) m_{k} 〈 v_{k}^{2} 〉 = const

, (

3 μ = 0.25

) and of equal velocity dispersion for all species

〈 v_{k}^{2} 〉 = const

, (

3 μ = 1

). These give the limiting slopes of

q = 2.75

and

q = 3.5

. Physically, we expect that an intermediate dependence between these two limiting cases may follow from a better understanding of the behaviors of threshold energies. This would imply a power-law size distribution with an exponent in the range

2.75 \leq q \leq 3.5

.

Observed variations of spectral properties of ring particles (39, 40) may indicate differences in the surface properties and, thus, in their elasticity and sticking efficiency. This implies differences in the velocity dispersion

〈 v_{k}^{2} 〉

and its dependence on k, resulting in different values of the exponent q. Moreover, variations in particle sizes among different parts of Saturn’s ring system have been inferred from Cassini data (16, 41). For our model, a different average particle size, or monomer size, implies different values of

E_{frag}

and

E_{agg}

as well as different values of the upper cutoff radius

R_{c}

.

Our results essentially depend on three basic assumptions: (i) Ring particles are aggregates composed from primary grains that are kept together by adhesive (or gravitational) forces; (ii) the aggregate sizes change due to binary collisions, which are aggregative, bouncing, or disruptive (including collisions with erosion); and (iii) the collision rates and type of impacts are determined by sizes and velocities of colliding particles. We stress that the power-law distribution with a cutoff is a direct mathematical consequence of the above assumptions only; that is, there is no need to suppose a power-law distribution and search for an additional mechanism for a cutoff as in previous semiquantitative approaches (9).

The agreement between observations and predictions of our model for the size distribution indicates that dense planetary rings are indeed mainly composed of aggregates similar to the dynamic ephemeral bodies suggested three decades ago (6, 7). This means that the (ice) aggregates constituting the cosmic disks permanently change their mass due to collision-caused aggregation and fragmentation, whereas their distribution of sizes remains stationary. Thus, our results provide another (quantitative) proof that the particle size distribution of Saturn’s rings is not primordial. The same size distributions are expected for other collision-dominated rings, such as the rings of Uranus (42, 43), Chariklo (44, 45), and Chiron (46, 47) and possibly around extrasolar objects (48–50).

The predictive power of the kinetic model further emphasizes the role of adhesive contact forces between constituents that dominate for aggregate sizes up to the observed cutoff radius of

R_{c} \sim 5 - 10 m

. The model does not describe the largest constituents in the rings, with sizes beyond

R_{c}

. These are the propeller moonlets in the A and B rings of Saturn, which may be remnants of a catastrophic disruption (12, 20). We also do not discuss the nature of the smallest constituents, that is, of the primary grains. These particles are probably themselves composed of still smaller entities and correspond to the least-size free particles observed in the rings (13).

Recently, cohesion was studied for dense planetary rings in terms of N-body simulations (51, 52). This model is similar to ours in that the authors use critical velocities for merging and fragmentation whereas we use threshold energies

E_{agg}

and

E_{frag}

; both criteria are based on the cohesion model. In these simulations a power-law distribution for the aggregates size

F (R) \sim R^{- q}

was obtained with slopes

2.75 \leq q \leq 3

for reasonable values of the cohesive parameter, consistent with our theoretical result. Moreover, the critical velocities for merging and fragmentation differ in most of the simulations by a factor of 2, which is in reasonable agreement with our model, where we estimated

E_{frag}

to be roughly twice

E_{agg}

(SI Text).

Materials and Methods

Boltzmann Equation.

The general equations [2] for the concentrations

n_{k}

have been derived from the Boltzmann kinetic equation. Here we consider a simplified case of a force-free and spatially uniform system. It is possible to take into account the effects of nonhomogeneity, as it is observed in self-gravity wakes, and gravitational interactions between particles. These, however, do not alter the form of resulting rate equations [2], which may be then formulated for the space-averaged values (SI Text).

Let

f_{i} \equiv f (m_{i}, v_{i}, t)

be the mass-velocity distribution function that gives the number density of particles of mass

m_{i}

with the velocity

v_{i}

at time t. In the homogeneous setting, the distribution function evolves according to the Boltzmann equation,

\frac{\partial}{\partial t} f_{k} (v_{k}, t) = I_{k}^{b} + I_{k}^{heat} + I_{k}^{agg} + I_{k}^{frag},

[21]

where the right-hand side accounts for particle collisions. The first term

I_{k}^{b}

accounts for bouncing collisions of particles (24); the second term

I_{k}^{heat}

describes the viscous heating caused by the Keplerian shear (53); and the terms

I_{k}^{agg}

and

I_{k}^{frag}

account, respectively, for the aggregative and disruptive impacts (explicit expressions for these terms are given in SI Text). To derive the rate equations [2] for the concentrations of the species,

n_{k} (t) = \int d v_{k} f_{k} (v_{k}, t)

, one needs to integrate Eq. 21 over

v_{k}

. Assuming that all species have a Maxwell velocity distribution function with average velocity

〈 v_{k} 〉 = 0

and velocity dispersion

〈 v_{k}^{2} 〉

, we obtain the rate equations [2] and the rate coefficients [1] (details in SI Text).

Generating Function Techniques.

To solve Eq. 6 we use the generating function

N (z) = \sum_{k \geq 1} n_{k} z^{k}

, which allows us to convert these equations into the single algebraic equation

N {(z)}^{2} - 2 (1 + λ) N N (z) + 2 (1 + λ) N n_{1} z = 0 .

[22]

Its solution reads

N (z) = (1 + λ) N [1 - \sqrt{1 - \frac{2 n_{1}}{(1 + λ) N} z}] .

[23]

Expanding

N (z)

, we arrive at the distribution [7].

Efficient Numerical Algorithm.

The numerical solution of Smoluchowski-type equations [2] is challenging as one has to solve infinitely many coupled nonlinear equations. We developed an efficient and fast numerical algorithm dealing with a large number of such equations. The application of our algorithm requires that the condition

\sum_{i = 1}^{k} C_{i, k + 1} n_{i} ≫ \sum_{i = k + 1}^{N} C_{i, k + 1} n_{i}

[24]

is obeyed for

k ≫ 1

and

N ≫ 1

, where

n_{i}

are the steady-state concentrations. For the case of interest this condition is fulfilled. Our algorithm first solves a relatively small set (

\sim 100

) of equations using the standard technique and then obtains other concentrations of a much larger set (

\sim 10,000

), using an iterative procedure (SI Text).

Acknowledgments

We thank Larry Esposito, Heikki Salo, Martin Seiß, and Miodrag Sremčević for fruitful discussions. Numerical calculations were performed using the Chebyshev supercomputer of Moscow State University. This work was supported by the Deutsches Zentrum für Luft und Raumfahrt, the Deutsche Forschungsgemeinschaft, and the Russian Foundation for Basic Research (project 12-02-31351). The authors also acknowledge partial support through the Marie Curie Action “International Research Staff Exchange Scheme” Dynamics and Cooperative Phenomena in Complex Physical and Biological Media (project 269139).

Supporting Information

Supporting Information (PDF)

Supporting Information

Download
288.81 KB

References

1

JN Cuzzi, RH Durisen, Bombardment of planetary rings by meteoroids - General formulation and effects of Oort Cloud projectiles. Icarus 84, 467–501 (1990).

Crossref

Google Scholar

2

JN Cuzzi, PR Estrada, Compositional evolution of Saturn’s rings due to meteoroid bombardment. Icarus 132, 1–35 (1998).

Crossref

Google Scholar

3

MS Tiscareno, et al., Observations of ejecta clouds produced by impacts onto Saturn’s rings. Science 340, 460–464 (2013).

Crossref

PubMed

Google Scholar

4

J Cuzzi, et al., An evolving view of Saturn’s dynamic rings. Science 327, 1470–1475 (2010).

Crossref

PubMed

Google Scholar

5

AW Harris, Collisional breakup of particles in a planetary ring. Icarus 24, 190–192 (1975).

Crossref

Google Scholar

6

DR Davis, SJ Weidenschilling, CR Chapman, R Greenberg, Saturn ring particles as dynamic ephemeral bodies. Science 224, 744–747 (1984).

Crossref

PubMed

Google Scholar

7

SJ Weidenschilling, CR Chapman, DR Davis, R Greenberg, Ring particles - Collisional interactions and physical nature. Planetary Rings, eds Greenberg R, Brahic A (University of Arizona Press, Tucson, AZ), pp. 367–415. (1984).

Google Scholar

8

NN Gorkavyi, AM Fridman, The rings of Uranus as resonances with unseen satellites. Sov Astron Lett 11, 302–303 (1985).

Google Scholar

9

PY Longaretti, Saturn’s main ring particle size distribution: An analytic approach. Icarus 81, 51–73 (1989).

Crossref

Google Scholar

10

RM Canup, LW Esposito, Accretion in the Roche zone: Coexistence of rings and ring moons. Icarus 113, 331–352 (1995).

Crossref

Google Scholar

11

F Spahn, N Albers, M Sremcevic, C Thornton, Kinetic description of coagulation and fragmentation in dilute granular particle ensembles. Europhys Lett 67, 545–551 (2004).

Crossref

Google Scholar

12

L Esposito, Planetary Rings: A Post-Equinox View, Cambridge Planetary Science (Cambridge Univ Press, Cambridge, UK), 2nd Ed. (2014).

Google Scholar

13

A Bodrova, J Schmidt, F Spahn, N Brilliantov, Adhesion and collisional release of particles in dense planetary rings. Icarus 218, 60–68 (2012).

Crossref

Google Scholar

14

EA Marouf, GL Tyler, HA Zebker, RA Simpson, VR Eshleman, Particle size distributions in Saturn’s rings from Voyager 1 radio occultation. Icarus 54, 189–211 (1983).

Crossref

Google Scholar

15

HA Zebker, EA Marouf, GL Tyler, Saturn’s rings - Particle size distributions for thin layer model. Icarus 64, 531–548 (1985).

Crossref

Google Scholar

16

J Cuzzi, et al., Ring particle composition and size distribution. Saturn from Cassini-Huygens, eds MK Dougherty, LW Esposito, SM Krimigis (Springer, Dordrecht, The Netherlands), pp. 459–509 (2009).

Crossref

Google Scholar

17

RG French, PD Nicholson, Saturn’s Rings II. Particle sizes inferred from stellar occultation data. Icarus 145, 502–523 (2000).

Crossref

Google Scholar

18

HA Zebker, GL Tyler, EA Marouf, On obtaining the forward phase functions of Saturn ring features from radio occultation observations. Icarus 56, 209–228 (1983).

Crossref

Google Scholar

19

MS Tiscareno, et al., 100-metre-diameter moonlets in Saturn’s A ring from observations of ‘propeller’ structures. Nature 440, 648–650 (2006).

Crossref

PubMed

Google Scholar

20

M Sremcevic, et al., A belt of moonlets in Saturn’s A ring. Nature 449, 1019–1021 (2007).

Crossref

PubMed

Google Scholar

21

MS Tiscareno, JA Burns, MM Hedman, CC Porco, The population of propellers in Saturn’s A ring. Astron J 135, 1083–1091 (2008).

Crossref

Google Scholar

22

AHF Guimaraes, et al., Aggregates in the strength and gravity regime: Particles sizes in Saturn’s rings. Icarus 220, 660–678 (2012).

Crossref

Google Scholar

23

N Borderies, P Goldreich, S Tremaine, A granular flow model for dense planetary rings. Icarus 63, 406–420 (1985).

Crossref

Google Scholar

24

NV Brilliantov, T Pöschel Kinetic Theory of Granular Gases (Oxford Univ Press, Oxford, 2004).

Crossref

Google Scholar

25

S Araki, S Tremaine, The dynamics of dense particle disks. Icarus 65, 83–109 (1986).

Crossref

Google Scholar

26

S Araki, The dynamics of particle disks. II. Effects of spin degrees of freedom. Icarus 76, 182–198 (1988).

Crossref

Google Scholar

27

S Araki, The dynamics of particle disks III. Dense and spinning particle disks. Icarus 90, 139–171 (1991).

Crossref

Google Scholar

28

C Dominik, AGGM Tielens, The physics of dust coagulation and the structure of dust aggregates in space. Astrophys J 480, 647–673 (1997).

Crossref

Google Scholar

29

K Wada, Collisional growth conditions for dust aggregates. Astrophys J 702, 1490–1501 (2009).

Crossref

Google Scholar

30

JA Astrom, Statistical models of brittle fragmentation. Adv Phys 55, 247–278 (2006).

Crossref

Google Scholar

31

C Guettler, J Blum, A Zsom, CW Ormel, CP Dullemond, The outcome of protoplanetary dust growth: Pebbles, boulders, or planetesimals? 1. Mapping the zoo of laboratory collision experiments. Astron Astrophys 513:A56. (2010).

Google Scholar

32

RC Srivastava, A simple model of particle coalescence and breakup. J Atmos Sci 39, 1317–1322 (1982).

Crossref

Google Scholar

33

H Salo, Numerical simulations of dense collisional systems: II. Extended distribution of particle size. Icarus 96, 85–106 (1992).

Crossref

Google Scholar

34

A Bodrova, D Levchenko, N Brilliantov, Universality of temperature distribution in granular gas mixtures with a steep particle size distribution. Europhys Lett 106, 14001 (2014).

Crossref

Google Scholar

35

F Leyvraz, Scaling theory and exactly solved models in the kinetics of irreversible aggregation. Phys Rep 383, 95–212 (2003).

Crossref

Google Scholar

36

PL Krapivsky, A Redner, E Ben-Naim A Kinetic View of Statistical Physics (Cambridge Univ Press, Cambridge, UK, 2010).

Crossref

Google Scholar

37

R Schrapler, J Blum, The physics of protoplanetesimal dust agglomerates. VI. Erosion of large aggregates as a source of micrometer-sized particles. Astrophys J 734, 108 (2011).

Crossref

Google Scholar

38

PL Krapivsky, E Ben-Naim, Shattering transitions in collision-induced fragmentation. Phys Rev E 68, 021102 (2003).

Crossref

Google Scholar

39

PD Nicholson, et al., A close look at Saturn’s rings with Cassini VIMS. Icarus 193, 182–212 (2008).

Crossref

Google Scholar

40

G Filacchione, et al., Saturn’s icy satellites and rings investigated by Cassini-VIMS. III. Radial compositional variability. Icarus 220, 1064–1096 (2012).

Crossref

Google Scholar

41

JE Colwell, J Cooney, LW Esposito, M Sremcevic, Saturn’s rings particle and clump sizes from Cassini UVIS occultation statistics. AGU Fall Meeting Abstracts 1. Available at abstractsearch.agu.org/meetings/2013/FM/P21E-01.html. Accessed July 2, 2015. (2013).

Google Scholar

42

JL Elliot, PD Nicholson, The rings of Uranus. Planetary Rings, eds R Greenberg, A Brahic (University of Arizona Press, Tucson, AZ), pp. 25–72 (1984).

Google Scholar

43

JL Elliot, RG French, KJ Meech, JH Elias, Structure of the Uranian rings. I - Square-well model and particle-size constraints. Astron J 89, 1587–1603 (1984).

Crossref

Google Scholar

44

F Braga-Ribas, et al., A ring system detected around the Centaur(10199) Chariklo. Nature 508, 72–75 (2014).

Crossref

PubMed

Google Scholar

45

R Duffard, et al., Photometric and spectroscopic evidence for a dense ring system around Centaur Chariklo. A&A 568, A79 (2014).

Crossref

Google Scholar

46

JL Ortiz, et al., Possible ring material around Centaur (2060) Chiron. A&A 576, A18 (2015).

Crossref

Google Scholar

47

JD Ruprecht, et al., 29 November 2011 stellar occultation by 2060 Chiron: Symmetric jet-like features. Icarus 252, 271–276 (2015).

Crossref

Google Scholar

48

Y Ohta, A Taruya, Y Suto, Predicting photometric and spectroscopic signatures of rings around transiting extrasolar planets. Astrophys J 690, 1–12 (2009).

Crossref

Google Scholar

49

EE Mamajek, et al., Planetary construction zones in occultation: Discovery of an extrasolar ring system transiting a young Sun-like star and future prospects for detecting eclipses by circumsecondary and circumplanetary disks. Astron J 143, 72 (2012).

Crossref

Google Scholar

50

MA Kenworthy, et al., Mass and period limits on the ringed companion transiting the young star J1407. Monthly Notices RAS 446, 411–427 (2015).

Crossref

Google Scholar

51

RP Perrine, DC Richardson, DJ Scheeres, A numerical model of cohesion in planetary rings. Icarus 212, 719–735 (2011).

Crossref

Google Scholar

52

RP Perrine, DC Richardson, N-body simulations of cohesion in dense planetary rings: A study of cohesion parameters. Icarus 219, 515–533 (2012).

Crossref

Google Scholar

53

J Schmidt, K Ohtsuki, N Rappaport, H Salo, F Spahn, Dynamics of Saturn's dense rings. Saturn from Cassini-Huygens, eds MK Dougherty, LW Esposito, SM Krimigis (Springer, Dordrecht, The Netherlands), pp. 413–458 (2009).

Crossref

Google Scholar

54

V Garzo, JW Dufty, Dense fluid transport for inelastic hard spheres. Phys Rev E 59, 5895–5911 (1999).

Crossref

Google Scholar

55

NV Brilliantov, F Spahn, Dust coagulation in equilibrium molecular gas. Math Comput Simul 72, 93 (2006).

Crossref

Google Scholar

56

V Garzo, CM Hrenya, JW Dufty, Enskog theory for polydisperse granular mixtures. ii. Sonine polynomial approximation. Phys Rev E 76, 031304 (2007).

Crossref

Google Scholar

57

J Piasecki, E Trizac, M Droz, Dynamics of ballistic annihilation. Phys Rev E 66, 066111 (2002).

Crossref

Google Scholar

58

JE Colwell, et al., The structure of Saturn's rings. Saturn from Cassini-Huygens, eds MK Dougherty, LW Esposito, SM Krimigis (Springer, Dordrecht, The Netherlands), pp. 375–412 (2009).

Crossref

Google Scholar

59

H Salo, Gravitational wakes in Saturn’s rings. Nature 359, 619–621 (1992).

Crossref

Google Scholar

60

JE Colwell, LW Esposito, M Sremcevic, Self-gravity wakes in Saturn’s A ring measured by stellar occultations from Cassini. Geophys Res Lett 33, L07201 (2006).

Crossref

Google Scholar

61

M Hedman, et al., Self-gravity wake structures in Saturn’s A ring revealed by Cassini VIMS. Astron J 133, 2624–2629 (2007).

Crossref

Google Scholar

62

RG French, H Salo, CA McGhee, LH Dones, HST observations of azimuthal asymmetry in Saturn’s rings. Icarus 189, 493–522 (2007).

Crossref

Google Scholar

63

A Toomre, On the gravitational stability of a disk of stars. Astrophys J 139, 1217–1238 (1964).

Crossref

Google Scholar

64

P Résibois, M De Leener Classical Kinetic Theory of Fluids (Wiley, New York, 1977).

Google Scholar

65

EA Peters, M Kollmann, TM Barenbrug, AP Philipse, Caging of a d-dimensional sphere and its relevance for the random dense sphere packing. Phys Rev E 63, 021404 (2001).

Crossref

Google Scholar

66

AP Hatzes, F Bridges, DNC Lin, S Sachtjen, Coagulation of particles in Saturn’s rings: Measurements of the cohesive force of water frost. Icarus 89, 113–121 (1991).

Crossref

Google Scholar

67

Bridges FG, Supulver KD, Lin DNC, Knight K, Zafra MD (1996) Energy loss and sticking mechanisms in particle aggregation in planetesimal formation. Icarus 123(2):422–435.

Google Scholar

68

NV Brilliantov, N Albers, F Spahn, T Pöschel, Collision dynamics of granular particles with adhesion. Phys Rev E 76, 051302 (2007).

Crossref

Google Scholar

69

P Flajolet, R Sedgewick Analytic Combinatorics (Cambridge Univ Press, New York, 2009).

Crossref

Google Scholar

Information & Authors

Information

Published in

Proceedings of the National Academy of Sciences

Vol. 112 | No. 31
August 4, 2015

PubMed: 26183228

Classifications

Submission history

Published online: July 16, 2015

Published in issue: August 4, 2015

Keywords

Acknowledgments

We thank Larry Esposito, Heikki Salo, Martin Seiß, and Miodrag Sremčević for fruitful discussions. Numerical calculations were performed using the Chebyshev supercomputer of Moscow State University. This work was supported by the Deutsches Zentrum für Luft und Raumfahrt, the Deutsche Forschungsgemeinschaft, and the Russian Foundation for Basic Research (project 12-02-31351). The authors also acknowledge partial support through the Marie Curie Action “International Research Staff Exchange Scheme” Dynamics and Cooperative Phenomena in Complex Physical and Biological Media (project 269139).

Notes

This article is a PNAS Direct Submission.

Authors

Affiliations

Nikolai Brilliantov

Department of Mathematics, University of Leicester, Leicester, LE1 7RH, United Kingdom;

View all articles by this author

P. L. Krapivsky

Department of Physics, Boston University, Boston, MA 02215;

View all articles by this author

Anna Bodrova

Faculty of Physics, Moscow State University, Moscow, 119991, Russia;

Institute of Physics and Astronomy, University of Potsdam, 14476 Potsdam, Germany;

View all articles by this author

Frank Spahn

Institute of Physics and Astronomy, University of Potsdam, 14476 Potsdam, Germany;

View all articles by this author

Hisao Hayakawa

Yukawa Institute for Theoretical Physics, Kyoto University, Kyoto 606-8502, Japan;

View all articles by this author

Vladimir Stadnichuk

Faculty of Physics, Moscow State University, Moscow, 119991, Russia;

View all articles by this author

Jürgen Schmidt¹ [email protected]

Astronomy and Space Physics, University of Oulu, PL 3000, 90014, Oulu, Finland

View all articles by this author

Notes

1

To whom correspondence should be addressed. Email: [email protected].

Author contributions: N.B. designed research; N.B., P.L.K., A.B., H.H., V.S., and J.S. performed research; N.B., P.L.K., A.B., H.H., V.S., and J.S. contributed new reagents/analytic tools; N.B., P.L.K., A.B., V.S., and J.S. analyzed data; and N.B., P.L.K., F.S., and J.S. wrote the paper.

Competing Interests

The authors declare no conflict of interest.

Metrics & Citations

Metrics

Note: The article usage is presented with a three- to four-day delay and will update daily once available. Due to ths delay, usage data will not appear immediately following publication. Citation information is sourced from Crossref Cited-by service.

Altmetrics

Citations

Cite this article

N. Brilliantov,
P.L. Krapivsky,
A. Bodrova,
F. Spahn,
H. Hayakawa,
V. Stadnichuk,
& J. Schmidt,

Size distribution of particles in Saturn’s rings from aggregation and fragmentation, Proc. Natl. Acad. Sci. U.S.A. 112 (31) 9536-9541, https://doi.org/10.1073/pnas.1503957112 (2015).

Copied!

Copying failed.

Export the article citation data by selecting a format from the list below and clicking Export.

Cited by

View Options

View options

PDF format

Download this article as a PDF file

DOWNLOAD PDF

Login options

Check if you have access through your login credentials or your institution to get full access on this article.

Personal login Institutional Login

Recommend to a librarian

Recommend PNAS to a Librarian

Save for later

Purchase options

Purchase this article to access the full text.

Single Article Purchase

Size distribution of particles in Saturn’s rings from aggregation and fragmentation

Proceedings of the National Academy of Sciences

Vol. 112
No. 31
pp. 9491-E4339

Restore content access

Restore content access for purchases made as a guest

Featured Topics

Articles By Topic

Featured Topics

Articles By Topic

Featured Topic

Articles By Topic

Significance

Abstract

Sign up for PNAS alerts.

Results and Discussion

Model.

Decomposition into Monomers.

Constant rate coefficients.

Size-dependent rate coefficients.

Power-Law Collisional Decomposition and Erosion.

Power-law decomposition.

Collisions with erosion.

Universality of the steady-state distribution.

Size Distribution of the Ring Particles.

Collision Integrals

Maxwell Approximation for the Velocity Distribution Functions

Derivation of the Rate Coefficients

Self-Gravity Wakes and Averaged Kinetic Equations

Estimates of the Characteristic Energies and the Aggregates Strength

Estimate of the Fragmentation Energy and of the Aggregates Strength.

Estimate of the Aggregation Energy.

Universality of Particles Size Distribution for Steep Distribution of Debris

Particle Size Distribution in the Presence of Collisions with Erosion

Efficient Numerical Solution of the Kinetic Equations

Numerical Calculation of the Distribution of Fragments

Conclusion and Outlook

Materials and Methods

Boltzmann Equation.

Generating Function Techniques.

Efficient Numerical Algorithm.

Acknowledgments

Supporting Information

References

Information

Published in

Classifications

Submission history

Keywords

Acknowledgments

Notes

Authors

Affiliations

Notes

Competing Interests

Metrics

Altmetrics

Citations

Cite this article

Cited by

View options

PDF format

Login options

Recommend to a librarian

Purchase options

Restore content access

Share

Share article link

Share on social media

Further reading in this issue

A role for the mevalonate pathway in early plant symbiotic signaling

Force generation by the growth of amyloid aggregates

IL-15 receptor α signaling constrains the development of IL-17–producing γδ T cells

Republicans are flagged more often than Democrats for sharing misinformation on X’s Community Notes

Integer partitions detect the primes

Bodily maps of emotions

Sign up for thePNAS Highlights newsletter