{'fraction_non_alphanumeric': 0.047396814523034764, 'fraction_numerical': 0.014942627162185305, 'mean_word_length': 4.309615821777677, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Anomalous aggregation regimes of temperature-dependent Smoluchowski equations

Temperature-dependent Smoluchowski equations describe the ballistic agglomeration. In contrast to the standard Smoluchowski equations for the evolution of cluster densities with constant rate coefficients, the temperature-dependent equations describe both -- the evolution of the densities as well as cluster temperatures, which determine the aggregation rates. To solve these equations, we develop a novel Monte Carlo technique based on the low-rank approximation for the aggregation kernel. Using this highly effective approach, we perform a comprehensive study of the phase diagram of the system and reveal a few surprising regimes, including permanent temperature growth and"density separation", with a large gap in the size distribution for middle-size clusters. We perform classification of the aggregation kernels for the temperature-dependent equations and conjecture the lack of gelation. The results of our scaling analysis agree well with the simulation data.

Temperature-dependent Smoluchowski equations describe the ballistic agglomeration. In contrast to the standard Smoluchowski equations for the evolution of cluster densities with constant rate coefficients, the temperature-dependent equations describe both -the evolution of the densities as well as cluster temperatures, which determine the aggregation rates. To solve these equations, we develop a novel Monte Carlo technique based on the low-rank approximation for the aggregation kernel. Using this highly effective approach, we perform a comprehensive study of the phase diagram of the system and reveal a few surprising regimes, including permanent temperature growth and "density separation", with a large gap in the size distribution for middle-size clusters. We perform classification of the aggregation kernels for the temperature-dependent equations and conjecture the lack of gelation. The results of our scaling analysis agree well with the simulation data.
Introduction. Aggregation processes are very ubiquitous in nature at different time and space scales, e.g. [1][2][3][4][5][6][7][8][9][10]. The classical tool to describe the aggregation kinetics is the celebrated Smoluchowski equations [6], which deal with the density of aggregates n k (t). Here the subscript k specifies the size of the aggregate, comprised of k monomers -the elementary units; these equations read for 1 ≤ k ≤ ∞ [7,8]: The kinetic coefficients C ij quantify the reaction rates between the aggregates of size i and j. They may either follow from a microscopic model or be constructed as empirical expressions [7,8]. Importantly, the classical Smoluchowski theory treats these coefficients as timeindependent so that the infinite set (1) forms a closed system of equations. It has been recently shown that for the ballistic agglomeration, the rate coefficients are time-dependent, C ij = C ij (t), as they are functions of the energy density (temperature T ) of the system [10,11] or of partial energy densities associated with the aggregates of size i (partial temperatures T i ) [12]. In the course of time, these quantities vary. Hence to make the system of equations closed, one needs to supplement Eqs. (1) for densities by the set of equations for temperatures T k (t) [12]: Here θ k = T k /m k and m k = m 1 k is the mass of aggregates of size k, with the diameter σ k = σ 1 k 1/3 . Eqs. (1) and (2) form a closed set; all of the rate coefficients C ij , B ij and D ij depend on the partial temperatures T i (t) and T j (t), see the Supplementary Material (SM). Moreover a microscopic analysis of particles collisions shows that besides of temperatures, the reaction rates sensitively depend on the interaction potential between particles, which may be put, for a wide class of interactions, into the form: where the constant a specifies the interaction energy, while λ 1 and λ 2 quantify the dependence of W ij on the size of particles i and j. For instance, λ 1 = λ 2 = 4/3 corresponds to the adhesive surface interactions, λ 1 = λ 2 = 3 stands for the dipole-dipole interactions and λ 1 = 3, λ 2 = 1 refers to the gravitational or Coulomb interaction, when the particles charges scale as their masses [12]. Still, the aggregation rate is determined not directly by W ij , but by its ratio to the characteristic kinetic energy of the colliding particles, that is, by the dimensionless quantity, where µ ij = m i m j /(m i + m j ) and ε is the restitution coefficient; it quantifies the dissipative losses at particles collisions. Hence C ij (t) = C ij (q ij (t)) and similarly B ij and D ij , see SM for details. Temperature-dependent Monte Carlo. The classical Smoluchowski system (1) may be analytically solved only for a few kernels C ij [7,8,13]. Generally, however, it requires a numerical analysis, e.g. [14][15][16][17][18][19]. Still, even the numerical solution is rather challenging, as the system of equations is infinite. The solution of the complete set (1) and (2) brings further complications. To address this problem, we develop a novel temperature-dependent Monte Carlo (MC) method, which is extremely efficient. It allows to investigate the behavior of huge systems in a wide range of parameters, which was not possible with the previous methods [20][21][22][23][24]. The main idea of the method is to exploit the low-rank approximation for the kinetic kernels, which has been successfully applied to solve classical Smoluchowski equations [25]. We adopt this approach to the MC scheme, with the extension for the temperature dependence, that is, for Eqs. (1) and (2). The method explicitly recalculates temperatures for arXiv:2110.00465v1 [cond-mat.stat-mech] 1 Oct 2021 each aggregate size without generating particle velocities. For our low-rank MC simulations only O (r log M ) operations are needed for each collision. For the system rank r ≡ 3, it allows performing ∼ 10 6 collisions every second without any use of parallel computation. Here M is the maximum cluster mass, and r is the rank of the kernel approximation. The implementation and detail of the new low-rank MC approach to our problem is discussed in SM.
Numerically-obtained phase diagram. We observe that the system obeying the temperature-dependent Smoluchiowski equations demonstrates extremely rich behavior: various temperature and aggregation regimes, including regimes of temporal and permanent temperature growth, the so-called "separation" regime, as well as classical aggregation with cooling. We vary two main parameters: Λ = λ 1 = λ 2 , which specifies, how the kinetic rates depend on the aggregates' size and q = q 11 (0) = a T1(0) 1 ε 2 2 Λ , which quantifies the (initial) ratio of the potential and kinetic energy of monomers. In simulations we use mono-disperse initial conditions (only monomers are available at t = 0) with the initial dimensionless density n 1 (0) = 0.0955 and temperature T 1 (0) = 1. We also use ε = 0.99, m 1 = 1 and σ 1 = 1.  Fig. 3.
The phase diagram in Fig. 1 illustrates the areas in the parametric space (q, Λ) corresponding to different evolution regimes. The most surprising is the aggregation "with separation" and aggregation with permanent temperature growth. In the former case, the cluster size distribution demonstrates an impressive density gap between small and large aggregates. That is, the density of intermediate-size clusters can be several orders of magnitude lower than that of small and large clusters, see Fig.  2. This may be explained by very large reaction rates for these clusters. In the latter case, the aggregation takes place in such a way that the rate of energy loss due to the agglomeration is lower than the aggregation rate. This results in the increasing energy per particle, that is, in growing temperature. Next, one can classify the regimes by evolution of temperature -it may follow the Haff's law with a continuous decay of temperature, or may alter the temperature regime, from the decay to growth and then back to decay; finally, a permanent growth from the very beginning is also possible.  Initial behavior. We start from classifying the initial behavior. For the qualitative analysis we use the average temperature T and ignore the difference of partial temperatures, T i = T and θ i = T /i (recall that m i = i). Then the densities n k follow Eqs. (1), while the equation for the average temperature reads, Here n(t) =  coefficients simplify to (see [12]): Initially monomers strongly dominate, that is, n ≈ n 1 , which yields from Eqs. (1) and (5) From the last equations and Eqs. (6) for C 11 and P 11 follows, where we take into account that in this caseṅ(0) Ṫ (0). This regime corresponds to the Haff's law of nonaggregative cooling, where c is a constant [12].
(b) For 1 − ε 2 q 1 we find: The solution to the last equation reads, where we also use Eqs. (7) and (6). The transition from the initial decay of temperature, as in (9), to the initial growth, as in (10), happens whenṪ (0) = 0, yielding for ε = 0.99, in agreement with the simulations, Fig. 1.
(c) Large q (low initial temperature), e −q (1 + q)) 1, give rise to the equation, Tṅ/n, and T = n 1/3 . This, together with Eqs. (7) and (6) leads to the solution: Unfortunately, it is not possible to observe this behavior for a sufficiently long time. Nevertheless, Eq. (12) shows that the temperature starts to decrease immediately. The value of q, demarcating the initial increase and decrease of temperature follows fromṪ (0) = 0, applied to Eq. (8): The result reads, again, in agreement with the simulation data, see Fig. 1. Scaling regimes and "separation". After a transient period when large clusters emerge, an aggregating system often enters a scaling regime, where the cluster size distribution may be described by a scaling function, n k (t) = s −2 Φ(k/s). Here the typical cluster size s = M 2 (t)/M 1 (M l = k n k k l ) is large, s 1. It scales as a power-law, s ∼ t z , for t → ∞ [7,8]. For temperaturedependent aggregation an additional temperature scaling T ∼ t −β may emerge. The scaling exponents z and β are determined by the characteristic values q ij for i, j ∼ s and Λ.
For small q ss 1 (high temperature) the rates C ij are homogeneous functions of i, j; the expansion of these coefficients (6) yields, for b 1. Similar expressions may be found for the rates P ij , yielding the following exponents, reported in [12]: This surprising behavior with the increasing average temperature T has been first reported in Ref. [12] as an intermediate regime for a finite-size system. Here we confirm it for different values of Λ (see Fig. 3) and demonstrate that it is permanent and stable in the thermodynamic limit, in agreement with the findings of [26]. The permanent growth of T is observed for Λ ≤ 5/4, see Fig.  1.
When Λ > 5/4, which corresponds to µ > 1, the classical criterion for gelation is fulfilled [7,27]. In this case, all mass is accumulated in a single gigantic gel, adsorbing all particles; the cluster densities become zero. For our system the situation is more involved -for very large clusters the condition q ij 1 is violated and converts into the opposite one, q ij 1, that is f ij 0, see Eqs. (6). The latter implies a purely ballistic agglomeration: with µ b = 1/6 < 1, which is strictly non-gelling. Hence the appearance of the pronounced gap in the cluster size distribution for the middle-size aggregates (we call it "separation") is the result of the "quasi-gelation" process. For small clusters, the scaling analysis is not applicable. The lack of "true"gelling in these systems is argued below. Large q ss (low temperature) with e −qss (1 + q ss ) 1 imply aggregation with cooling [12]: In simulations we observe, however, different exponent, T ∼ t −0.22 , which is caused by the failure of the approximation T i = T in this limit, see Fig. 3 (e), (f).
The lack of gelation. Qualitatively, the lack of gelation for temperature-dependent Smuluchowski equations follows from the conversion of the aggregation rates C ij for very large clusters into the ballistic, non-gelling case. These arguments exploit, however, the gelling criterion for the classical case. Hence it is desirable to prove the lack of gelation directly for the complete set of equations, (1) and (2). To this end, we need to prove that the moments of the density distribution would always be bounded [27]. Indeed, when the gelation happens, the second moment M 2 (t) = k n k k 2 becomes infinite [27]. To estimate M 2 (t) we also need the second moment for temperatures M 2,θ (t) = k n k θ k k 2 . Multiplying Eqs. (1) and (2) with k 2 and summing them over all k we arrive at (see SM for the derivation): Here E ij is the symmetric positive function, which describes the rate of energy exchange in bouncing collisions (see SM). For non-aggregating granular mixtures it was theoretically and numerically shown that T i grows with the size of the aggregates i, for all distribution with a density dominance of monomers [28]; we apply this for the studied aggregating systems. In all our simulations we always observed that either T i = T j , or T i > T j for i > j (except for a period of time with separation). Hence the conjecture T i > T j for i > j is justified. In this case the last term in Eqs. (19) is negative, yielding the result, that is, the moment M 2 (t) is bounded at any time instant and gelation is not possible. Classification of the rate kernels. In the classical temperature-independent theory, the aggregation rate kernels C ij are classified to be of three types, corresponding to qualitatively different behavior [7,27]. The temperature-dependent theory possesses an extended set of kernels: C ij , B ij and D ij . For practical reasons, it would be worth to have the according classification aligned with the classical one.
For small q 1 one can apply the expansion leading to Eqs. (14) and (15). From these equations follows that Class I kernel is realized for Λ > 3/4, which corresponds to ν < µ. This class is characterized by the power-law size dependence of the density of small clusters, Φ(x) ∼ x −(1+µ−Λ/2) for x 1 (Fig. 3a), where T (x) ∼ x Λ/3 for small x. The case Λ = 3/4 corresponding to ν = µ is classified as Class II kernels. It is characterized by similar power-law dependence for small clusters, Φ(x) ∼ x −τ , as for the Class I, but the exponent τ is not universal. Finally, Class III is realized for Λ < 3/4, that is for ν > µ. It is characterized by the exponential disappearance of small clusters, Φ(x) ∼ exp(−ax µ−ν−Λ/2 ) for x 1 (Fig. 3b) [7,27].
In the opposite case of large q 1 temperature always decreases, f ij 0 and kernels tend to the ballistic one, given by Eq. (17), with µ b = 1/6. This regime is stable and again corresponds to the Class III kernel. Note that in aggregation with separation small-size clusters eventually disappear, leading to the same scaling solution as for q 1. Therefore, all systems locating on the phase diagram outside permanently increasing temperature region always converges to the same Class III scaling with µ = µ b .
Conclusion. We investigate the aggregation kinetics for temperature-dependent Smoluchowski equations. In contrast to standard Smoluchowski equations, which describe the evolution of clusters densities of different sizes with time-independent rates, we consider two coupled sets of equations -one for the cluster densities and another for the partial temperatures of the aggregates, which define the reaction rates. For the numerical solution of these sets of equations, we develop a novel, highly efficient Monte Carlo approach based on the low-rank approximation of the reaction kernels. It allows for a fast solution of huge systems of equations. We explore the system's behaviour for a wide range of parameters and obtain a complete evolution phase diagram. It possesses several surprising regimes, including the regime of permanent temperature growth and density "separation". We make a classification of the aggregation kernels and conjecture a lack of gelation for temperature-dependent aggregation. The results of our scaling analysis are in good agreement with the simulation results.

SUPPLEMENTARY MATERIAL
Below we present the detail of the new Monte Carlo method for the solution of the temperature-dependent Smoluchowski equations. We also give the complete ex-pressions for the rate kernels and some derivation details. References to the equations from the main text are given in bold.

Kinetic rate coefficients for the temperature-dependent Smoluchowski equations
Here we write the temperature-dependent Smoluchowski equations (1) with the explicit indication of the temperature dependence of the rate coefficients C ij : Here T i and T j are temperatures of clusters of size i and j, which generally differ for clusters of different size. These temperatures obey the following equations, written for the reduced variables, θ i = T i /m i , where m i = im 1 is the mass of clusters of size i [12]: D kj n k n j , k = 1, ∞.
(22) Equations (1)-(2) have been obtained for the ballistic agglomeration of particles, interacting with the energy W ij , Eq. (3), which depends on the size of the particles and on the nature of the inter-particle forces (see the main text). The according kinetic coefficients read [12]: , W ij = a i 1/3 j 1/3 λ1 i 1/3 + j 1/3 λ2 , Under the approximation of equal temperatures, T i = T for all i (valid in some regimes), Eqs. (22) reduce to Eqs. (5)-(6).

Temperature-dependent Monte Carlo
In order to make our simulations equivalent to the solution of temperature-dependent Smoluchowski equations, we assume the distribution of speeds to be Maxwellian for each cluster size [12]. Naturally, we need a large number of clusters of each size to justify this assumption. But even if clusters of a particular size become scanty or disappear, this will not noticeably impact the overall solution, as the assumption remains true for the most of collisions.
Hence we can use for the Monte Carlo simulations the kernels C ij , B ij and D ij , as defined in (23), to determine the collision frequencies and the corresponding temperature variation.
The time step τ between collisions can be determined using equation (21). Note that the factor 1 − f ij in the kernel C ij defines the aggregation probability. Therefore, to determine the time between any (not only aggregative) collisions, we can use the kernelĈ ij = C ij / (1 − f ij ). Let us define the system volume as V = N/n = N (t = 0)/n(t = 0), where N is the total number of clusters and n is the total cluster density. If only one particle disappears during the time τ , then Hence, where N i is the number of particles of size i. Since Eq. (25) allows for aggregation of the particle with itself, we need to exclude such collisions, applying the rejection step with the probability 1/N i whenever i = j. The coefficient 1 2 in Eq. (25) prevents from a double counting of each colliding pair.
If we fix i and j we can calculate average time between collisions τ i,j for some fixed pair of sizes. Similarly, taking into account only aggregative collisions of particles of size i and j, we can calculate the average time between the according aggregation events, which we denote τ agg i,j : To determine the average change of temperatures during the aggregative collision of particles of size i and j, we can use kernels B ij and D ij . We find that Similarly we can find the set of equations for the temperature variation during the restitutive collisions using the other part D res ij = D ij − D agg ij of the kernel D ij : (29) Note that with the above procedure, we can compute the average temperature change directly. This is significantly more efficient than applying the procedure of randomly choosing cluster speed from a Maxwellian velocity distribution. In other words, we effectively apply the mean-field approximation for temperatures. Also, note that in this version of MC, we use effective virtual particles corresponding to an ensemble of real particles. Hence we can straightforwardly double the system size any number of times.
If we now apply Gillespie [20] or inverse [21] method, we need at least O(M ) operations after every collision, where M is the maximum cluster size. Indeed, when θ i changes, we need to recalculate i-th row and column of the collision kernelĈ. However, if we use a low-rank approximation of the coagulation kernel [25,26,29], only O(r) elements under this approximation change; here r is the rank of the system.
Next, we give a brief description of the algorithm in the context of temperature-dependent ballistic aggregation. The paper with the detailed description of the low-rank Monte Carlo algorithm for the general case is in preparation.
To approximate the collision frequenciesĈ ij we use the following matrix A of rank 3: Eq. (30) shows that we can pick the probabilities up from the matrix 2A instead ofĈ (we can replace A ji by A ij because of the collision symmetry). We use rejection sampling for the cases, where we have an overestimate, which happens with the probability This idea is similar to the majorant kernels approach exploited in Ref. [23]. After the multiplication of the matrix A by the vectors composed of n i and n j , we still get the sum of 3 rank-1 matrices. Indeed, if then A ij n i n j = A • nn T ij , where • denotes the element-wise product. It can be also defined as In order to keep track of the sums and quickly choose the sizes of the colliding particles, we construct a segment tree on each of the vectors u (k) and v (k) . Then the update of the full structure requires O(r log M ) operations and in case of ballistic agglomeration r = 3. In particular, segment trees contain the sums of all elements of u (k) and v (k) and thus we also know the sums of the u (k) v T (k) elements. To choose a pair of colliding particles, we firstly pick one of 3 rank-1 matrices with the probability, proportional to the sum of its elements (O(r) operations) then the element of u (k) with the probability, proportional to its value (O(log M )) and finally a column of v (k) (another O(log M ) operations). Therefore, the total complexity of the method is O(r log M ).
Derivation of the upper bounds M2(t) and M 2,θ (t) Here we derive the differential equations (19) for the second moment of number density M 2 = k n k k 2 and the second moment of n k θ k , equal to M 2,θ = k n k θ k k 2 .
First, we estimate the upper bound of the sum k n k T k = nT avg , where T avg is the average temperature. Multiplying the system (22) by k and summing over all k leads to The sum of expressions containing 1 − g ij is always negative. Bouncing term also cannot lead to increase of the kinetic energy, so it must be negative too. Therefore, the sum k n k T k always decreases and thus can be bounded as n(t)T avg (t) n(0)T (0) = const . Here we used the Cauchy-Schwarz inequality i n i f (i)g(i) i n i f (i) j n j g(j), along with Eq. (31) to limit the sum i n i T i and the equation for the total mass of the system i n i i = M 1 = const. Now, we find the equation for the moment M 2,θ (t).
Note, that the expression, containing 1 − g ij is always negative and so can be neglected. Next, in the positive term of the bouncing part, we can replace 1 + ε by 2, which will always increase the sum. Then