Increasing temperature of cooling granular gases.

The kinetic energy of a force-free granular gas decays monotonously due to inelastic collisions of the particles. For a homogeneous granular gas of identical particles, the corresponding decay of granular temperature is quantified by Haff's law. Here, we report that for a granular gas of aggregating particles, the granular temperature does not necessarily decay but may even increase. Surprisingly, the increase of temperature is accompanied by the continuous loss of total gas energy. This stunning effect arises from a subtle interplay between decaying kinetic energy and gradual reduction of the number of degrees of freedom associated with the particles' dynamics. We derive a set of kinetic equations of Smoluchowski type for the concentrations of aggregates of different sizes and their energies. We find scaling solutions to these equations and a condition for the aggregation mechanism predicting growth of temperature. Numerical direct simulation Monte Carlo results confirm the theoretical predictions.

Introduction

Dilute systems of inelastically colliding particles, also called granular gases, belong to the most intensively studied systems in the physics of granular matter[1,2]. Granular gases may be investigated using the tools of kinetic theory and hydrodynamics of molecular gases, with some modifications. These modifications are necessary to take into account the loss of mechanical energy due to particle collisions, quantified by the coefficient of restitution, defined as the ratio of post-collisional and pre-collisional relative normal velocity. Granular gases are inherently unstable as they tend to develop self-organized clusters[3,4] and other instabilities[5]. Although these instabilities inspired physicists to study granular gases, even the first stage in the evolution of a granular gas when it is homogeneous and isotropic reveals a rich phenomenology, which are fundamentally different from molecular gases. These may be exemplified by various types of correlations[4,6-9] and deviations from the Maxwell velocity distribution[10-16]. Based on an enormous amount of scientific literature in the past 25 years, by now the homogeneous state of granular gases is well understood[1,17-19].

An important feature of natural granular gases is the agglomeration of particles due to forces of different type. Such phenomena can be found, e.g., in form of soot-agglomeration in smoke gas[20], in aerosols[21,22], as well as in astrophysical systems such as dust in interplanetary space, in planetary rings consisting of aggregative ice particles[23] and protoplanetary disks[24-28]. While simple approaches assume that particle collision are always aggregative[28-32], more detailed models[27,28,33] take into account that only collisions at small impact rate lead to aggregation, while for larger impact velocity particles rebound or shatter[28,33]. Naturally, the mechanism of particle interactions determines conditions under which the particles aggregate, thus, a realistic model of aggregation should consider the physics of particle interactions, such as surface adhesion[27,28,33], or van der Waals interactions[34], electrostatics of charged or dipolar particles[35], or gravity[34].

Although it has been shown that the incorporation of aggregation changes the properties of granular gases considerably[28,30,31,33,35-38], by now the role of the physical aggregation mechanism was not considered. Moreover, while the size distribution of aggregates has been analyzed[27,28,33,35], a constant mean kinetic energy was assumed; a coupling between the evolution of the size distribution and that of the mean kinetic energy of the aggregates has not been investigated. It is however well known that the mean kinetic energy of the system is not constant but decreases with time[30,31,36-38], which strongly influences the kinetics of ballistic aggregation.

For an initially uniform gas of monomers, we analyze the evolution of the size distribution of aggregates and the average kinetic energy of these species for different realistic mechanisms of ballistic aggregation. From the Boltzmann equation for a multi-component granular gas[1], we derive a set of coupled equations of Smoluchowski type for the concentration of aggregates of different sizes and the corresponding average kinetic energies. We obtain different regimes of the system’s evolution characterized by either decreasing or increasing temperature. The increase of temperature is a surprising result given that all collisions of particles are dissipative.

Results

Kinetic theory

We consider a dilute and uniform gas of granular particles and introduce the mass–velocity distribution function, f ≡ f (m, v, t), which is the concentration of particles of i monomer masses at velocity v at time t. Its evolution obeys the Boltzmann equation[28,33]:where and are, respectively, the collision integrals for the aggregative and restitutive collisions. Let us introduce the aggregation model and explain the structure of Iagg and Ires. When non-adhesive granular particles of mass m and m collide, the energy of their relative motion, with v = v − v, is reduced by the factor ε2. The coefficient of restitution, ε ≤ 1, thus, characterizes the inelastic nature of particle collisions[1]. For adhesive particles considered here, we assume that the particles agglomerate if the post-collisional relative kinetic energy is smaller than the energy needed to overcome the energy barrier W due to attractive forces, that is, . The energy, W, depends on the size or the particles and the nature of the attractive forces. For instance, for adhesive forces, , where the constant A depends on material parameters and r are the particles’ radii[39]. Let a particle of mass m consist of i monomers of mass m0 and radius r0. Then an aggregate of mass m = m0i is of radius r ~ i1/3. The general formcharacterizes the dependence of W on the size of colliding particles, where a is a constant of the dimension of energy. The choice of λ1/2 describes a variety of attractive interactions. For example, λ1 = λ2 = 4/3 corresponds to the adhesive surface interactions introduced above; λ1 = 3, λ2 = 1 characterizes gravitational, or Coulomb interaction when the particles’ charges scale as their masses. Similarly, λ1 = λ2 = 3 stands for dipole–dipole interactions[35], etc. Note that in the aggregation condition, we do not consider explicitly the particle rotation; we assume that its impact on the agglomeration kinetics may be effectively accounted for by the factor a in Eq. (2).

For the case that particles aggregate, that is, ε2E ≤ W, the corresponding collision integral reads[28,33]Here Θ ≡ Θ(W − ε2E), with the Heaviside step-function, Θ(x), guarantees that the aggregation condition is fulfilled; m = m + m and mv = vm + mv, due to the mass and momentum conservation, σ = r + r and the inter-center unit vector e at the collision instant. The factors in the integrand in Eq. (3) have their usual meaning[1]: is the volume of the collision cylinder and Θ (−v ⋅ e) selects only approaching particles. The first sum in the right-hand side of Eq. (3) refers to collisions, where a particle of size k is formed from smaller particles of size i and j, while the second sum describes the collisions of k-particles with all other aggregates. In case of restitutive collisions, the collision integral has its usual form[1],with the additional factor Θ (ε2E − W) to exclude an aggregation at the impact. Here, , with and being the velocities of the inverse collision, that results with the post-collision velocities v and v[1].

Let us introduce the number density, n, of particles of size i, their partial temperature T, the total number density, N, and the average temperature, T[40]:We assume that the distribution function may be approximated by a Maxwellian,where is the mean thermal velocity of the particles of size i. Multiplying Eq. (1) with unity and with , using Eqs. (3) and (4) and integrating over v, we obtainwhere θ = T/m and the coefficients C, B, and D depend on σ, W, m, m, and θ, θ, see “Methods” section for details. The set of Eq. (7) for the concentrations of species (i.e., of the agglomerates of size k) and their average kinetic energy are the extended Smoluchowski-type equations: They describe the aggregation processes when different species have time-dependent individual temperatures. For non-aggregative particles, that is for W = 0, we obtain C = B = 0 and D = ξ(θ + θ)/(nθ), where ξ are the cooling coefficients[41,42]. The set of Eq. (7) then reduces to in agreement with refs.[41,42].

From numerical simulations (see below) we find that the partial temperatures evolve as T(t) = ϕT(t), where the constants ϕ weakly depend on i. For the qualitative analysis we assume T(t) ≈ T(t), where T(t) is the characteristic temperature of the system. Then Eq. (7) reduce towithand we abbreviate and μ = mm/(m + m). Below, we will consider some limiting cases.

Hot gas limit

For , that is, for all particles, the aggregation barrier is much smaller than the average kinetic energy of the particles. From Eq. (10) then follows C ~ (a/T)2 and P contains two terms, one ~(1 − ε2) and the other one ~(a/T)3. This entails different regimes of gas behavior:

Non-aggregative cooling occurs when (a/T) → 0, thus, C ≈ 0, Eq. (8) yields , which implies the equationwith Haff’s solution[43] for a non-aggregative cooling gas, T ~ (1 + t/τ0)−2, where .

Partial aggregation with cooling takes place for smaller temperature. In this case, we take into account aggregation but disregard the cooling due to aggregation. Then C ~ (a/T)2 and P ~ (1 − ε2); in more detail, and , with the dimensionless kernelsand constants b1 = k0(a/ε2)2, b2 = (2/3)k0(1 − ε2) and . The kernels are homogeneous functions of i and j:where μ = 1/6 and μ = 2(2λ1 − λ2)/3 + 1/6 = (2/3)Λ + 1/6. Here we introduce Λ = 2λ1 − λ2, where Λ/3 is the homogeneity degree of the attraction barrier (2).

For large time, t, Smoluchowski equations with homogeneous kernels give rise to scaling solutions[44-46], therefore, we seek solution to Eqs. (8) and (9) in the scaling formFrom Eq. (15) follows N(t) ~ t = −. Substituting Eq. (15) into Eqs. (8) and (9), using the homogeneity of kernels and approximating the discrete variables by continuous variables, , one obtains the equation for the scaling function Φ(x) (see more details in “Methods” section),where w is the separation constant[45,46]. Equation (16) is to be supplemented by the relations for the exponents z and β, which follow from the requirement that the power of t should be the same on both sides of Eqs. (8) and (9):Solving for z and β and using the expressions for μ and μ, we obtain for this regime,Note that the scaling analysis is valid if z > 0, therefore, 5 − Λ > 0. If this condition is violated a more subtle analysis is needed; this may indicate a gelation[44-46] or run-away growth[35].

Aggregation with temperature growth may take place for a molecular gas (ε = 1) or nearly elastic granular gas. In this case, and we obtain, , with given by Eq. (13), while P ~ (a/T)3. It may also be written as , where b3 = (2/9)k0(a/ε2)3 and readsthat is, is a homogeneous function with the homogeneity constant μ = Λ + 1/6.

The scaling analysis given above yields again Eq. (16) for the scaling function and the relations for the exponents z and β:From Eq. (20), we find

The above result is surprising: If Λ > 0, which corresponds to the increase of the attractive barrier with aggregates’ size, the exponent β is negative, that is, , and the characteristic temperature of an aggregating gas increases! This effect may be understood as follows: For an elastic or nearly elastic gas, the aggregation conditions (defined by the constants λ1 and λ2) may lead to a regime, when the total number of particles N decays faster than the total energy of the system. In this regime, the increase of the gas temperature occurs during its cooling, hence it behaves like a system with a negative-specific heat[47,48], see the discussion below.

Cold gas limit

Complete aggregation with cooling occurs in the limit , such that for all species. The aggregation barrier is much larger than the average kinetic energy of the particles, therefore, almost all the collisions lead to aggregation. Then from Eq. (10) follows, and wherethat is, both kernels and are homogeneous with the homogeneity constant μ = μ = 1/6. From Eq. (22) follows , then Eqs. (8) and (9) entail the equation:The scaling analysis of Eqs. (8) and (23) yields again Eq. (16) for Φ(x) and the relationswhich finally give the exponents,This is the aggregation regime of ordinary Smoluchowski type, when all collisions lead to agglomeration of the particles.

Numerical confirmation

To check the predictions of the previous section, we performed large-scale Monte Carlo simulations. The details of the numerical method are given in “Methods” section below. Figure 1 shows the evolution of the partial temperatures, T, of i-mers.

Fig. 1

Partial temperatures of i-mers. Evolution of the partial temperatures, T, of i-mers for a granular gas of 107 particles for λ1 = λ2 = 4/3, a = 0.1, and ε = 0.99. Initially, the gas of monomers has the temperature, T1(0) = 1. The dashed lines show the limiting cases of a hot gas (regime of non-aggregative cooling) when the temperature follows Haff’s law, T ~ t−2, and cold gas when almost all collisions are aggregative and T ~ t−1/3. Both cases are in agreement with the theory, see Eqs. (11) and (25)

While for an ordinary granular gas, temperature would monotonously decay due to Haff’s law, T(t) ~ t−2, in case of a gas of aggregative particles, temperature evolves non-monotonously by revealing a number of different regimes. We start with a homogeneous gas of monomers at temperature T1(0) = 1. Then initially, the condition is satisfied and the gas cools according to Haff’s law. This corresponds to the regime of the non-aggregative cooling for a hot gas described by Eq. (11). In the course of time, temperature drops down, the aggregation becomes non-negligible such that the evolution of temperature deviates from Haff’s law; the system now evolves according to the regime of partial aggregation with cooling. With further decrease of temperature and increase of the ratio (a/T), aggregation dominates over dissipative loss for the evolution of temperature. This corresponds to the regime of aggregation with temperature growth, where the characteristic temperature of the system increases as , see Eq. (21). Eventually, the latter regime terminates because the condition needed for this regime is violated with increasing temperature. The larger the size i of a particle, the larger the value of W and hence the large temperature T, when the change of the regimes takes place; this is clearly visible in Fig. 1. In the final state of its evolution, the gas approaches the limit of a cold gas, when almost all collisions are aggregative. In this regime, the temperature decays with time as T ~ t−1/3, in agreement with the theoretical prediction, Eq. (25).

Figure 2 shows the regime of aggregation with temperature growth in more detail. In agreement with the theoretical prediction, Eq. (21), the increase of temperature in this regime becomes steeper with increasing value of Λ, which is proportional to the homogeneity degree Λ/3 of the aggregation barrier.

Fig. 2

Evolution of temperature for different aggregation mechanisms. The rate of temperature growth in the regime of increasing temperature depends on the aggregation mechanism, quantified by the parameter Λ. This parameter characterizes the dependence of the aggregation barrier on the agglomerate size. a N = 107, ε = 0.99, a = 0.1, Λ = 0.4, thus β = −0.173, see Eq. (21); b same but Λ = 1.6, thus β = −0.941. With increasing value of Λ, the increase of temperature becomes steeper, in agreement with the theoretical predictions, , of Eq. (21)

It is also interesting to estimate the onset time, t, when the granular temperature starts to increase. This regime occurs in a hot granular gas when . If initially (a/T)3 is much smaller than (1 − ε2), the regime of non-aggregative cooling takes place and temperature decreases according to Haff’s law, T(t) ~ 1/(t/τ0)2. Here,  =  with η being the packing fraction of the gas (see the “Methods” section) and τ is the initial mean collision time. In the course of time, temperature decrease so that (a/T)3 grows. For the order of magnitude estimate of the onset time t one can use the condition (1 − ε2) ≈ (a/T(t))3. For the parameters ε = 0.99, a = 0.1, T(0) = 1, σ = 1, m = 1, and η = 0.05 used in the simulations depicted in Figs. 1 and 2, we obtain . This qualitatively agrees with the numerical result, of t ≈ 1000, which corresponds to approximately 600 collisions with the duration of , see Figs. 1 and 2.

Conditions for experimental observations

We estimate the range of parameters where the surprising phenomenon of increasing temperature of a cooling granular gas may be observed. As it follows from the above analysis, this is expected for a hot gas, when . Note, however, that while the condition allows a simplified analytical treatment, a milder condition (1 − ε2) < (a/T)3 ensures the increase of temperature as well.

Let us consider the case of adhesive interactions associated with the van der Waals forces[34,49]. In this case, the adhesion energy W defined in Eq. (2) reads[39],where γ is the adhesion coefficient, which is twice the surface free energy per unit area of a solid in vacuum, R = rr/(r + r), , with ν and Y being, respectively, Poisson ratio and Young modulus, and q0 = 1.457 is a numerical constant. With r = r0i1/3, where r0 is the monomer radius, we estimate,which implies the homogeneity exponents λ1 = λ2 = 4/3 and the constant a in Eq. (2):Since the homogeneity degree of W is positive, Λ/3 = (2λ1 − λ2)/3 = 4/9 > 0, the regime of aggregation with temperature growth is possible.

Using the data available in the literature for the restitution coefficient[50,51] and material parameters (γ, Y, ν, and the material density ρ) we obtain, that for a granular gas of very hard ceramic particles (Y = 370 × 109 GPa) of diameter σ = 1 mm moving at the characteristic velocity v0 = 0.16 ms−1 , the value of is indeed significantly smaller than . For the opposite case of rather soft acrylic particles (Y = 3 × 109 GPa), similar estimates indicate that such regime may be realized for a granular gas of particles of the same size but for the characteristic velocity v0 = 0.27 ms−1; in this case, and . Hence, we conclude that the regime of increasing temperature of a cooling granular gas may be observed for many realistic systems.

Let us estimate the onset time t for real systems discussed above, which are prepared at the hot gas conditions, and then freely evolve. Suppose that initially (a/T(0))3 = 0.001. Then using the estimates for t outlined above, we obtain that t = 17 s for the first system of (hard) ceramic particles of diameter σ = 1 mm; the gas is to be prepared with the packing fraction η = 0.05 and initial characteristic velocity v0(0) = 3 ms−1. Similarly, we obtain t = 72 s for the second system of (soft) acrylic particles of the same diameter, prepared at the same packing fraction with the characteristic velocity v0(0) = 6 ms−1.

Experimental investigations of cooling granular gases are not simple since to assure force-free conditions, the action of gravity has to be suppressed. Such experiments have been performed, either under true microgravity conditions aboard the space station[52], in parabolic flights[53,54], or in sounding rockets[55,56], or by means of magnetic levitation[57]. By now, however, the focus of these experiments was to check experimentally the cooling law theoretically predicted by Haff and others[43]. Hence, the experiments have been designed to assure purely repulsive interaction of the particles, that is, to keep attractive forces as small as possible. This might be the reason why the temporary increase of temperature reported in our paper was not mentioned in any of these studies. As it follows from the estimates given above, to ensure the emergence of the regime of interest, one needs to prepare a system at particular conditions.

Analogy with systems of negative-specific heat

It is worth to note that the behavior of an aggregating granular gas with increasing temperature may be interpreted as a manifestation of a negative-specific heat: The energy of the system decreases in dissipative collisions, while the granular temperature grows, that is, the smaller the energy of the gas, the larger its temperature. The negative-specific heat characterizes equilibrium systems with long-range (e.g., gravitational) interactions, see e.g., refs. [47,48]. Such analogy is very interesting and tempting. Still however, the direct application of this concept to non-equilibrium systems, addressed in our study, is to be justified by further analysis. In particular, one needs to investigate the opposite process, when temperature decreases with increasing system energy. This may be possible for an extended model, which includes shattering collisions and mechanisms of energy input into the system. Then, if the number of particles, emerging in shattering impacts, increases faster than energy of the system, thanks to the energy input, the temperature of the gas would decrease, manifesting a negative heat capacity. To apply the concept of negative heat capacity to aggregating gases, a detailed analysis of the thermodynamic additivity[47,48] in their cooling and heating states will also be needed. We leave these fascinating problems for future studies.

Discussion

By means of kinetic theory, we show that a force-free gas of aggregative particles behaves fundamentally different than an ordinary granular gas of purely repulsive particles. In both cases, the kinetic energy decreases monotonously due to dissipative collisions of the particles. However, while the temperature of ordinary granular gases decreases monotonously as well, according to Haff’s law, T(t) ~ t−2, the evolution of granular gases of aggregating particles demonstrates subsequent regimes of very different nature. These regimes are characterized by either decreasing or increasing temperature. The astonishing fact that the temperature of a gas of particles, which dissipate energy increases, may be understood as follows: although the total energy of the system decreases, the total number of particles diminishes faster due to agglomeration, yielding a boost of energy per particle. This surprising effect is possible if the aggregation barrier increases with the size of the aggregates (i.e., Λ > 0), and, as we have shown above, may be observed for many realistic systems. Interestingly, the increase of temperature with decaying energy corresponds to a negative-specific heat for equilibrium systems[47,48]; the application of this concept to non-equilibrium systems requires however further analysis.

Technically, we derived an extended set of Smoluchowski equations: equations for concentrations of different species, n, which correspond to standard Smoluchowski equations and equations for the partial temperatures of the species, T. Using the approximation T(t) ≈ T(t), we elaborate a scaling theory and reveal several regimes of temperature evolution, including the surprising regime of increasing temperature of a cooling gas. Numerical Monte Carlo simulations confirm the predictions of our theory.

Methods

Kinetic kernels of extended Smoluchowski equations

To derive the first set of rate Eq. (7) for n(t), we integrate the Boltzmann equation, Eq. (1), over v. Since , the left-hand side of the Boltzmann equation turns then into dn/dt and gives the rate of change of the concentrations n. The right-hand side gives the contributions to dn/dt from different parts of the collision integral. Since the restitutive (non-aggregating) collisions preserve the number of particles, we findTo compute , we assume that the velocity distribution function may be approximated by a Maxwellian, Eq. (6). For a non-aggregating granular gas, the velocity distribution function deviates from Maxwellian; the deviation may be expressed in terms of Sonine polynomials expansion[10-14,16]. These deviations are however small and, as it has been shown in refs. [41,42,58], may be neglected when the cooling coefficients ξ are computed. Therefore, the approximation (6) is justified.

The collision integral for aggregative collisions, given by Eq. (3) may be written as a sum of two parts,First, we compute ,where we use Eq. (6) for the velocity distribution function. 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 substitutewherethe above integral with respect to v, v, and e may be generally written aswhere , , and α = β = γ = 0, for the particular case of the integrals in Eq. (30). We also take into account that the Jacobian of transformation from (v, v) to (u, w) is equal to unity. Integration over u gives (π/a)3/2. Integration over the unit vector gives 4π and we are left with the integral over w. Integration over directions of the vector w gives π, so finally we need to compute the remaining integral:As the result we obtain:Similarly, one can find the integral in Eq. (33) for other values of α ≠ 0, β ≠ 0, and γ ≠ 0.

Let us now calculate . First, we notice that , since the other part of the integrand does not depend on . Then the remaining integration is exactly the same as for , which have been already performed, therefore we find:Hence, after some algebra, we obtain,where

To derive the second set of rate equations in the system, Eq. (7), for n(t)θ(t), we multiply the Boltzmann equation, Eq. (3) with and integrate over v. Since , the left-hand side of the equation turns into 3md/dt(nθ). In the right-hand side of this equation we again encounter the Gaussian integrals, as in Eq. (33), but with different values of α, β, and γ. These integrals may be computed exactly in the same way as it has been shown above for the particular case, α = β = γ = 0. Finally, we arrive at the set of equations for θ = T/m,with the coefficients:where Q has been defined above and Δ = m/(m + m).

For W = 0 we obtain C = B = 0 andwhere ξ are the cooling coefficients for a non-aggregating granular mixture[41,42,58]. Then Eq. (39) takes the form,in agreement with refs. [41,42,58].

If T = T for all i, we find and arrive at Eq. (11) for C. The coefficients P of Eq. (10) may be obtained noticing thatUsing then Eqs. (40) for B and D, we arrive for T = T at P, given by Eq. (10).

Scaling analysis of the extended Smoluchowski equations

Generally, it is not possible to solve analytically the infinite set of Eqs. (8) and (9), however, one can find a scaling solution, which reflects the most prominent features of the exact solution. To understand the nature of the scaling approach, we first consider a non-rigorous, qualitative analysis and then focus on its rigorous counterpart. Exploiting the basic scaling hypothesis[45,46] we assume that for large time and large k a scaling solution to Eqs. (8) and (9) is attained. Hence, we seek the solution to these equations in the formwhere s(t) is the characteristic aggregate size. It may be defined through the second moment as s = M2/M1, where and M1 = const. is the total mass which is conserved. Multiplying Eq. (8) with k2 we obtainwhere we use the condition of the mass conservation along with the result of the main text for the dimensionless kernel:Taking into account that the total number of aggregates and its average size are related as N ~ M1/s ~ 1/s, we recast Eq. (9) into the formwhere we again use the relation for the dimensionless kernel:  ~ . If we now assume the scaling behavior for temperature, T ~ t−, and substitute it into Eqs. (45) and (47), we find from the condition of consistency the relations (17).

Turn now to more rigorous analysis. As it has been shown above, in the limit of hot granular gas as well as in the cold gas limit, the kinetic kernels of Eq. (8) may be written as:where c, p, and ν, ν are constants and and are the dimensionless homogeneous functions of i and j with the homogeneity exponents μ and μ, respectively. We also assume that the evolution of temperature obeys a power law, T ~ t−, while the concentrations n obey the scaling form,Following the approach of ref.[45], we write the left-hand side of Eq. (9) as:where x = k/t. This is equal to the right-hand side, which we write, changing from the discrete variables i, j to the continuous ones, C → C(i, j) as:where we take into account that k = xt, j = yt and useand similar relations that hold true for homogeneous kernels. Hence, Eq. (9) transforms intoIf we divide both sides of the above equation by , we observe that one side of the equation depends only on t while the other one on x and y. This implies the following condition[45]where w is the separation constant. This condition yields Eq. (16) for Φ(x) and the relationFor the scaling solution (49), one can write for :Then one can write for Eq. (10):which gives another scaling relation:The scaling relations, Eqs. (55) and (58), correspond to the pairs of scaling relations, Eqs. (18), (20), and (24) for different values of ν and ν, that one has for different aggregation regimes.

Direct simulation Monte Carlo

To confirm the theoretical results, we perform numerical direct Monte Carlo simulations (DSMC). DSMC is a numerical technique used to directly solve the Boltzmann equation. It was first elaborated by Bird[59] for the simulation of molecular gases and later generalized and applied to the Enskog equation for dissipative granular gases, see e.g., refs. [60-64]. As there is an extended introductory literature on the application of DSMC to granular gases, e.g. refs. [65,66], here we describe only the details of our approach where we deviate from the standard DSMC.

Laboratory time for non-aggregating particles

For reasons which will be clear below, our DSMC algorithm does not rely on the laboratory (real) time but on the number of performed collisions. For the simple case of non-aggregating particles where the total number of particles is preserved, we compute the laboratory time à posteriori from the number of collisions via the collision frequency, ν:That is, during a short interval of time, Δt, where the temperature, T, does not change noticeably in a homogeneous system of N particles in volume V, there occur CΔ collisions. Here, we use the notation,where m and σ are respectively the mass and diameter of the particles. Then an interval of the laboratory time Δt may be computed through the number of collisions CΔ as:where we introduce the notation,with ρ being the density of the particles’ material. In our simulations we use σ = 1, m = 1, ρ = 6/π, T(0) = 1 and the packing fraction η = (π/6)σ3(N/V) = 0.05, which guarantees the accuracy of DSMC. The initial mean collision time in these units is τ = ν−1 = 1.48.

DSMC simulation of aggregative particles

When the particles of a dissipative gas agglomerate, obviously, their total number is not preserved. In our simulations, we start with typically 108 monomers and in the course of time larger and larger particles emerge, such that after some time not monomers, but much larger particles dominate. Simultaneously, the total number of particles decreases persistently. After a relatively short period of time it will not be possible to obtain reliable data due to the poor statistics. Therefore, we keep the number of particles in the system approximately constant by regularly expanding the system: we start the simulation with N monomers, but when arriving by agglomeration at particle number N/2, we duplicate all particles. That is, each particle is replaced by two particles of identical velocity and mass. Effectively, this operation corresponds to doubling the system volume, V, which defines a scaling variablewhere V0 is the initial system volume. After the first doubling, F = 2, after the second, F = 4, etc. Consequently, the data obtained as result of the simulation have to be re-scaled before applying Eqs. (61) and (62):Here, NMC is the number of particles, simulated at the current time corresponding to a certain temperature. Similarly, all other variables with index MC. This way, the number of particles in the simulation is always in the interval N = [N0/2 + 1, N0], where N0 is the initial number of particles. This expansion of the system size allows to obtain good statistical data independently of the aggregation process.

Laboratory time for agglomerating particles

The composition of a granular gas of aggregating particles changes in time. That is, at a given time, the system consists of N(1) monomers, N(2) dimers,…,N( i-mers. The abundance of monomers can be considered as a subsystem of particles, the same for dimers and, in general, for i-mers. For the case of agglomerating particles, we use the fact that in dilute systems, subsystems of particles behave independent of one another. Obviously, since these subsystems coexist, there are many ways to determine the laboratory time from the number of collisions via Eq. (61), where we use the data from subsystem i: N(, m(, T(. Consequently, from the simulation results we obtain many different laboratory times, which are, of course, theoretically identical, but not numerically due to the evolution of the system and the finite system size: near the beginning of the simulation, most of the particles are monomers, therefore, the statistics of the monomers delivers the most reliable results while the statistics based on 20-mers is unreliable since the abundance of 20-mers is yet small. For a given particle size, i, there is a certain interval of time when the abundance of i-mers is the largest subsystem from which we obtain the most reliable data. For earlier time, the statistics is poor since most of the particles are yet smaller than i, for later times, the statistics is poor as well since most of the i-mers vanished from the system by agglomeration. Consequently, in our simulation we determine the laboratory time from the most reliable species, I, consisting of the largest number of particles at the given time. Therefore, we compute

Data availability

All relevant data are available from the authors.

Peer Review File