The Smoluchowski Ensemble—Statistical Mechanics of Aggregation

We present a rigorous thermodynamic treatment of irreversible binary aggregation. We construct the Smoluchowski ensemble as the set of discrete finite distributions that are reached in fixed number of merging events and define a probability measure on this ensemble, such that the mean distribution in the mean-field approximation is governed by the Smoluchowski equation. In the scaling limit this ensemble gives rise to a set of relationships identical to those of familiar statistical thermodynamics. The central element of the thermodynamic treatment is the selection functional, a functional of feasible distributions that connects the probability of distribution to the details of the aggregation model. We obtain scaling expressions for general kernels and closed-form results for the special case of the constant, sum and product kernel. We study the stability of the most probable distribution, provide criteria for the sol-gel transition and obtain the distribution in the post-gel region by simple thermodynamic arguments.


Introduction
Aggregation is the process of forming larger clusters through the merging of smaller clusters.This generic process is encountered in a large variety of systems, from polymerization and colloidal aggregation to the clustering of social groups and the merging of galaxies.The mathematical foundations of aggregation were set by Smoluchowski (von Smoluchowski, 1906), whose particular problem interest was in Brownian coagulation.The aggregation equation, more commonly known as Smoluchowski equation, is rate equation on the distribution of particles whose size (mass) changes via binary aggregation events.For a discrete population of clusters with integer masses in multiples of a unit mass ("monomer" or "primary particle") it takes the form where c i is the number concentration of clusters with mass k and K i,j the aggregation kernel, a rate constant for the merging of masses i and j.A large body of literature has focused on the theory of the Smoluchowski equation, the existence of analytic solution and the scaling limit (Leyvraz, 2003).Of particular interest is gelling, a condition that arises under the product kernel K i,j = ij and refers to the formation of a giant structure as in polymer gels and which is manifested by the failure of the Smoluchowski equation to conserve mass.This process is commonly described as a phase transition and suggests the possibility that statistical thermodynamics, a theory developed for equilibrium states of interacting particles, may perhaps be applicable in this problem.Studies of Smoluchowski aggregation fall in one of two categories, kinetic and stochastic.The kinetic approach is based on Eq. ( 1) and its solution.Stable solutions of the Smoluchowski equation conserve mass.Gelling is identified as the point where mass conservation breaks down (Ziff et al., 1982;Hendriks et al., 1983).Post-gel solutions require additional assumptions as to how the gel and dispersed phases interact (Ziff and Stell, 1980).The stochastic approach views clusters as entities that merge with probabilities proportional to the aggregation kernel.It was first formulated by Marcus (Marcus, 1968) for a discrete finite population and its formal mathematical treatment was developed by Lushnikov, who obtained solutions for certain special cases including gelation (Lushnikov, 1978(Lushnikov, , 2011(Lushnikov, , 2005b(Lushnikov, ,a, 2004)).In Lushnikov's method each possible distribution is given a probability whose evolution in time is tracked via a generating functional.The approach is explicitly probabilistic and views the Smoluchowski equation as the deterministic approximation of the underlying stochastic process (Aldous, 1999).A different approach within the probabilistic realm makes use of combinatorial methods.This treatment originated with Stockmayer (1943) and was further explored by Spouge (Spouge, 1985(Spouge, , 1983b;;Hendriks et al., 1985;Spouge, 1983a).The combinatorial approach considers the number of ways to form a particular distributions of clusters and assigns probabilities in proportion to the multiplicity of distribution.The ensemble of distributions is then reduced to the most probable distribution, which is identified by maximizing the combinatorial multiplicity.This approach has two appealing advantages.It treats a time-free ensemble in which time appears implicitly via the mean cluster mass, which acts as the progress variable.Secondly, it brings the problem closer to the viewpoint of statistical mechanics and the notion that an ensemble may be represented in the scaling limit by its most probable element.Stockmayer saw this connection and his treatment of gelation is replete with references to the theory of phase transitions (Stockmayer, 1943).The analogy between aggregation and thermodynamics was not formalized, however.Stockmayer obtained the gel point by mathematical, not thermodynamic tools, and arrived at a post-gel solution that is not consistent with the kinetics of aggregation (Ziff and Stell, 1980).
We have previously shown that gelation can be indeed treated as a formal phase transition and presented solutions for the product kernel in the pre-and post-gel regions using a thermodynamic approach (Matsoukas, 2015b) building on our earlier work on the cluster ensemble (Matsoukas, 2015a(Matsoukas, , 2014)).Here we generalize the methodology to formulate a rigorous thermodynamic theory of Smoluchowski aggregation.The treatment begins with a finite population that starts from a well defined state and construct the set of all possible distributions that can be reached in a fixed number of elementary transitions.The probability of distribution in this ensemble is governed by the kinetics of the elementary processes that act on the population.In the thermodynamic limit the most probable distribution is overwhelmingly more probable than all others and is governed by a set of mathematical relationships that we recognize as thermodynamics.The work is organized as follows.In Section 2 we define the Smoluchowski ensemble of distributions and their probabilities.In Section 3 we formulate the probability in terms of a special functional W that introduces the partition function and the Shannon entropy of distribution.In Section 4 we treat the scaling limit and derive the thermodynamic relationships of the Smoluchowski ensemble.In Section 5 we obtain solutions of the Gibbs form for the classical kernels, constant, sum and product.We analyze the stability and phase behavior of the ensemble in Section 6 and treat the sol-gel process as a phase transition.In Section 7 we express the results in the continuous domain and finally offer concluding remarks in Section 8.

The Smoluchowski Ensemble
We consider a population of clusters composed of i = 1, 2 • • • units (monomers).In binary aggregation two clusters merge to form a new cluster in which the total number of members is conserved via the schematic reaction (2) The merging of a pair constitutes an elementary stochastic event whose probability depends on the size of the reactant clusters via the aggregation kernel K i,j .At the initial state the population consists of N 0 = M single members (monomers).This distribution constitutes generation g = 0.The next generation is constructed by implementing every possible aggregation event in the distribution of generation g = 0.
The set of distributions formed constitutes the microcanonical ensemble of generation g = 1.We continue recursively to form the ensemble of distributions in generation g by implementing all possible aggregation events, one at a time, in all distributions of the parent ensemble.We represent a distribution of clusters by the vector n = (n 1 , n 2 • • • ), where n i is the number of clusters with i members.All distributions in generation g satisfy the conditions The first condition expresses the fact each elementary event decreases the number of clusters by 1 according to the stoichiometry of the transition and the second expresses the fact that the number of members is conserved.Conversely, all distributions that satisfy the conditions in Eq. ( 3) are members of the ensemble of generation g because they can be formed in g steps starting from a distribution of M monomers.We view the two equations in Eq. ( 3) as the constraints that define the ensemble of feasible distributions.We call this ensemble microcanonical to indicate that it is conditioned by two extensive constraints that fix the mean cluster mass M/N = x in all distributions of the ensemble.The progression of the ensemble in the direction of advancing generation may be represented in the form of layered graph (Fig. 1).Vertices represent distributions and edges represent connections via elementary transitions.Edges are directed from parent in generation g − 1 to offspring in generation g.Layers are organized by generation and contain all distributions in the generation.The graph is begins with the initial distribution of monomers and is connected until it reaches the terminal state in which all particles have joined the same cluster.Aggregation may be viewed as a random walk on this graph.A trajectory is a possible sequence of connected edges from the initial state (N = M) to the final state in which all mass is contained in a singe cluster (N = 1).Our goal is to establish the probability of any distribution in any generation in terms the aggregation kernel K i,j .This probability also depends on the initial state.For us, this will be fixed to a state of all monomers.

Kinetics
When cluster masses i and j − i in distribution n ′ of generation g − 1 merge to form a cluster of mass i parent distribution n ′ is transformed to offspring distribution n ′ via the transition This transition is represented by an edge in the graph of Fig. 1.The rate of the transition is proportional to the number of ways to choose the reactants and the aggregation kernel, The total rate by which parent n ′ produces offspring is where These kinetic equations summarize the Smoluchowski aggregation model.

Probabilities
We assign a probability to each distribution n within each generation g and formulate the propagation of probability between generations as follows: Here n is a distribution in generation g, n ′ is a parent of n in generation g − 1 via the transition shown in Eq. ( 4) and T(n ′ → n) is the transition rate from parent to offspring, defined as where R g−1 is the mean propagation rate from parent generation g − 1 to the next: with the summation over all distributions n ′ in generation g − 1. Expressing the transition rate in terms of the aggregation kernel we obtain with We may confirm that P(n) as defined in Eq. ( 8) satisfies normalization over all distributions in the same generation.Beginning with P(n 0 ) = 1 at the initial state, Eq. ( 8) uniquely determines the probabilities of all distributions in all future generations.

Smoluchowski Equation
The mean number of cluster with mass k in generation g is with the summation going over all distributions in the same generation.We will derive the evolution of the mean distribution from parent generation g − 1 to generation g.The probability of distribution P(n) is given by Eq. ( 8) and is expressed as a summation over its parents n ′ .By the stoichiometry of the transition the parent and offspring distributions satisfy Combining this relationship with ( 13) and (8) the result is (see Supplementary Material) The left-hand side is the change in the mean number of k-mers between generations; the right-hand side is the ensemble average of the generation and depletion of k-mers within the distributions of the ensemble.Define the mean time increment ∆t from parent to offspring generation as then Eq. ( 15) reads In the mean field approximation we set n = n to obtain This is the Smoluchowski equation for discrete finite binary aggregation.Here we have obtained it as the mean field approximation of the stochastic process defined in Section 2.2.The exact result is Eq. ( 17).The mean field approximation implies that a single distribution is representative of the entire ensemble.We examine the conditions under which this is true in Sections 4 and 6.

Partition Function and Selection Functional
We express the probability P(n) in generation g in the form where n! is the multinomial coefficient of vector n, g is the number of clusters in all distributions of generation g, W(n) is a functional of distribution n, to be determined, and Ω G,N is the partition function and ensures proper normalization of P(n), with the summation over all distributions in generation g = M − N. By expressing P(n) in the form of Eq. ( 19) the calculation of P(n) becomes a problem of determining W(n) Ω M,N .

Shannon Entropy
The multinomial coefficient represents the "natural" multiplicity of distribution n, namely, the number of ways to order the clusters in the distribution if clusters with the same number of particles are treated as indistinguishable.The log of the multinomial coefficient in the Stirling approximation log It is a concave functional of n with functional derivatives It is homogeneous in n with degree 1 and satisfies the Euler condition We set q i = n i /N and apply H to vector q.By homogeneity the result is written as In this form H reverts to the familiar entropy functional, historically associated with Boltzmann, Gibbs and Shannon.We will call it Shannon functional and avoid the term "entropy," whose precise meaning varies across disciplines.For our purposes the Shannon functional is defined as and may be applied to any vector a with non-negative elements regardless of normalization.

The Selection Functional
Functional W(n) biases the statistical weight of distribution n relative to its natural multiplicity.We call it selection functional because it effectively selects distributions relative to each other.The functional derivatives of log W is and define the cluster factors w i;n , a property cluster mass i in distribution n.In general log W the cluster bias w i;n depends not only on i as well on the distribution n on which these factors are evaluated.In the special case that log W linear functional of n the cluster factors are functions i alone and they are the same in all distributions.This special condition is associated with the Gibbs distributions discussed is Section 5.
If W(n) = 1 for all distributions, then the probability of distribution is proportional to its multiplicity n!.If this special condition is met we will call the ensemble unbiased.The partition function of the unbiased ensemble can be easily determined by a combinatorial argument: it is equal to number of ways to assign M objects in N groups and is given by Accordingly, the probability of distribution in this special case is In a population undergoing transformations the selection functional is determined by the kinetic details of the mechanisms that produce the transformations; in the case of aggregation it is determined by the aggregation kernel K i,j .Whether the unbiased ensemble is a possible solution under some kernel remains to be seen.

Propagation Equations
At the initial state all members are at the monomeric state, n i,0 = Mδ i,1 .We set W(n 0 ) = 1 and since n 0 != 1 we also have Ω M,M = 1.We insert Eq. ( 19) into Eq.( 8) and express the summation over parents of n as a summation over all all pairs (i − j, j) that produce mass i in distribution n.The result is (see Supplementary Material) Here N is the number of clusters in distribution n of generation g = M − N, n ′ is the parent distribution via the transition (i − j) + (j) → (i) and K M,N+1 is the mean kernel in the parent generation g ′ = g − 1.
The left-hand side of Eq. ( 30) depends solely on M and N whereas the second term on the right-hand side contains functionals of distribution n.This term must be the same for all distributions n in the same generation in order to produce a result that is a function of M and N alone.From Eq. ( 19) it is clear that W and Ω may be defined within a proportionality constant α M,N ; as long as this constant is common for all distributions it has no effect on the probabilities of the ensemble and may be chosen arbitrarily.We choose it to satisfy the following criterion: if W = 1 in some generation for all distributions of the ensemble, we require it to be 1 in the next generation as well.The choice that satisfies this condition is to set the double summation in Eq. ( 30) to 1. Equation ( 30) now splits into two separate recursions, one for the partition function, and one for the selection functional, The recursion for the partition function is readily solved to produce the partition function in generation g = M − N: Accordingly, the partition function is equal to the unbiased partition function times the product of the mean kernels from generation 0 up to the parent generation of the current ensemble.The recursion for the selection functional may be expressed as and in this form it gives the selection functional of the offspring as a linear combination of that of its parents.In principle this can be solved recursively for any distribution in any generation.For certain special cases the recursion can be solved in closed form.These are discussed in Section 5.

Most Probable Distribution
We define the scaling limit by the condition M, N → ∞ at fixed M/N = x.The expectation is that in this limit the intensive mean distribution n k /N must converge to a limiting distribution pk that is independent of M and N and depends only x only: We further anticipate that the probability of distribution P(n) becomes infinitely sharp around a single distribution, n * = Np * , such that p * k is not merely the most probable distribution, it is overwhelmingly more probable than any other distribution in the ensemble.This further that mean and most probable distributions converge to each other: This convergence is an implicit requirement for the validity of the Smoluchowski equation: the mean-field approximation is exact if a single distribution is representative of the entire ensemble.This is possible only if P(n) peaks very sharply about the most probable distribution.When a single term dominates the summation that defines the partition function in Eq. ( 21), the log of the sum converges to the log of the maximum term, As a further consequence of the intensive convergence in (35) we have the Euler relationship for log W: where log w * i is the functional derivative of log W(n * ), Equation ( 38) expresses the fact that log W must be homogeneous functional of the MPD.This condition follows from Eq. ( 37) and the homogeneity properties of H(n * ) and log Ω M,N .The most probable distribution (MPD) maximizes the probability in Eq. ( 19) among all distributions that satisfy the constraints in Eq. (3).By Lagrange maximization we obtain the MPD in the form and q and β are parameters related to the Lagrange multipliers.We insert the MPD into Eq.( 37) to obtain log Ω M,N = βM + (log q)N. (41) This fundamental equation relates the partition function to the primary variables of the ensemble: the macroscopic variables (M, N) that define the ensemble and the Lagrange multipliers that appear in the MPD.The convergence of n * k /N to intensive limit p * k implies that β and q are intensive, i.e., the are functions of x = M/N but not of M or N individually.This further implies that Eq. ( 41) is homogeneous function of M and N with degree 1 and thus must satisfy Euler's theorem, Direct comparison with Eq. ( 41) leads to: Thus the Lagrange multipliers that appear in the MPD are the partial derivatives of the partition function.Differentiation of Eq. ( 41) with respect to all variables that appear on the right-hand side gives This is the Gibbs-Duhem equation associated with the Euler equation for log Ω M,N in Eq. ( 41).It may be written as In this form its expresses the relationship between β, q and x.
The MPD maximizes the log of the microcanonical weight, H(n) + log W(n) and its maximum is log Ω M,N .Therefore we have the inequality: It is satisfied by all distributions n in the (M, N) ensemble with the equal sign only for n = n * .This is the fundamental variational principle of the ensemble: it defines the MPD and generates all relationships of this section.

Thermodynamics
The set of equations obtained above form a network of relationships all of which appear in statistical thermodynamics.Equation ( 40) is the generalized canonical distribution, a member of the exponential family, whose parameters β and q are related to the microcanonical partition function via Eqs.( 41), ( 43) and ( 44).We define the extensive form of the canonical partition function Q(β, N) via the Legendre transformation of log Ω: and thus we recognize q = Q 1/N as the intensive form of the canonical partition function.
The variational condition that produces the set of thermodynamic relationships is the inequality in Eq. ( 47), which defines the MPD as the distribution that maximizes the microcanonical weight.Expressing H(n * ) and log W(n * ) in terms of the Euler relationships ( 22) and ( 38), respectively, this inequality takes the form Eq. ( 40) Partition Function Ω M,N = βM + (log q)N Eq. ( 41) Eq. ( 43) Eq. ( 44) where p i = n i /N.The inequality is satisfied by all distributions p i with mean x = M/N and the equality applies only to p i = p * i .With w * i = 1 it reduces the second law: the log of the microcanonical partition function is equal to the entropy of the equilibrium (most probable) distribution and this is larger than the entropy of any other distribution with the same mean.We will see in Section (5) that the case w * i = 1 corresponds to the constant kernel.
Table 1 summarizes the results of this section.These relationships follow from the maximization of the probability is expressed in the form in Eq. ( 19) and are independent of the details of aggregation.These details enter through Eqs. ( 33) and (34) that define the partition function and the selection functional of aggregation as a function of the aggregation kernel.

Gibbs Distributions
In the thermodynamic formulation the probability of distribution is determined by the selection functional and this in turn depends on the aggregation kernel via the propagation equation Eq. (34).A special type of functional is of the form whose log is linear in with functional derivatives log w i .Here the w i are functions of i alone and do do not depend on n.If the selection functional is is given by Eq. ( 50) the probability of distribution in Eq. ( 19) takes the form Probability distributions of this type are called Gibbs distributions (Berestycki and Pitman, 2007) and are frequently encountered in stochastic processes (Kelly, 2011).Several important results can be obtained in analytic form.In particular, the mean distribution is (Matsoukas, 2019b): The result is exact for all We apply this selection functional of Eq. ( 50) to the transition (i − j) + (j) → (i) that converts parent distribution n ′ to offspring n.By the stoichiometry of the transition we have and with this Eq.( 32) gives One possible solution that satisfies this equation for all distributions n is This is not the only possible solution for W in Eq. ( 32), however, it may or may not be acceptable; if it is, we have obtained a Gibbs distribution and the kernel is a Gibbs kernel.
We have identified three kernels for which Eq. ( 56) is the correct solution.These are the constant kernel, the sum kernel and their linear combinations.The product kernel, a kernel of theoretical importance because its relationship to gelation, is a quasi-Gibbs kernel and is discussed in Section 5.3.We will provide detailed solutions for the constant and sum kernels.We will not discuss the linear combination in part because the results are more involved but mainly because this kernel reverts to the sum kernel when cluster masses are large thus it does not contribute to our understanding of aggregation beyond what we learn by studying the constant and sum kernels separately.

Constant Kernel
With K i,j = 1 Eq.( 56) gives w i = 1 for all i.The selection functional is W(n) = 1 for all n and we conclude that the ensemble is unbiased whose partition function is given by Eq. ( 28): The mean distribution follows from Eq. ( 53) and is given by To obtain the most probable distribution we calculate the parameters β and q from Eqs. ( 43) and ( 44) along with (59).The differentiations may be done by first replacing the factorials in the partition function with 5.2.Sum Kernel the Stirling expression.Alternatively we may obtain these parameters by the discrete difference form of these derivatives and apply the asymptotic conditions M, N ≫ 1.The latter method is simpler: The MPD is assembled according to Eq. ( 40) with w * k = w k : For large x this goes over to the exponential distribution which is the well known result for the constant kernel.Here x stands for the continuous cluster mass.

Sum Kernel
The ensemble average of the sum kernel is We obtain the partition function from Eq. ( 33).The result is The factors w i that satisfy Eq. ( 56) are and the mean distribution is obtained from (53), This is an exact result for all The parameters β and q are obtained similarly to those for the constant kernel: Combining with Eq. ( 40) we obtain the MPD in the form with θ = 1 − 1/ x.We use the the Stirling formula for the factorial the MPD in the continuous limit takes the form 5. Gibbs Distributions Figure 2: Approach to scaling limit for the sum kernel at x = 10.The MPD is calculated from Eq. ( 71) and the dashed lines from Eq. ( 68) at M/N = 10.
Figure (2) shows the MPD for x = 10 and the mean distribution from Eq. ( 68) at fixed M/N = 10 for various values of M and N.In the scaling limit the mean distribution converges to the MPD.

Quasi-Gibbs Kernels -The Product Kernel
We are able to obtain closed-form expressions for the partition function of the constant and sum kernels and heir linear combinations because they all satisfy the condition for all n.Simply put this states that the mean kernel is the same in all distributions of the ensemble, therefore also equal to the ensemble average kernel.In this case the calculation of the ensemble average kernel is trivial and does not require knowledge of the probabilities n.The constant kernel, sum kernel and their linear combinations are the only kernels that satisfy (73) in the strictest sense, i.e., for all n that satisfy the two constraints in (3).We will refer to Eq. ( 73) as the Gibbs condition because it generates Gibbs distributions.We may relax the requirement that all distributions obey the Gibbs condition with the milder requirement that it be obeyed by most distributions.We consider the product kernel, whose mean within distribution n is Here i = M/N and i 2 are the normalized first moment and second moments of n, respectively.In the limit N → ∞, M/N = fixed this scales as 5.3.Quasi-Gibbs Kernels -The Product Kernel in most distributions except those that contain clusters of the order M.1 According to Eq. ( 76) is a quasi-Gibbs kernel: it satisfies the Gibbs condition in Eq. ( 73) asymptotically in most but not all feasible distributions.We proceed to obtain the Gibbs distribution of the product kernel and test its validity.Inserting ( 76) into (33) we obtain the partition function: We complete the solution by evaluating w i from Eq. ( 56), The mean distribution is obtained by inserting these results into Eq.( 53): Unlike Eq. ( 63) and (68) this result is not exact.This can be demonstrated numerically by the fact this distribution is not normalized to unity and its mean is not M/N but approaches proper normalization in the asymptotic limit.This failure arises from the fact that Eq. ( 53) requires a Gibbs probability distribution that strictly applies to all distributions of the (M, N) ensemble.72) and the mean distribution (dashed lines) from Eq. ( 79).The distributions for x = 4 are not stable and do not represent the correct state of the population.

Gibbs Distributions
We obtain β and log q from Eqs. ( 43) and ( 44): and in the continuous limit These results are summarized in Table 2 along with those for the constant and sum kernels.
The relationship between the mean and the most probable distribution of the product kernel is shown in Fig. 3 for two values of the mean cluster, x = 1.75 and x = 4.At x = 1.75 the mean distribution calculated from Eq. ( 79) is not exact but its moments asymptotically approach the correct values as M and N are increased at fixed M/N = x.At x = 4 the behavior is different.A peak develops at the long tail of the distribution.It is pushed to ever larger sizes but never vanishes.This suggests the emergence of a gel phase but the mean distribution from Eq. ( 79) is not correct: its mean does not converge to x as M and N are increased but to a smaller value.As we discuss in the next section, the MPD at x = 4 is unstable and does not represent the state of the population.In the shaded region the system is stable and is represented by its MPD.The unshaded region is unstable and the system is split into two phases, a sol phase and a gel phase, each represented by its own MPD.Both graphs provide equivalent criteria of stability.
6 Phase Behavior

Stability
The fundamental inequality of the ensemble is Eq. ( 47) that defines the most probable distribution.This condition implies that the microcanonical functional is concave and in turn this implies that log Ω is a concave function of M and N and requires (see Supplementary Material) These equivalent conditions guarantee the existence of the MPD in the form of Eq. ( 40).In thermodynamic language they ensure that the MPD represents a stable state.The parameters β and q of the constant, sum and product kernel are plotted in Figs.4a and 4b, respectively, as a function of the progress variable θ = 1 − 1/ x.According to Eq. ( 84) stability requires β to be decreasing function of x and q increasing function of x.The constant kernel is stable at all θ: β decreases and q increases monotonically over the entire range of θ.The sum kernel is also stable at all θ but reaches the limit of stability at θ = 1 or x = ∞.This kernel is borderline stable: it is stable for all finite times and reaches instability at t = ∞.The product kernel is stable up to x = 0.5 beyond which point both β and q violate the stability criteria.
To survey the stability landscape of aggregation we employ the power-law kernel with arbitrary exponent ν ≥ 0. This is a homogeneous kernel with degree ν that reverts to the constant and product kernel for ν = 0 and ν = 2, respectively.We treat this as a quasi-Gibbs kernel by analogy to the product kernel.If we take the ensemble average power-law kernel to scale as it is a straightforward matter to show that the parameters β and q are Notice that with ν = 1 recover the sum kernel.This turns the power-law kernel into a useful tool, a homogeneous kernel that reproduces the correct (β, q) values for the constant, sum and product kernels, and which may be used to interpolate (and cautiously extrapolate) to other homogeneous kernels by varying ν.The stability limit in power-law aggregation is reached at accordingly, the MPD is stable in 0 ≤ θ * and unstable in θ * < θ ≤ 1.In Figs.
(4)a and (4)b the stable region is indicated by shading.For ν ≤ 1 the system is stable at all θ from 0 to 1.For ν = 1 the limit of stability appears at θ = 1, which is reached in infinite time.In practice the system is stable at all finite times.For ν > 1 the stability border appears in finite time when the mean size reaches the critical value For ν = 2 (product kernel) the limit of stability is reached at x * = 2.We see from Fig. 4 that both β and q reach the limit of stability at the same time.

Phase Splitting -The Sol-Gel Transition
When the system crosses into the unstable region its state is no longer represented by the MPD but by a mixture of two phases each with its own MPD.What are these phases?To answer this question we must recognize that the elements of the ensemble are fundamentally discrete distributions; the apparent continuity in the scaling limit is a mathematical artifact, a convenience but not a fundamental quality of the ensemble.We must consider therefore a finite system.Given a distribution of M particles arranged into N clusters, the maximum cluster mass that is possible is k max = M − N + 1 and is found in a single distribution of the ensemble, one in which one cluster contains M − N + 1 particles and the remaining N − 1 clusters contain one unit mass each.The range k max /2 ≤ k max is special: it may contain no more than a single cluster.This is simply because there is not enough mass in the distribution to have two clusters that are both larger than k max + 1.In the asymptotic limit M, N → ∞ at M/N = fixed the mass of a cluster in this range is of the same order as the total mass in the distribution.This means that the mass in the region k > k max is of the same order as that in k < k max .A cluster in k > k max represents a giant component, a single element of the population that carries a finite fraction of the total mass contained in the distribution.
The distributions that contain a giant cluster form a set that we call gel phase.The set of all other distributions forms the sol phase.Given an individual distribution n of the ensemble, a certain fraction of mass is contained in the sol region with the rest in the gel region.The ensemble averages of these fractions define the sol fraction φ sol and gel fraction φ gel in the ensemble: If P(n) is such that in the scaling limit φ gel → 0, the ensemble consists of a single phase, the sol, and is represented by the MPD in Eq. ( 40).If φ gel > 0 the ensemble is represented by a mixture of the two phases.We will determine their distributions and construct the tie line between the two phases.We suppose that the state at (M, N) consists of a sol phase with M sol , N sol = N − 1, and a gel phase with M gel = M − M sol .The evolution of the sol phase is governed by Eq. ( 31), which we now write as and interpret it to mean the the sol at all times satisfies this equation with θ sol = 1 − 1/ xsol .If the state is in the stable region in contains a single phase, the sol, with θ = M/N and q-β parameters from Eq. ( 87).
If the system is in the stable region 0 ≤ θ < θ * it consists of a single phase, the sol, with θ = 1 − N/M.In the unstable region the system contains a sol phase with M sol , N sol = N − 1 and a gel phase with 6.2.Phase Splitting -The Sol-Gel Transition M gel = −M sol , N gel = 1.The sol phase is obtained from Eq. ( 91) with θ sol = 1 − M sol /(N − 1) and its β-q parameters are given by Eq. ( 87) with θ = θ sol .The mass of the gel phase is then obtained from the conservation conditions M gel = M − M sol .This procedure is summarized below.IN all cases the overall state is described equivalently by (M, N), We detail below the phase behavior and the construction of the tie line in the two phase region.
Stable region 0 ≤ θ < θ * The system consists of a sol phase and its MPD is with Unstable region θ * ≤ θ < 1 The system consists of a sol phase with mass fraction φ sol and a gel phase with fraction φ gel = 1 − φ sol .
1. Obtain θ sol by solving with q(θ) from Eq. ( 87) and 2. Obtain φ sol and xsol as 3. Obtain the gel fraction from mass balance: The mean size of the gel cluster is φ gel M, where M is the total mass in the system.In the scaling limit the size of the gel cluster is ∞.
The evolution of the gel fraction and the mean cluster size for the product kernel (ν = 2) are shown in Fig .(5) as a function of θ.The gel fraction is zero up until the gel point (θ * = 0.5) and increases according to Eq. ( 96) once in the sol-gel regime.The mean cluster size increases in the pre-gel region but decreases in the post-gel region as clusters of the sol are lost by reaction with the gel cluster.At θ → 1 (t → ∞) all mass have is found in the gel phase except for a single sol particle with unit mass.This is is to be interpreted to mean that as long as there is a sol phase, mass is transferred from the sol to the gel.If we insist on the presence of the sol, its mass will be reduced to the smallest possible value, a single particle with until mass.
The evolution of xsol past the gel point retraces its pre-gel history.This is a consequence of Eq. ( 94) that resolves the sol phase (θ sol ) given the overall state of the system (θ).For θ > θ * the q value of the sol is equal to the q value of the overall system but the θ sol value of the sol must be in the stable region.The symmetry of q(θ) about θ * = 0.5 in the case of the product kernel produces a correspondingly symmetric evolution of xsol .The dashed lines are Monte Carlo simulations with M = 200 particles and are shown for comparison.The deviation from theory near the gel point is due to finite size effects.In these simulations a relatively small number of particles was used to permit the collection of a large number of realizations within reasonable computational time. 12./,0&%

Monte Carlo Simulations
The predictions of the theory can be placed into context via Monte Carlo (MC) simulation.In MC we generate an ensemble of independent trajectories via Markov transitions of the form starting from the same initial state (all monomers).By averaging the sampled distributions at fixed number of steps g = M − N we obtain the mean distribution n in generation g.The simulation begins with a distribution of M monomers and performs a random walk on the aggregation graph (Fig. 1) until it reaches a state with a single cluster.Therefore the simulation covers the entire range from θ = 0 to 1.
Figure ( 6) shows the mean distribution obtained by MC simulation with the product kernel using M = 200.Up until the gel point is reached the state is a single sol phase.It is characterized by a population of clusters whose tail decays fast enough that its moments are finite.Above the gel point a gel peak emerges.It becomes more pronounced and moves to larger sizes as aggregation progresses.The sol distribution contracts past the gel point.More specifically, it retraces its steps back to the monomeric state as θ increases.For example, the sol distribution at θ = 0.9 is identical to that at θ = 0.1 except that it carries less mass.In the Smoluchowski literature this is known as the Flory solution to gelation (Ziff and Stell, 1980).A competing solution by Stockmayer (Stockmayer, 1943) predicts that that the intensive distribution of the sol phase remains constant past the gel point except for the fact that its mass gradually decreases as it is transferred to the gel.As it turns out, Stockmayer solution implicitly assumes that P(n) is strictly a Gibbs distribution.In this case the sol-gel tie line is obtained by equating the temperatures of the two phases and the sol distribution is indeed found to be constant throughout the post gel region (Matsoukas, 2014).A discussion of the Stockmayer solution is beyond the scope of this work but a commentary is The gel phase emerges at θ * = 0.5 and moves towards ever larger sizes (arrows mark the theoretical predictions).The distribution of the sol grows in the pre-gel region range 0 < θ < 0.5 but contracts once past the post-gel point (θ > 0.5).

Continuous Limit
We define the continuous limit by the conditions Thus in addition to the scaling limit we require the mean cluster size to be much larger than the unit mass such that the cluster mass may be treated as a continuous variable, which we denote as x.Equations ( 64), ( 72) and (83) refer to this limit.We present the corresponding expressions for the partition function and the selection functional.
In the continuous domain all intensive properties of the ensemble are functions of the mean cluster size x.Thus we write β = β( x), q = q( x), w * i = w(x; x), and express the partition function in intensive form log ω( x) = (log Ω M,N )/N.We write the MPD as and satisfies the normalizations Since the microcanonical probability peaks sharply about the MPD (we are assuming a stable single-phase state) all ensemble averages revert to averages over the continuous MPD.Strictly, the lower limit in integrations over the MPD is 1 but in most cases it may be shifted to 0.
Here and in the results that follow we maintain the lower limit as 1 (the unit mass is the smallest mass in the ensemble) with the understanding that in most practical cases it can be shifted to 0. The cluster function w(x; x) is related to the functional derivative of the selection functional at the MPD: and the notation w(x; x) indicates this function of x will generally depend on x as well since the functional derivative of non linear functionals depend on the function on which the derivative is evaluated.The ensemble average kernel is equal to the mean kernel within the MPD The log of the intensive partition function lg These are the intensive forms of Eqs. ( 41) and ( 43), respectively.The partition function of aggregation is obtained from Eq. ( 33) by expressing the summation over log K g and an integral over K( x): The parameter β is obtained from Eq. ( 102) and log q from (101): (105) Finally we obtain the selection functional from Eq. ( 34) setting n → n * and converting the summations into integrals and write the result in the form where n * at x = M/N and n * ′ is its parent at x′ = M/(N − 1).The left-hand side is the change in log W( f ) upon d x = x − x′ , therefore the final result when x ≫ 1 is Equations ( 97), (103) ( 104), ( 105) and ( 107) provide an equivalent mathematical description of Smoluchowski aggregation in the continuous limit.These are accompanied by the variational condition which is the continuous form of Eq. ( 49) and is satisfied by all distributions p(x) with mean x.The equality defines the solution to the Smoluchowski process, the MPD, f (x).Special Case: Constant Kernel As a demonstration we apply these results to the constant kernel.The right-hand side of Eq. ( 107) is zero and we obtain W[ f ] = 1 at all times.It follows that w(x; x) = 1.The parameters β and q are β = 1/ x, q = x, (109) and the MPD becomes This is the well-known solution of the constant kernel in the continuous domain.The partition function of the constant kernel is and the inequality in Eq. ( 108) becomes whose right-hand side is the Shannon entropy.For fixed x, x > 0 it is maximized by the exponential distribution whose entropy is 1 + log x: the inequality is indeed satisfied.

Summary
With the results obtained here we have made contact with several previous works in the literature.The mean distribution for the constant kernel in Eq. ( 60) was given in (Hendriks, 1984).The same authors obtained a recursion for the partition function similar to that in Eq. ( 31) and presented a treatment with many 8. Summary common elements to ours but lacking the thermodynamic component of this work.The recursion for the cluster weights in Eq. ( 56) has appeared in treatments of aggregation, both deterministic (Lushnikov, 1973;Leyvraz, 2003) and stochastic (Spouge, 1983a,b;Hendriks, 1984).The mean distributions in the continuous limit for the mean and sum kernels and for the product kernel in the pre-gel region are well known results in the literature (Leyvraz, 2003).The instability of power-law kernels with ν ≥ 2 has been discussed by Ziff et al. (1982).All of these serve to validate the theory presented here and demonstrate that the thermodynamic treatment provides a unified theory of aggregation that brings previously disconnected results under a stochastic treatment we we call the Smoluchowski ensemble.The Smoluchowski ensemble is a probability space of distributions that are feasible under the rules of binary aggregation.The structure of this space, i.e., the connectivity of the graph in Fig.
(1), is solely determined by the condition that aggregation is a binary event; the probability measure over this space is determined by the rate expression prescribed by the aggregation model.In Smoluchowski aggregation the rate is directly proportional to the number of clusters that appear on the reactant side of the aggregation reaction and on the aggregation kernel.Under certain conditions that we review below the probability measure in the scaling limit is sharply peaked around one particular distribution of the ensemble, the most probable distribution (MPD).In this limit all ensemble averages reduce to averages of the MPD and this distribution alone suffices to produce all properties of the ensemble.The Smoluchowski coagulation equation is the time evolution of the most probable distribution in the asymptotic limit.
The step that turns the Smoluchowski ensemble into a thermodynamic ensemble is Eq. ( 19), which expresses the probability of distribution in the ensemble in terms of two special functionals, the multinomial coefficient n! and the selection functional W(n).This formulation introduces the partition function Ω M,N as the central property of the ensemble to which al other properties are connected.The thermodynamic calculus, summarized by the equations in Table 1, is a mathematical consequence of the variational condition that defines the most probable distribution in Eq. ( 40) as the solution to the constrained maximization of the probability P(n) in Eq. ( 19).The constraints are given by Eqs.(3) that fix the zeroth and first order moments of the distribution.These constraints define a microcanonical ensemble of distributions with fixed mean x = M/N.
The MPD obtained by the method of constrained maximization is stable provided that that partition function is concave in its independent variables.In extensive terms, log Ω M,N must be concave in M, N; in intensive terms, log ω( x) must be concave in x.The two conditions are equivalent and define the stability criterion of the MPD.As in regular thermodynamics, when the stability criterion is violated the system experiences phase splitting and exists a mixture of two phases -mathematically, as a linear combination of two independent MPD's.In aggregation these phases are the sol phase, which is represented by the MPD in Eq. ( 40) and the gel phase (giant component), which in the scaling limit is represented by a delta function at ∞.The splitting into a sol and gel phase is treated by the theory in a very natural and rigorous manner.
Notably, entropy in this treatment plays no special role.The Shannon entropy of distribution is the log of the multinomial coefficient.In the scaling limit, entropy is a component of the partition function through equation ( 37 In this form we have recovered the inequality of the second law as stated in in statistical thermodynamics: the entropy of the equilibrium distribution (MPD) is at maximum with respect to all feasible distributions, namely, all distributions with the same mean.As is well known this distribution is exponential.The constant kernel is special.With W(n) = 1 the probability of distribution is proportional to n!; accordingly, all ordered sequences of N clusters with total mass M are equally probable.The ordered sequence of Bibliography cluster masses in this case is analogous to microstate in statistical mechanics and the condition W = 1 analogous to the postulate if equal a priori probabilities.In the general case the Shannon entropy and the log of the microcanonical partition function are not the same.The fundamental functional that is maximized is the microcanonical weight n!W(n), whose log is The selection functional incorporates the effect of the aggregation kernel and in this sense it the point of contact between thermodynamics and the mathematical model of the stochastic process that gives rise to the probability space of interest.In Smoluchowski aggregation the model is defined by the transition rate in Eq. ( 11) and the corresponding governing equation for W is Eq. ( 34).
The thermodynamic formalism developed here is not limited to aggregation.Two alternative derivations that make no reference to stochastic process that gives rise to the probability space have been given in (Matsoukas, 2014) and (Matsoukas, 2019a).As long as log W is a homogeneous functional with degree 1, the thermodynamic relationships follow as a direct consequence of the maximization of the microcanonical probability in Eq. ( 19) under the constraints in Eq. ( 3).The details of aggregation enter through Eqs.Eq. ( 33) and ( 34) that give the the partition function and selection functional in terms of the aggregation kernel.The approach may be generalized to other processes including growth by monomer addition and breakup.These will be treated elsewhere.

Figure 1 :
Figure 1: The aggregation graph for M = 7.Each layer contains all feasible distributions in that generation.

Figure 3 :
Figure 3: Approach to the scaling limit for the product kernel in the pre-gel region with (a) x = 1.75 and (b) x = 4.The MPD is calculated from Eq. (72) and the mean distribution (dashed lines) from Eq. (79).The distributions for x = 4 are not stable and do not represent the correct state of the population.

Figure 4 :
Figure4: Phase diagram of power-law kernels: In the shaded region the system is stable and is represented by its MPD.The unshaded region is unstable and the system is split into two phases, a sol phase and a gel phase, each represented by its own MPD.Both graphs provide equivalent criteria of stability.

Figure 5 :
Figure 5: (a) Gel fraction and (b) mean sol cluster size as a function of the progress variable θ.Past the gel point the mean size in the sol retraces its pre-gel history back to its initial size xsol = 1.The dashed lines are MC simulations with M = 200 particles.

Figure 6 :
Figure 6: Monte Carlo snapshots of the mean distribution of the product kernel with M = 200 particles (open circles are MC results, solid lines are calculated from theory).The gel phase emerges at θ * = 0.5 and moves towards ever larger sizes (arrows mark the theoretical predictions).The distribution of the sol grows in the pre-gel region range 0 < θ < 0.5 but contracts once past the post-gel point (θ > 0.5).
), log Ω M,N = H(n * ) + log W(n * ),where H(n * ) is the Shannon functional evaluated at the MPD.In the special case of the constant kernel W(n * ) = 1.In this case the partition function reduces to the Shannon entropy of the MPD,log Ω M,N = −N ∑ ) ≤ H(n * ) = Ω M,N ; (constant kernel).

Table 1 :
Summary of thermodynamic relationships

Table 2 :
Summary of Constant, Sum and Product Kernel; in all cases θ