Supercluster states and phase transitions in aggregation-fragmentation processes

We study the evolution of aggregates triggered by collisions with monomers that either lead to the attachment of monomers or the break-up of aggregates into constituting monomers. Depending on parameters quantifying addition and break-up rates, the system falls into a jammed or a steady state. Supercluster states (SCSs) are very peculiar non-extensive jammed states that also arise in some models. Fluctuations underlie the formation of the SCSs. Conventional tools, such as the van Kampen expansion, apply to small fluctuations. We go beyond the van Kampen expansion and determine a set of critical exponents quantifying SCSs. We observe continuous and discontinuous phase transitions between the states. Our theoretical predictions are in good agreement with numerical results.

Addition and disintegration mimic social phenomena, e.g., users joining forums which may eventually disintegrate (partially or completely).Aggregating and disintegrating objects in social networks may be also firms, enterprises, etc. [29][30][31].
The addition process is symbolically represented by the reaction scheme (see also Fig. 1) Thus an elementary object (a monomer denoted by M ) collides with another object (a monomer or a cluster) to form a cluster of larger mass.Here I k denotes a cluster composed of k monomers (k-mer), and A k is the merging rate.The Becker-Döring equations [32][33][34][35] and their continuum counterpart, Lifshitz-Slyozov-Wagner model [36,37], rely on aggregation with addition mechanism.
Fragmentation works concurrently with aggregation.In fragmentation, clusters break into smaller clusters.Fragmentation can occur spontaneously or via collisions.We consider fragmentation processes caused by collisions with monomers, i.e., dual to addition that is also caused by collisions with monomers.Schematically where ℓ i=1 j i = k.We thus tacitly assume that a monomer always remains separated after collision, mimicking the situation when an energetic monomer hits a cluster and retains its identity.The collision-induced disintegration processs (1) is described in the framework of Oort-Hulst-Safranov-Dubovski models [38][39][40][41].In complete disintegration, or shattering, clusters break into monomers [25,26,[42][43][44].Symbolically In this work, we explore systems with addition and fragmentation and report a rich set of behaviors, including continuous and discontinuous phase transitions, the formation of steady and jammed states, and the emergence of supercluster states (SCSs) which are peculiar jammed states.In the SCS, the cluster densities vanish in the thermodynamic limit.To shed light on the SCSs, we study finite systems and show a non-extensive nature of the SCSs manifested by sub-linear scaling of the total number of clusters with system size.To describe the SCSs, we develop a framework extending the van Kampen expansion applicable to extensive systems.
The SCSs have been detected in our Letter [45].Here we present a more detailed analysis of the SCSs.We consider a model with complete disintegration (shattering), and also with partial disintegration, and demonstrate the emergence of SCSs in both cases.Hence we conjecture that SCSs are generic for systems with addition and fragmentation.We also present a detailed analysis and classification of phase transitions that were only briefly mentioned in [45].
Our analytical treatment is focused on addition with shattering.Each such model is characterized by a set of rates A k and S k , and some models are analytically tractable as we demonstrate in Sects.II-IV.

II. ADDITION AND SHATTERING PROCESSES
The governing equations describing addition and shattering read where Eqs.(3b) apply to all k ≥ 2. The first and second terms on the right-hand side of Eq. (3a) describe the loss of monomers via addition due to monomer-monomer and monomer-cluster interactions.The last term in (3a) describes the gain of monomers due to shattering.Similarly, the last two terms on the right-hand of Eqs.(3b) represent the loss of k-mers due to addition and shattering, while the first term gives the gain due to addition.In writing Eqs.(3) we tacitly assume that the system is well-mixed, spatially homogeneous, and dilute.Even when these assumptions are satisfied, the description provided by Eqs.(3) is mean-field in nature, so it may be erroneous in low spatial dimensions [44,46].
For many practical applications the rates of addition and shattering depend algebraically on the cluster size: Here we set the amplitude of the addition rate to unity by using the appropriate time units; the constant λ quantifies the shattering intensity.The dependence A k ∼ k a is natural since the aggregation rate is often proportional to the surface area [35,[47][48][49].This additionally implies that a ≤ 1. Models with rates growing faster than mass, a > 1, are also ill-defined in the thermodynamic limit as an infinite cluster forms at time t = +0, see e.g.[50][51][52][53]).For networks, growth with exponent a > 0 reminds preferential attachment [29], and the behavior also drastically changes when a > 1, see [54].The exponent s is determined by the shattering mechanism and it usually satisfies the constraint s ≤ 1.By re-scaling densities, we set the mass density to unity if not stated otherwise: Using the modified time we linearize Eqs. ( 3) Hereinafter M ν = k≥1 k ν c k denotes the ν th moment.Equation (7a) is not closed.If the system of equations for c 1 , M a , M 1+s is closed, one can proceed analytically.The moment M ν evolves according to The term k≥1 (k + 1) ν k a c k on the right-hand side of (8) can be expressed through the moments only when ν is a non-negative integer.Thus closed equations for the moments emerge when a and 1 + s are non-negative integers.In the physically acceptable range a ≤ 1 and s ≤ 1, there are five possibilities: (a, s) = (1, 1), (a, s) = (0, 1), (a, s) = (1, 0), (a, s) = (0, 0) and (a, s) = (0, −1); these models admit analytical treatment.The one-parameter class of models with exponents (a, s) = (a, a − 1) is also partly tractable as we shall show below.
If not stated otherwise, we always consider the most natural mono-disperse initial condition Solving (9a)-(9b) subject to (10) gives Different behaviors emerge for λ < 1, λ = 1 and λ > 1.In the subcritical region, λ < 1, the monomer density vanishes at τ = τ max (λ).Setting c 1 (τ max ) = 0 in (11b) and solving for τ max yields implying that τ max (λ) is an increasing function of λ.The final modified time increases from τ max (0) = ln(2) to τ max (1) = ∞.The quantity τ max (λ c ) remained finite in models with (a, s) = (0, 0) and (a, s) = (1, 1) studied in [55], and also in the model with (a, s) = (0, 1).This is a mathematical reason for the peculiarity of the critical regime in the present model.Combining (11a) and (12) gives the total final density (valid for the mono-disperse initial condition): Equations (9c) can be solved recurrently starting with the monomer density (11b).One finds etc. Specifying to τ max given by Eq. ( 12) one establishes c k (∞).The results become more and more cumbersome, and we have not found a general compact formula valid for all k.We mention exact results only in the extreme case of λ = 0 when shattering is absent.This pure addition process was solved in [6].The densities are and In the other extreme, namely in the proximity of the critical point, 0 < 1 − λ ≪ 1, the final densities for k ≥ 2. Thus the final densities vanish linearly, c k (∞) ∼ 1 − λ.The total cluster density N ∞ also vanishes, but in a different manner, namely as √ 1 − λ in the λ ↑ 1 limit.This follows from Eq. (13).
The critical regime occurs at λ = λ c = 1.Solving Eqs.(9) subject to (10) yields Inverting the definition (6) we obtain allowing us to re-write (17) as The final densities c k (∞) vanish at λ c = 1, yet the mass density is conserved.The same happens in pure aggregation where c k (t) → 0 as t → ∞ yet the mass distribution widens, and the sum k≥1 kc k (t) = 1 remains constant.The distinction with pure aggregation become clear if we compare the final outcomes in a finite system.All mass is engulfed by a single cluster in an aggregation process.
In the addition and shattering process with rates A k = k and S k = 1, the final state is very different as we show below.Monomers also disappear, but the overall number of clusters diverges with total mass, albeit sub-linearly.We call such final outcome a supercluster state, as clusters are predominantly large.
In the supercritical regime (λ > 1), the cluster densities relax to the steady state following from Eqs. (9).Summarizing with τ max (λ) depending on the initial conditions (it is given by (12) for the mono-disperse initial condition).
In the vicinity of the critical point for k ≥ 2. Thus we have a continuous phase transition from a jammed state to a steady state occurring through the critical supercluster state, see Fig. 3.
B. The model with (a, s) = (0, −1) In this model the evolution is governed by The critical shattering parameter λ c = 1 demarcates different regimes (Fig. 4).In a subcritical regime, λ < 1, the system eventually arrives at a jammed state with vanishing monomer density, c 1 (τ max ) = 0, and cluster densities c k (τ max ), determined by initial conditions.
When λ = λ c = 1, one solves (22) iteratively to give etc. Equivalently etc.The asymptotic behavior for t ≫ 1 and fixed k reads All final densities vanish at the critical point (see Fig. 4).
To determine the asymptotic behavior of the total density N we return to (22b) with λ = 1.Treating k as a continuous variable we arrive at a wave equation which is solved to give The dominant contribution to the sum N = k≥1 c k is gathered near k * ≈ τ where c k (τ ) has a sharp maximum.
This observation allows us to compute N relying only on the sharpness of the maximum and mass conservation In the general case of arbitrary λ we apply the Laplace transform to Eq. (22b) and obtain which is iterated to find Applying the Laplace transform to k≥1 kc k (τ ) = 1 gives k≥1 k c k (p) = 1/p.Using Eq. ( 27) and which is the definition of the (ordinary) hypergeometric function [56] we obtain Substituting ( 28) into ( 27) we arrive at Using (29) we express the Laplace transform of the total cluster density through the ratio of hypergeometric functions: The τ → ∞ behaviors (corresponding to t → ∞) are encoded in the p → 0 behaviors of the Laplace transforms.To extract such behaviors we use the integral representations of the hypergeometric functions appearing in Eqs. ( 29)- (30): If λ > 2, these hypergeometric functions are regular at p = 0 (that is, at ϵ = 1) and equal to Thus all c k (p) have a simple pole, c k (p) → c k (∞)/p, indicating the existence of a steady state size when λ > 2.
Using (29) we find Summing c k (∞), or using (30), we find the total cluster density in the steady state (λ > 2): Albeit all cluster densities relax to stationary values when λ > 2, the moments M ν (τ ) with ν ≥ λ − 1 grow indefinitely.Indeed, Eq. (33) gives For instance, cluster densities are stationary in the range 2 < λ ≤ 3, but the second moment M 2 (∞) diverges, so M 2 (τ ) grows indefinitely.We now show that for τ ≫ 1.To establish the asymptotic behaviors (36) we rely on the Laplace transform of the second moment derived using Eqs.(29).As usual, the large τ behavior of M 2 (τ ) is encoded in the p → 0 behavior of M 2 (p).The hypergeometric function in the denominator in the right-hand side of Eq. ( 37) is regular at p = 0 [its value is given by (32b)].The hypergeometric function in the nominator is defined via [56] When p ↓ 0, we have ϵ ↑ 1 and observe that the sum is dominated by the large n behavior.Hence summation can be replaced by integration.Using the asymptotic Γ(1+λϵ+n) ≃ n −λϵ and ϵ n−1 ≃ e −pn , and setting ϵ = 1 in regular terms, we obtain Plugging (32b) and ( 38) into (37) we get from which we deduce an algebraic growth of the second moment given in (36) in the 2 < λ < 3 range.In the marginal case of λ = 3 we similarly derive from which we deduce a logarithmic growth given in (36).
The time dependence (36) suggests that when τ ≫ 1 the densities are stationary and given by Eq. ( 33) up to a crossover size k * ∼ τ , while for k > k * the densities quickly vanish.This implies that the moments diverge algebraically when ν > λ − 1, viz. as and logarithmically in the marginal case: In terms of the physical time we thus have when λ > 2 and t ≫ 1.
We now turn to the critical region 1 ≤ λ ≤ 2. Inside this region, F [2, 2; 2 + λϵ; ϵ] diverges at p → 0 implying that c k → 0 for τ → ∞, see Eq. ( 29).We conclude that the system possesses a critical interval with the lower λ c,l = 1 and upper λ c,up = 2 critical points.Again, vanishing cluster densities in conjunction with mass conservation indicate the formation of the supercluster state discussed below.Hence the final densities of clusters and monomers read: with δ k,1 ensuring that c 1 (∞) = 0 in the jammed regime.
The densities c k (τ max ) depend on the initial conditions.Thus the system undergoes a continuous phase transition from a jammed state into a supercluster state at λ c,l = 1, and a continuous phase transition from a supercluster state to a steady state at λ c,up = 2 (see Fig. 4).
Analyzing the singularity of F [2, 2; 2 + λϵ; ϵ] as p → 0, one can find the temporal relaxation of the cluster densities in the 1 < λ < 2 range.With ϵ = (1 + p) −1 and p → 0 we can replace ϵ by 1 in the right-hand side of Eq. (31b), apart from the denominator (1 and analyzing the integral, we find that its dominant part gathers in the region 1 Computing the integral we obtain Inserting ( 43) into ( 28) we find from which we extract the asymptotic The total cluster density follows from Eqs. (22a) and ( 45): Using we can re-express previous results via physical time, e.g.
Here we consider a class of models with A k = k a and S k = λk a−1 .The governing equations read In Sects.II A and II B we discussed two representatives of this class of models: (a, s) = (1, 0) and (a, s) = (0, −1).
FIG. 4. The model with (a, s) = (0, −1).A continuous phase transition from a jammed state at λ < 1 to the supercluster state for 1 < λ < 2 is followed by a continuous phase transition from the supercluster state to a steady state at λc,up = 2.The final densities in a jammed state depend on the initial conditions: Shown are results for the mono-disperse initial conditions (solid lines) and the monomer-dimer initial conditions (dashed lines) with c1(0) = 0.2 and c2(0) = 0.4.Curves: The analytical and numerical solution of the ODE.Dots: Monte Carlo results (in simulations, the total number of monomers was M = 10 6 ).
Since τ max = ∞ above the critical point, we can apply the Laplace transform for λ ≥ 1.Then from Eqs. (52b) we iteratively obtain Using again k≥1 k c k (p) = 1/p following from mass conservation, we find p c The Laplace transform c 1 (p) has a simple pole at p = 0 giving the steady-state density c 1 (∞) = 1/G(0, λ).Setting p = 0 on the right-hand side of Eq. ( 56) and massaging the sum we obtain The summand behaves as Γ(2 + λ) k −(λ+a−1) in the large k limit, so the sum on the right-hand side of (57) converges when λ > 2 − a.If λ ≤ 2 − a, the sum in (57) diverges yielding vanishing final densities.Hence λ c,up = 2 − a as stated in (53).
The above steady-state densities do not depend on the initial conditions.Summarizing A continuous phase transition from the jammed state to the supercluster state occurs at the lower critical point, λ c,l = 1.Then a continuous phase transition from the supercluster state to the steady state takes place at the upper critical point, λ c,up = 2 − a.
Consider the critical region 1 < λ < 2 − a. First, we re-write G(p, λ) given by Eq. ( 56) as Re-writing the product as and expanding in the p → 0 limit gives with the second asymptotic valid when k ≫ 1.We also use the asymptotic valid when k ≫ 1.In the p → 0 limit, the main contribution to C(p, λ) is gathered when k ≫ 1.This allows us to replace summation over k by integration and use (59a)-(59b).When a < 1 we get where Λ ≡ λ−1 1−a .This parameter varies in the range 0 < Λ < 1 in the critical region 1 < λ < 2 − a.
We have c 1 (p) = 1/pG(p, λ) and take its inverse Laplace transform to extract the large time asymptotic In terms of the physical time The total cluster density exhibits the same temporal behavior as the density of monomers.Asymptotically, For a = 0, we recover R(0, λ) = (λ + 1)/(λ − 1) in the critical region 1 < λ < 2.
The asymptotic behaviors ( 61)-( 63) are valid inside the critical region 1 < λ < 2 − a.More peculiar behaviors occur at the boundaries λ = 1 and λ = 2 − a.In Eqs. ( 60)-( 62), we have also assumed that 0 < a < 1.The behaviors at the boundaries again require more careful treatment.We have analyzed these behaviors: the model with a = 0 (Sec.II B) and the model with a = 1 (Sec.II A).
At the lower critical point, the ratio converges when a > 0 and diverges when a ≤ 0. The decay of the density of monomers at λ = 1 can be extracted from ( 62) by taking the λ → 1+0 limit.Using . This asymptotic disagrees with c 1 ∼ t −1 decay at λ = 1 and a = 0 due to already mentioned peculiarities at the extreme values of the parameters.For instance, the amplitude in ( 62) is singular when a = 0 reflecting these peculiarities.
At the upper critical point, Λ → 1 and (62) suggests a logarithmic decay c 1 ∼ (ln t) −1 .This logarithmic decay agrees with (50) that was carefully derived at the upper critical point λ = 2 in the model with a = 1.

D. Models with arbitrary exponents (a, s)
Addition-and-shattering processes with algebraic rates A k = k a and S k = λk s appear analytically intractable when the exponents (a, s) are arbitrary.The only exception is the steady-state regime.Below we combine analytical results for the steady states with simulations and general expectations gained from analytically tractable addition-shattering models studied earlier.
The governing equations for the monomer density c 1 (τ ), cluster densities c k (τ ), and moments M a (τ ) and M 1+s (τ ) are given by (7a), (7b) and (8).Numerically, we observe only two regimes-the system reaches a jammed state or a steady state.These regimes are demarcated by the critical shattering rate λ c that depends on the exponents (a, s) and initial conditions.
As expected, the critical shattering rate is an increasing function of the exponent a and a decreasing function of the exponent s, see Fig. 5(a).The maximal modified time τ max is a decreasing function of the exponent s but has a more complicated dependence on the exponent a: For large s, the time τ max is an increasing function of a; for small s, it may be a non-monotonous function of a, see Fig. 5(b).Our simulations show that the critical shattering rate λ c depends on the initial conditions.The stationary densities c 1 and c k (∞) are found from (7a)-(7b): Here we shortly write Combining mass conservation, k≥1 kc k (∞) = 1, with (65), we fix the final density of monomers: The product (66) exhibits qualitatively different large k behaviors depending on whether σ < 1 or σ > 1: If s > a − 1, we should use (68a) and then (67) gives c 1 (∞) > 0. Thus when s > a − 1, the steady state is possible.When s < a − 1, we should use (68c) implying that the sum in (67) diverges [recall that a ≤ 1].Thus the system falls into a jammed state if s < a − 1.
These simple arguments explain why interesting behaviors in models with algebraic rates occur when the exponents are related via s = a − 1. Equation (68b) implies that when λ > 2 − a the sum in (67) converges; c 1 (∞) > 0 and the system is in the steady state when s = a − 1 and λ > 2 − a.When λ < 2 − a, the sum in (67) diverges implying that c 1 (∞) = 0 and the system is either in the jammed state or in the SCS.The analysis for the class of models with s = a − 1 presented before shows that the jammed state arises when λ < 1 while the SCS emerges in the range 1 < λ < 2 − a.
The final densities in the steady state for the model with s = a − 1 simplify when a is an integer.We already know the answers when (a, s) = (1, 0) and (0, −1).In the next example (a, s) = (−1, −2), the densities in the steady state (λ > 3) are Cluster densities undergo a discontinuous (first order) phase transition at λ = λ c so that the jammed densities c k (∞), for λ = λ c , differ from the equilibrium densities c k (∞) for λ = λ c + 0. This is illustrated in Fig. 6.

III. CONTINUOUS, DISCONTINUOUS AND WEAK PHASE TRANSITIONS
For addition-and-shattering processes with algebraic rates, A k = k a and S k = λk s , two regimes depending on the parameters (a, s, λ) generically arise.In the jamming regime, monomers disappear, and evolution stops.In the steady-state regime, the evolution continues forever, but the densities are stationary (in the infinite size limit).For instance, if s > a − 1, the jamming regime occurs when λ < λ c while the steady-state regime occurs in the complimentary λ > λ c range.The transition between these regimes is discontinuous.The magnitude of the critical shattering parameter λ c is non-universal (it depends on the initial condition).
Models with (a, s) = (a, a − 1) exhibit particularly rich behaviors.There are three different regimes: 1.The jammed regime in the range λ < λ c,l .The jammed regime is non-universal as it depends on the initial condition.The magnitude of the lower critical shattering parameter is universal: λ c,l = 1.
2. The supercluster states (SCSs) occur in the range 3. The steady state regime in the range λ > λ c,up .The steady state densities are universal (independent on the initial condition).
We now briefly discuss weak phase transitions in models with (a, s) = (a, a − 1) occurring on the boundary between jammed and supercluster states at the low critical point λ = λ c,l = 1.There are infinitely many critical values a p where transitions occur.The critical values are defined by It is convenient to set a 1 = 1, so a 2 is defined by The To appreciate the emergence of phase transitions at these values we start with the governing equations for the model with exponents (a, a − 1) and λ = 1: Equation (72a) yields c 1 = e −2τ , and then from (72b) with k = 2 one gets If 2 a + 2 a−1 > 2, equivalently a > a 2 , the asymptotic behavior of the density of dimers is c 2 ≃ B 2 e −2τ with . All densities decay similarly and only amplitudes vary: when k ≥ 1 and a > a 2 .To determine the amplitudes we specialize the Laplace transform (55) to λ = 1 and use the Laplace transform c 1 (p) = 1 p+2 of c 1 = e −2τ .We get This Laplace transform has a simple pole p = −2 with residue which is the amplitude in (74).When a = a 2 , the density of dimers is c 2 = τ e −2τ , and generally Similarly to (55) we derive where we have used the Laplace transform c 2 (p) = 1 (p+2) 2 of c 2 = τ e −2τ .Thus c k (p) has a pole of order 2 at p = −2 with amplitude which yields the amplitude appearing in (77).Similarly in the a 3 < a < a 2 range for k ≥ 2 with amplitudes Continuing these calculations we find the decay laws for the densities.In terms of the physical time with α p = (p ap + p ap−1 )/2.In contrast to the phase transition between jamming and steady-state regimes where the final densities are finite and undergo a jump across the transition point, the cluster densities vanish in the present case, and only the decay exponent jumps from α p−1 to α p when the transition point a p is crossed.Therefore we call these phase transitions weak.Finally, we consider the total cluster density at the lower critical point (λ = 1).When a = 1, the individual cluster densities decay algebraically according to (19b), and the total cluster density also decays algebraically, Eq. (19a).An algebraic decay (82) of cluster densities when a < 1 suggests a similar behavior of N (t), yet much slower logarithmic decay occurs for all a < 1.In the a ↑ 1 limit, the decay law (83) reduces to the exact solution for the cluster density, N = (1 + 2t) −1/2 , of the model with a = 1 and λ = 1, see Fig. 8.When a = 0, the asymptotic (83) reduces to (26).In the a < 1 range, the logarithmic decay law (83) is derived below, Eq. (114), using essentially the same arguments as in the derivation of (26).

IV. THE NATURE OF THE SUPERCLUSTER STATES
We have shown that addition-and-shattering processes with rates A k = k a and S k = λk a−1 exhibit intriguing behaviors in the range 1 ≤ λ ≤ 2−a.All cluster densities decay to zero, so it is neither a steady state where final densities remain positive nor a jammed state where only monomers disappear in the final state.The above results refer to an infinite system.To shed light on the nature of supercluster states we analyze a finite system initially composed of M ≫ 1 monomers.The evolution stops when the last monomer disappears.A naive criterion gives an estimate of the time t * when the last monomer disappears.Equation ( 63) tells us that N (t)/c 1 (t) remains finite when 1 < λ < 2 − a.This apparently implies that the total number of clusters N ∞ remains finite: Simulations disagree with (85) and indicate that N ∞ diverges with system size (see Fig. 9) : In the jamming and steady-state regimes N ∞ ∼ M, while in the SCS the growth is sub-linear: δ < 1.The final mass distribution in the CSC has a scaling form from which we express through δ the exponents characterizing the scaled final mass distribution: γ = 1 − δ and α = 2δ − 1. Therefore one anticipates more tie bounds on the exponent δ, viz. 1 2 ≤ δ ≤ 1, which were indeed observed in simulations (see Fig. 9).
The time t * when the last monomer disappears also scales algebraically with system size To deduce t * one should not use the naive criterion (84).The evolution of monomers just before the system reaches the SCS (Fig. 10) hints at the flaws in reasoning based on (84).Significant fluctuations in the number of monomers close to t * indicate that relying on the deterministic average number of monomers Mc 1 (t) is questionable.Secondly, just before the system drops into the SCS, the number of monomers is still very large instead of being of the order of one as posited by Eq. (84).
To account for fluctuations in finite systems it is customary to employ the van Kampen expansion [58].Van Kampen expansions have been used in several areas (see, e.g., [44,[58][59][60][61]) including aggregation and annihilation processes [62][63][64][65][66][67].The idea is to decompose quantities of interest into extensive deterministic components and sub-extensive stochastic components.In the present case The terms linear in M are deterministic: the densities c k (t) obey the rate equations describing infinite systems.
The terms proportional to √ M are are stochastic: η k (t) are (evolving) random variables.The magnitude √ M of stochastic terms agrees with the basic tenets of statistical physics and probability theory.
Fluctuations are negligible in the thermodynamic limit (see Fig. 4 comparing simulation results for a large system with the mean-field solution corresponding to the infinite system).In our system, the van Kampen expansion is consistent in the jamming regime.The formation of the SCS, however, is dominated by fluctuations.The van Kampen expansion is quantitatively incorrect at the late stage, but it predicts qualitative features such as the exponents δ and β.The steady-state regime is quasisteady when M is finite.The van Kampen expansion is applicable for times that scale exponentially with M, but it does not describe a rare fluctuation that eventually leads to the extinction of monomers.
Analysis of finite systems is, in principle, straightforward but technically involved.Instead of rate equations describing the infinite system, we must rely on stochastic rules describing addition or shattering events.The state of the system is quantified by C = (C 1 , C 2 , C 3 , ...), where C k (t) is the number of clusters of size k.The quantities C k (t) are non-negative integers satisfying In an elementary reaction, C may transform as follows: We have shown only the components of C that changed and presented reaction rates.The reaction channel shown at the top describes the formation of dimers.The factor M −1 is necessary to recover Eqs.(3a)-(3b) in the M → ∞ limit.The following reaction channels in (92) describe addition and shattering involving clusters with k ≥ 2. To numerically investigate the supercluster state in finite systems, we use the Monte Carlo (MC) technique, namely, the Gillespie algorithm [57].
Figure 10 shows the evolution of the monomer density obtained by MC together with the mean-field (MF) behavior.At the beginning of the evolution, the MC and MF dependencies almost coincide.At latter times, c 1 (t) = C 1 (t)/M experiences large fluctuations and abruptly drops to zero.The intensity of fluctuations depends on the exponent a and the shattering coefficient λ.Large fluctuations reflect that close to the emergence of the SCS, clusters tend to be large.The shattering of such clusters produces a large number of monomers.This "duel" between addition and shattering eventually leads to the disappearance of monomers.
Using stochastic rules (92) we deduce Equations (93c) are valid for k ≥ 3. The evolution equations (93) depend on the second order moments.One can write exact equations for these moments but they depend on the third order moments.This hierarchical patterns continues making the system analytically intractable.Some progress is possible, however, upon relying on the van Kampen expansion.Consider the simplest moment ⟨C 2  1 ⟩.Using (92) we deduce the exact evolution equation Combining (93a) and ( 94) we find the exact equation for the variance Using the van Kampen expansion and the shorthand notations we re-write the terms appearing in (95) as where we have taken into account that ⟨η k ⟩ = 0. Plugging these expansions into (95) and equating the leading terms of the order O(M 2 ) we arrive at To close this equation one needs equations for covariances W 1k with all k ≥ 2. The only exception is the case of λ = 1 when the term with covariances vanishes.
For the model with (a, s) = (1, 0) [Sec.II A], the SCS occurs at λ = 1, so this is the only interesting value of the shattering parameter.Equation (96) becomes Using c 1 = e −2τ and M 2 = 2e τ − 1 following from the exact solution (17b) we reduce (97) to The variance diverges with time while the density vanishes: V 1 ∼ t 1/2 and c 1 ∼ t −1 .Thus fluctuations eventually dominate.
For the model with (a, s, λ) = (1, 0, 1) the total number of monomers is The deterministic part decreases with time while the stochastic part increases as rather than the naive criterion (84).Plugging V 1 ∼ t 1/2 and c 1 ∼ t −1 into (100) leads to (89) with β = 2/5.Recalling that N (t) = 1/ √ 1 + 2t, see (19a), we get i.e., δ = 4 5 .The exponents α and γ describing the scaled final mass distribution are therefore α = 3 5 and γ = 1 5 .Simulations suggest that when the system approaches the SCS, fluctuations rapidly drive it towards the final jammed state without monomers (see Fig. 10).Thus we estimate the final mass distribution in the SCS and final jammed state as we confirm the values of exponents and even get Φ(κ) = e −bκ with some unknown amplitude b.The prediction for the scaled average mass distribution is uncontrolled as it relies on the mean-field solution in the region where fluctuations become important, but it is in a fair agreement with simulations.
As another example, consider the model with exponents (a, s) = (0, −1).The SCS occurs [Sec.II B] when 1 ≤ λ ≤ 2. When λ = 1, Eq. (96) becomes Using c 1 = e −2τ and N ≃ τ −1 , see (26), we find Plugging this into (100) yields (89) with β = 1/2.Another exponent has the maximal value δ = 1, albeit the scaling law (86) acquires a logarithmic correction Generally if λ = 1 and a < 1, Eq. ( 96) simplifies to The moment M a+1 dominates the right-hand side when τ ≫ 1.To establish its asymptotic we return to the governing equation (72b) which we re-write in the form obtained by treating k as a continuous variable as we have done in deriving Eq. ( 24) which follows from Eq. ( 107) when a = 0. Introducing an auxiliary variable we re-write Eq. (107) as from which To determine the asymptotic behavior of M a+1 we again rely on a key feature of c k (τ ), namely that it has a sharp maximum at and hence Eq. (106) gives Plugging this into (100) leads to the scaling law (89) with β = 1/2 and a logarithmic correction: The same steps as in computing (111) give [cf.Eq. ( 83)] Hence V. PARTIAL DISINTEGRATION Here we demonstrate the emergence of the SCS in systems with partial fragmentation resulting in abundant production of monomers.One such model posits that a significant part of an aggregate (say a half) is shattered into monomers while the other part is kept whole.We analyze a more symmetric model of abundant incomplete disintegration.As previously, we postulate that a clustermonomer leads to absorption of monomer with the rate A k .Another possibility is fragmentation: A k−mer may break into k monomers, or a dimer and k − 2 monomers, or a trimer and k − 3 monomers, etc.All these breaking events are assumed equiprobable, so a cluster of size k may disintegrate into k − 1 equiprobable ways.This type of disintegration is schematically described by For the complete disintegration, we assume the homogeneous kernels, A k = k a and R k = λk r .Thus the probability that a cluster of size k breaks into a chunk of size j and k − j monomers is λk r /(k − 1) for all j.
The density of monomers obeys while for k ≥ 2 For instance, the first two terms in the right-hand side of (117a) describe the loss of monomers in addition process, the third terms represents the gain of monomers from incomplete disintegration and the last term gives the gain due to complete disintegration.In terms of the modified time, Eqs. (117) become Summing these equations we arrive at the evolution equation for the cluster density: As in the case of complete disintegration, there exists a critical value λ c such that for λ ≤ λ c , the system evolves to a jammed final state, while for λ > λ c it reaches an equilibrium state.The evolution of the monomer density is also the same: When λ < λ c , it decays to zero at τ = τ max , where the system falls into the jammed state.When λ > λ c , the monomer density is always positive, and the system reaches a steady state.When λ = λ c , the system arrives at the jammed state at τ = τ max , where c 1 (τ max ) = 0, and A. Supercluster states Similar to systems with shattering, the SCSs exist if the exponents shifts by one: r = a − 1.In this case τ max → ∞ for λ c,l ≤ λ ≤ λ c,up and c k (∞) = 0 for this range of λ.At the lower critical point, the transition from the jammed state with c 1 (∞) = 0 and c k (∞) ̸ = 0 for k ≥ 2 to the SCS with c k (∞) = 0 for all k takes place.The transition from the SCS to the steady state occurs at the upper critical point.The evolution of the monomer density for λ c,l ≤ λ ≤ λ c,up is similar to that for the case of complete disintegration (see Fig. 12).
To find the upper critical point, we recall that the phase transition is continuous so that c k (∞) vanish when λ → λ c,up .Consider clusters heavier than monomers, I = N − c 1 .Subtracting the first equation in (119) from the last we find that I satisfies In the steady state which we use to recast the second equation in (119) into for k ≥ 2. From this recurrence we deduce which together with mass conservation Hence we obtain the critical interval associated with the SCS: Simulations agree with these predictions (see Fig. 13).
Closed-form expressions can be derived for integer a.If (a, r) = (0, −1), the densities c k (∞) and N ∞ are In the subcritical regime, the monomer density vanishes; in the supercritical regime, it quickly reaches a positive value.These results agree with the MF predictions.In the critical regime, fluctuations dominate for large time.
In this model, the SCS occurs at a single point λ c = 2.The model with parameters (a, r, λ) = (1, 0, 2) appears most tractable for analytical exploration of the SCS in systems with abundant incomplete disintegration.The governing equations cannot be solved analytically, but appear amenable to asymptotic analysis leading to and suggesting the scaling form c k (t) = k −1 t −1 Φ(k/t).

B. Finite Systems
In a single reaction event in a finite system, the configuration C = {C 1 , C 2 , . ..} may transform into one of the following configurations: We show only the components that have been altered.The top channel represents merging of two monomers, the next describes addition of monomer to clusters of mass k ≥ 2, and the last channels describe disintegration of clusters of mass k ≥ 2. Figure 14 shows the behavior of the monomer density C 1 (t)/M in systems with M = 10 5 − 10 7 monomers.The dominance of fluctuations for λ inside the critical interval associated with the SCS is visible.Outside the critical interval, the cluster densities obtained from the MF rate equations agree with MC data.In the critical interval, fluctuations dominate over the MF predictions.
We have shown that systems with abundant incomplete disintegration demonstrate qualitatively similar behaviors to systems undergoing shattering.In particular, the SCSs are also non-extensive, with the final number of clusters exhibiting a sub-linear scaling in M, see Fig. 15.The transition time to the SCS also has a power-law dependence on the system size (Fig. 15).

VI. CONCLUSIONS
We have investigated the evolution of aggregates triggered by collisions with monomers.A collision may result in the attachment of a monomer or the break-up of an aggregate into constituting monomers.We have assumed that addition and shattering rates vary algebraically with size, A k = k a and S k = λk s .Without loss of generality, we set the amplitude in the addition rate to unity; the amplitude in the shattering rate is denoted by λ, so the reaction rates are parameterized by (a, s, λ).Depending on these parameters, we observed three types of behavior: 1.The system falls into a jammed state where monomers disappear and the evolution stops; other cluster densities depend on the initial condition.
2. The system reaches a steady state with cluster densities independent on the initial condition.
3. The system reaches a supercluster state, a peculiar jammed state where all densities vanish in the thermodynamic limit.
The SCSs occur when s = a − 1 and the shattering amplitude lies in the interval 1 ≤ λ ≤ 2 − a.The phase transitions at λ c,l = 1 and λ c,up = 2 − a are continuous.When λ = 1, there is also an infinite set of weak phase transitions at a n (n = 1, 2, . ..) manifested by an abrupt change of the exponents characterizing the power-law decay of the cluster densities.
The nature of the SCSs is best revealed in finite systems initially composed of M ≫ 1 monomers.The SCSs are peculiar jammed states.The time to reach the SCS scales as M β , sometimes with a logarithmic correction [cf.Eq. ( 113)].The final number of clusters in the SCS also scales algebraically, N ∞ ∼ M δ , sometimes with a logarithmic correction.The final number of clusters is sub-linear in system size, 0 < δ < 1.We have argued for stronger bounds 1  2 ≤ δ ≤ 1.The extreme value, δ = 1, supplemented by a logarithmic correction, is feasible, see Eq. ( 115), so N ∞ continues to scale sub-linearly with M.
The SCSs are born in a fluctuation-dominated process.This results in non-self-averaging and non-extensive characteristics of the SCSs.The van Kampen expansion becomes dubious close to the birth of the SCS, albeit we showed how to use it to probe the basic features of the SCS, e.g., the exponents β and δ.
We have argued that systems with abundant incomplete disintegration exhibit qualitatively similar behaviors to systems undergoing shattering.The SCSs again emerge with exponents shifted by one: r = a − 1.We have little theoretical understanding of the SCSs in systems with abundant incomplete disintegration since we do not have analytical solutions for the time-dependent densities c k (t).The model with (a, r) = (1, 0) in which the SCS occurs at a single point λ = 2 looks most interesting and tractable, and establishing the asymptotic behavior of c k (t) looks achievable [cf.(127)].However, even in that model, we do not see how to obtain a closed equation for the variance ⟨C 2 1 ⟩ − ⟨C 1 ⟩ 2 , like Eq. (97) for complete disintegration.We need such an equation for finding the exponents β and δ.
Another possible generalization is to postulate that monomers and a few other light cluster species are mobile and consider the processes triggered by collisions with them.The pure addition processes of this kind are already analytically intractable when there are two mobile species, e.g., monomers and dimers.However, the phenomenology is the same, viz., the system quickly reaches a jammed state.Exploration of SCSs in additionshattering processes of this type is a challenge for future work.
We emphasize that in a finite system, a steady state is quasi-steady as monomers eventually disappear in a rare fluctuation, and the system gets jammed.The (average) lifetime T of quasi-steady states is astronomically large, namely, it scales exponentially with system size ln T ≃ CM. (128) Giant adsorption times like (128) arise e.g. in population dynamics where they are known as extinction times.
Wentzel-Kramers-Brillouin (WKB) technique is a powerful toolbox for finding the controlling exponential behavior (128).Single-population models admit analytical treatment, see, e.g., [68][69][70][71][72][73] and a review [74].An intriguing challenge is to develop the WKB for additionshattering processes possessing many populations (cluster species).Systems of two interacting populations are generally intractable analytically [74].Therefore computing the amplitude C in (128) could be impossible but probing it numerically is perhaps feasible in the WKB framework.Supercluster states may arise in collision-controlled finite systems where jamming is inevitable, and the overall rates of competing merging and disintegration mechanisms are comparable.In addition-shattering processes, when the merging and disintegration rates are comparable, A k ∼ kS k , the exponents must obey s = a − 1 when rates are algebraic.Indeed, we found the supercluster states in this situation with an extra constraint on the shattering intensity: 1 ≤ λ ≤ 2 − a.In addition-chipping processes, the rates are comparable when A k ∼ C k .Therefore when the exponents are equal, s = a, the supercluster states could emerge in addition-chipping processes, and they do emerge [75] when the chipping intensity obeys λ = 1.
So far the supercluster states have been detected only in the mean-field framework.Finding them in finite dimensions is an intriguing challenge.

FIG. 5 .
FIG. 5. (a) The critical shattering rate λc as a function of a and s.(b) The maximal modified time in the critical regime, τmax(λc), as a function of a and s.The kinetic rates A k = k a and S k = λk s ; the mono-disperse initial condition is used.

FIG. 6 .
FIG. 6.The final cluster densities undergo a discontinuous phase transition at λ = λc(a, s) as illustrated for a = 0.6 and s = 0.2 (top panel) and for a = 0.2, s = 0.6 (bottom panel).The mono-disperse initial condition is used.

FIG. 7 .
FIG. 7. Illustration of the phase diagram of the infinite system with aggregation and shattering.Left panel: For general addition and shattering exponents (a, s), the system undergoes a discontinuous phase transition from a jammed to a steady state at λ = λc.On the s = a − 1 line, the system undergoes a continuous phase transition.The right panel details the behaviors on the line s = a − 1 of the left panel.The continuous phase transition from a jammed to a supercluster state occurs at λ = 1.The continuous phase transition from a supercluster state to the steady state occurs at λ = 2 − a.On the λ c,l = 1 line, there is an infinite set of weak phase transitions at a = a1, a2, . . .ap, . ...
) describing the average ⟨C k ⟩ of the total number of clusters of mass k.The total numbers C k (∞) = C k (t * ) are non-self-averaging random quantities (they significantly vary from realization to realization) in the SCS.Combining the scaling form (87) with (86) and mass conservation leads to relations α + γ = δ, α + 2γ = 1 (88)

FIG. 13 .
FIG. 13.The final densities of monomers (left panel) and clusters (right panel) for the model described by Eqs.(119) with exponents (a, r) = (0, −1).Curves from top to bottom correspond to the densities, c2, c3 and c4.At λ c,l = 2, the system undergoes a continuous phase transition from a jammed state to the SCS; at λc,up = 3, it undergoes a continuous phase transition from the SCS to the steady state.The final densities in the jammed state depend on the initial conditions.Solid lines: c1(0) = 1; dashed lines: c1(0) = 0.2, c2(0) = 0.4.Curves are solutions of rate equations; Monte Carlo results are shown by dots.

FIG. 15 .
FIG. 15.Left panel: Final cluster density, N∞, as a function of system size M. Dots show MC results; line is a fit to N∞ ∼ M δ with δ = 0.6777.Right panel: The transition time t * as a function of system size.Dots show MC results; line is a fit to t * ∼ M β with β = 0.7835.The parameters are (a, r, λ) = (0, −1, 2).