Stochastic Theory of Discrete Binary Fragmentation—Kinetics and Thermodynamics

We formulate binary fragmentation as a discrete stochastic process in which an integer mass k splits into two integer fragments j, k−j, with rate proportional to the fragmentation kernel Fj,k−j. We construct the ensemble of all distributions that can form in fixed number of steps from initial mass M and obtain their probabilities in terms of the fragmentation kernel. We obtain its partition function, the mean distribution and its evolution in time, and determine its stability using standard thermodynamic tools. We show that shattering is a phase transition that takes place when the stability conditions of the partition function are violated. We further discuss the close analogy between shattering and gelation, and between fragmentation and aggregation in general.


Introduction
Binary fragmentation is the splitting of a cluster into two. It is the inverse of binary aggregation, in which two clusters merge into a new one. Aggregation and fragmentation are two of the most common mechanisms in population balances. Both processes are irreversible, but in opposite directions: to ever increasing sizes in aggregation, and ever decreasing sizes in breakup. In addition, both processes are capable of exhibiting a remarkable dynamic behavior manifested in the breakdown of mean-field models that describe them. In aggregation, this behavior is known as gelation, and is characterized by the emergence of a giant cluster with infinite mass (the "gel") that accounts for a finite fraction of the total mass [1][2][3][4]. An analogous behavior is known to occur in fragmentation in the form of a ghost phase of clusters with zero mass ("dust") that contain a finite fraction of the total mass [5,6]. This process of "shattering" has evoked analogies to phase transitions and motivated numerous investigations of the kinetic equations of fragmentation and the conditions that may result to this unique behavior [7][8][9][10][11][12][13][14][15].
The usual approach to the study of fragmentation is through the fragmentation equation, an integro-differential equation for the mean cluster distribution [16]. In reality, fragmentation is a discrete stochastic processes: a discrete cluster disintegrates into discrete fragments. The distribution of clusters itself is a stochastic variable, whose evolution in time represents one of many possible trajectories in the space of feasible distributions. The classical treatment presumes that the solution to the fragmentation equation is the same distribution that would be obtained by averaging over all trajectories. Whether this assumption is correct or not requires a formal treatment of the stochastic problem, but one case where this is clearly not true is shattering. At its most elementary level, fragmentation is mass conserving; the loss of mass conservation in the kinetic equation is indication that this equation no longer tracks the true mean distribution. In this paper, we formulate fragmentation as a discrete stochastic process. The distribution of clusters is itself a random variable and transitions in time. The set of all possible distributions that can be reached in a fixed number of fragmentation events forms the fragmentation ensemble. The goal of the theory is to assign probabilities to the distributions of the ensemble based on the breakup law. Once these probabilities are assigned, the mean distribution may be obtained by direct averaging.
This approach is most general and encompasses the classical kinetic treatment as a special result. It has the further advantage that it makes contact with statistical mechanics. The set of feasible of distributions is a thermodynamic ensemble: its most probable distribution in the asymptotic limit is overwhelmingly more probable than all others and maximizes a functional similar to entropy. This connection offers the tools to treat shattering as a formal phase transition using the stability criteria of phase equilibrium. We have previously applied this approach to binary aggregation and shown that gelation is a formal phase transition in the thermodynamic sense [4,[17][18][19]. Here, we formulate the statistical mechanics of the fragmentation ensemble and show that shattering is a phase transition mathematically symmetric to gelation.
The paper is organized as follows. In Section 2, we summarize the elements of the ensemble theory, which forms the basis of the stochastic treatment. In Section 3, we construct the ensembles of discrete binary fragmentation and obtain the kinetic equations that relate transitions between generations. In Section 4, we apply the theory to random fragmentation and obtain the complete solution in closed form. We discuss shattering in Section 5, using stability analysis, and finally, we summarize the results in Section 6.

Statistical Thermodynamics of the Cluster Ensemble
We begin with a brief summary of the thermodynamic theory that forms the basis of the study. A detailed exposition is given elsewhere [17,20]. We consider a discrete distribution of clusters n = (n 1 , n 2 · · · ), where n i is the number of clusters with mass i = 1, 2 · · · . The number of cluster in n is N and their total mass is M: We construct the ensemble E M,N of every possible distribution that satisfies these conditions; this ensemble contains all possible ways to partition integer mass M into N clusters. We assign probabilities to the distributions of the ensemble according to where n! is the multinomial coefficient, W(n) is a functional of distribution n that embodies the details of the stochastic process and which will remain unspecified until we consider specific models of fragmentation; and Ω M,N is the partition function, a factor that ensures proper normalization of probabilities, We refer to P(n) as microcanonical probability and to Ω M,N as its microcanonical partition function. The most probable distribution n * in the ensemble is obtained by maximizing log P(n) with respect to n under the two constraints in Equation (1). The result of this maximization is with β (inverse "temperature"), q ("canonical" partition function) and w * k ("multiplicity" of cluster size k) given by In the asymptotic limit M > N 1 the most probable distribution becomes overwhelmingly more probable. All properties of the ensemble in the asymptotic limit reduce to the corresponding property of the most probable distribution. The fact that n * maximizes the microcanonical probability leads to the following inequality where H(n) is the Shannon entropy of distribution The inequality in Equation (7) is a generalized statement of the second law: it applies to all distributions in the ensemble and becomes an exact equality only for n = n * . Equations (2)-(7) provide a treatment in the language of statistical mechanics of the ensemble of distributions defined by the two constraints in Equation (1). Next, we apply this formulation to stochastic fragmentation.

Fragmentation Kernel
In binary fragmentation, mass k breaks into a pair of fragments with masses k − j and j. All masses are larger than or equal to 1 (monomers do not break any further). The process can be represented in the form of a reaction (k) where F k−j,j is the rate of formation of ordered pair of fragments (k − j, j) and satisfies the symmetry condition F j,k−j = F k−j,j . A number of related functions may be defined based on this kernel. The breakage rate of mass k is and the number of fragments with size j produced from particle mass k is The mean fragmentation kernel within distribution n is the mean value of F i,j among all breakup events in the distribution: where M − N is the total number of fragmentation events in distribution n. The ensemble average fragmentation kernel in ensemble E M,N is the mean value of the fragmentation kernel F i,j over all possible fragmentation events in the ensemble: where P(n) is the probability of distribution and the outer summation is over all distributions in ensemble E M,N .

Fragmentation Ensemble
To construct the ensemble of distributions in fragmentation we begin with a single particle with mass M and subject it to g fragmentation events, an event consisting of replacing a cluster with two fragments. After g events the distribution contains N = g + 1 elements (clusters). The set of all distributions that can be reached in g fragmentation events forms the fragmentation ensemble E M,N . All distribution in the ensemble satisfies the two constraints in Equation (1) and, conversely, all distributions that satisfy these constraints can be formed by binary breakage of mass M in g = N − 1 fragmentation steps.
We represent this process pictorially in the form of a layered graph. Figure 1 illustrates this graph for M = 7. Nodes represent fragment distributions and edges indicate transitions between distributions following a fragmentation event. The graph is layered the number of fragmentation events to reach a distribution and each layer contains the complete ensemble of distributions that can be reached in fixed number of fragmentation events. The graph is the exact reverse of that of aggregation: by reversing the direction of all arrows we obtain the graph for discrete binary aggregation of M monomers (bottom) advancing towards a single cluster (top) through M − 1 binary aggregation events [19].

Transitions between Distributions
When a cluster with mass i in distribution n of ensemble E M,N−1 breaks up into fragments i − j and j the result is a distribution n in the ensemble of the next generation, E M,N ; schematically, The probability of this transition is proportional to the number n i of clusters with mass i in the parent distribution and to the fragmentation rate F i−j,j : Here, F M,N−1 is the mean fragmentation kernel over all fragmentation events in the parent ensemble: where P(n ) is the probability of distribution in the parent ensemble. There are M − N + 1 ordered pairs of fragments in any distribution that contains N − 1 clusters with total mass M. Accordingly, (M − N + 1) F M,N−1 , which appears as the normalizing constant in the transition probability in Equation (15), is the total fragmentation rate in parent ensemble E M,N−1 .
The probability of distribution propagates from one generation to the next with transition probability P(n → n) according to the Master Equation, with the summation taken over all parents n of distribution n; in terms of Figure 1 the summation is over all edges of the fragmentation graph that lead to n. At the initial state distribution n 0 consists of a single cluster with mass M. In this state, we have With Equation (18) To obtain this result, we used the fact that the parent-offspring relationship is reversible (Figure 1): if parent distribution n produces offspring n via the fragmentation event (i) → (i − j) + (j), this parent can be constructed by applying the aggregation reaction (i − j) + (j) → (i) to the offspring distribution. In view of this observation, each element of the double summation in Equation (19) is a fragmentation parent of n that contributes towards the offspring distribution with rate F i−j,j . Accordingly, the double summation goes over all parents of offspring distribution n. Since the left-hand side is independent of distribution n, so must be the right-hand side. This can be true only if the quantity in braces is independent of distribution n, We fix Λ by imposing a normalization condition on the sampling functional: if W(n) is constant for all distributions, we require W(n) = 1. We will confirm in Section 4 that if F i,j = 1 then W(n) = const. This implies Λ = 1. Accordingly Equation (19) separates into two independent recursions, one for the partition function and one for the sampling functional. For the partition function, we have which is readily solved for Ω M,N starting with Ω M,1 = 1 in generation 0. The result is The corresponding expression for the sampling functional is which gives the sampling functional of distribution in terms of the sampling functional of its parents. The parameters β and q are and are obtained by applying the finite-difference form of the derivatives in Equation (6) to the partition function in Equation (22).

Kinetics and Thermodynamics
We may now summarize the key results of the stochastic theory of fragmentation. Equations (22) and (23) (5), are to be obtained from the derivatives of the sampling functional according to Equation (6). In the context of this theory, the mathematical problem reduces to the calculation of Ω and log W. We demonstrate this procedure in Section 4 with an example that produces analytic results for the partition function and the sampling functional.

Mean Distribution
The mean number of clusters with mass k is where n k is the kth element of distribution n, P(n) is the probability of distribution, and the summation is over all n in the ensemble. Using P(n) from Equation (17) and the stoichiometry of the fragmentation reaction we obtain Here primed terms refer to parent elements of element n k and the ensemble average on the right-hand side is over the parent ensemble E N,N−1 . Lastly, we use Equation (15) for the transition probability P and perform the summations over the Kronecker deltas. The final result is This is the governing equation for the evolution of the mean distribution: on the left-hand side, we have the change in the mean number of cluster of size k between generations; on the right-hand side, we have the average over all trajectories emanating from the the parent generation. If the ensemble is represented by a single distribution, the mean distributionn, the result can be written as as difference equation forn k with the ensemble averages dropped: Here, ∆N = 1 is the change in the number of clusters between successive generations and F is the mean kernel over all fragmentation events inn. To make full connection with the classical treatment we write the rate equation for the number of clusters N. Each binary fragmentation event produces one net new cluster, therefore, the rate at which the parent generation produces clusters is Dividing Equation (29) by (30) we obtain We recognize the result as the classical fragmentation equation for the mean distribution. It is the appropriate limit of of Equation (28) if the probability of distribution, P(n), is sharply peaked around a single distribution. In this case, and in this only, all average quantities over the ensemble can be replaced by the same quantity evaluated over a single distribution, the most probable distribution of its generation, also equal to the mean distribution in the ensemble.

Special Case: Random Fragmentation (F i−j,j = 1)
The constant kernel, F i−j,j = 1, represents random fragmentation, a process in which all possible fragmentation events are equally probable. The fragmentation rate is a k = k − 1, equal to the number of possible events within mass k, and the fragment distribution is uniform, b(i|k) = 2/(k − 1). For this model we obtain the full solution of the stochastic problem in closed form.
The mean kernel in any distribution n is F(n) = 1 and it follows that F M,N = 1 for all M, N. The partition function follows from Equation (22) and is equal to the number of ways to break mass M into N ordered fragments [21]. To obtain the sampling functional we note the identity We conclude from Equation (23) that W(n) must be the same for all distributions in all generations, and since W(n 0 ) in generation zero, we obtain with w i = 1. This is a special case of functionals that factorize in the form Such systems, also known as Gibbs ensembles, appear in various probabilistic contexts and their properties are well known [22][23][24]. In particular, the mean cluster distribution is given in terms of the partition function and the factors w k by [17] The derivatives β and q are and the most probable distribution is exponential: wherex = M/N is the mean cluster size. It is easy to show that the mean distribution in Equation (36) and the most probable distribution converge to each other. Equation (36) for the mean fragment distribution is exact for all integer masses M > N and fragment masses M − N + 1 ≥ k ≥ 1. The first few fragment distributions for general M are shown in Table 1. For N = 2 the distribution is uniform; for N = 3 it is linear in fragment size; in general the distribution is a polynomial of degree N − 1 and always monotonically decaying in the size range (0, M). It is interesting to point out that these distributions are the same as those produced by a single random fragmentation event into N pieces (see Equation (12) in [21]). That is, g consecutive random binary fragmentation events produce the same ensemble of fragments as as a single random fragmentation event into g − 1 pieces. These distributions are also the same as those in discrete binary aggregation with constant aggregation kernel K i,j = 1 [19]. That is, fragmentation with F i,j = 1 and aggregation with K i,j = 1 both assign the same probabilities to all distributions int he event space, even though the two processes move in opposite directions with different transition probabilities each. This symmetry between fragmentation and aggregation has not been recognized previously.
As a demonstration of the exactness of Equation (36) we perform a Monte Carlo simulation of random fragmentation with M = 50 (Figure 2). The agreement between simulation and the theoretical distribution is excellent. (Details on the Monte Carlo simulation are given in Section 5.3).  Table 1. Table 1. Mean distributions of random binary fragmentation with constant kernel F i,j = 1 for N = 2, 3, 4 and 5.

Power-Law Breakage
A special family of fragmentation kernels is of the power-law form where ν is a constant. This kernel represents binary breakage with power-law rate a k = (k − 1) ν+1 and uniform fragment distribution b(i|k) = 2/(k − 1). The power-law kernel has been studied extensively, especially since it leads to shattering when the exponent is negative [5][6][7][8]16]. From a kinetic perspective, shattering is the result of the accelerated fragmentation rate with decreasing size that produces an accumulation at the low end of the distribution. Here, we examine shattering under the criteria for thermodynamic stability. In molecular thermodynamics phase splitting occurs when the concavity of the microcanonical partition function is violated. The same condition identifies the presence of shattering, as we will see. The mean fragmentation rate of the power-law kernel in distribution n is when M N 1. Thus we obtain the scaling We may now obtain the partition function. Using Equation (22) we find Its parameters β and q are obtained from the derivatives of the partition function according to Equation (6). The result can be expressed as where θ is a scaled size variable, defined as The scaled variable maps the complete evolution of the system into the interval 0 ≤ θ ≤ 1: Fragmentation commences at θ = 1 − 1/M ∼ 1 and ends at θ = 0.

Stability and Phase Transitions
Stability requires log Ω M,N to be concave function of M and N. The condition ensures that the most probable distribution is the solution to ∂ log P(n * )/∂n * = 0 and is given by Equation (5). Concavity requires the second derivatives of log Ω M,N to be negative. In terms of the scaled variable θ the conditions reduce to the concave conditions are Using Equations (44) and (45), both inequalities in Equation (47) reduce to the same result: This condition guarantees concavity of the partition function. In the language of thermodynamics, it guarantees the stability of a single-phase system; if it is violated the system splits into two distinct coexisting phases. The stability maps of power-law fragmentation are shown in Figure 3 on the (β, θ) and (q, θ) planes. For ν > 0 the system is always stable and is represented by a single phase, a distribution of clusters given by Equation (5). When ν is negative stability depends the critical value θ * = 1/(1 − ν) falls within the interval (0, 1) and splits it into two regions. Below θ * the system is stable, above θ * is not. As fragmentation always starts at θ = 1, the system begins in the unstable regime: shattering is observed at all times. Figure 3. Stability map of power-law kernel in terms of the parameters β (left) and q (right). Stable branches are shown blue, unstable in red.

Monte Carlo Simulations of Shattering
The most direct way to illustrate the shattering instability is through simulation. We use Monte Carlo to sample trajectories in the phase space of fragmentation. The simulations are conducted as follows. Given N clusters we choose a cluster to break with probability proportional to (k − 1) ν , where k is the mass of the cluster. Mass k is then split into integer fragments i and k − i, which are chosen with equal probability among all possible k − 1 fragments. The simulation begins with a single cluster with mass M, which is then subjected to M − 1 fragmentation events until it forms M fragments with unit mass. This constitutes one trajectory that tracks the distribution of fragments from generation g = 0 to generation g = M − 1. The mean distribution is calculated by averaging 5000 trajectories. Figure 4 shows results for three different masses, M = 50, 100 and 200, to highlight the convergence to asymptotic behavior. Cluster distributions are plotted against the scaled mass z = i/M. In the limit M → ∞ this scaling amounts to breaking a particle of unit mass into fragments in the continuous interval from 0 and 1. For ν = 1 (random fragmentation), the system is stable at all times and its distribution is given by Equation (39). The results of the simulation in Figure 4a-c For ν = −3, the system is unstable and the distribution of clusters contains a new feature: the number of the smallest fragment size (unit mass) deviates markedly above the distribution of the rest of the fragments. This feature persists: it does not go away as M is increased, neither does is subside as fragmentation goes on. In fact, the divergence of this segment of the population becomes more pronounced as M → ∞. In this limit the divergent mass is z * = 1/M → 0: accordingly, a finite fraction of the mass of the system is accumulated into fragments of zero size.
McGrady and Ziff [6] pointed out the close connection between shattering and gelation. Aggregation with power-law kernel K i,j = (ij) α/2 is analogous to the power-law fragmentation kernel F i,j = (i + j) ν : it produces stable solutions for α ≤ 1 but leads to gelation for α > 1. This connection is not merely qualitative. The parameters β and q for power-law aggregation were obtained by [19]. We compare them to those for power-law fragmentation ( Table 2). With the substitution α = 1 − ν, the parameter q and the region of stability are the same in both processes. Gelation occurs when a cluster appears in the giant region (M − N + 1)/2 + 1 < kM − N + 1, a region defined by the condition that it may contain at most one cluster [4]. In the scaling limit of aggregation (M > N → ∞ at fixed M/N) the giant cluster goes over to infinity and disappears from the size distribution, giving rise to an apparent loss of mass and the breakdown of mass conservation. The giant cluster is also present in shattering. In terms of the scaled mass z = i/M, the giant cluster in fragmentation is asymptotically of the order 1. Under shattering conditions (Figure 4d-f) the neighborhood of z ≈ 1 remains populated at all times. By contrast, stable solutions are devoid of clusters in this region (Figure 4a-c). The presence of a giant cluster is a universal feature of instability, whether the process is aggregation or fragmentation. Table 2. Comparison between power-law fragmentation and power-law aggregation.

Fragmentation
Aggregation Fundamentally shattering and gelation represent the same condition: the presence of a single cluster in the giant region that coexists with a distribution of clusters in the sub-giant region. It is the different scaling between fragmentation and aggregation that sends the ghost phase to zero size (shattering) or infinite size (gelation), where it becomes invisible either in the form of an infinite concentration of clusters with zero mass (shattering) or as a zero concentration of particles with infinite mass. There is, however, one difference: in gelation the system begins in the stable branch and within finite time crosses into the unstable regime. As we have shown previously [4,19] in this case we can construct the tie line between the two phases by standard thermodynamic arguments to obtain the distribution of the dispersed phase and the gel fraction at all times. In shattering, the system begins in the unstable region and never crosses over to the stable regime. Though we believe that in this case also it ought to be possible to construct the tie line, we were not able to do it here. This requires a more detailed treatment of the unstable partition function, in particular, the asymptotic scaling of the mean fragmentation kernel in Equation (42) when the exponent is negative. In discrete finite systems instability cannot occur instantaneously. The initial state is stable and it must take a finite number of fragmentation events to reach instability. The scaling form of the mean kernel in Equation (42), while sufficient to map the stability boundary, it is not of sufficient accuracy to study the onset of phase splitting.

Conflicts of Interest:
The author declares no conflict of interest.