Morphology-dependent random binary fragmentation of in silico fractal-like agglomerates

Linear binary fragmentation of synthetic fractal-like agglomerates composed of spherical, equal-size, touching monomers is numerically investigated. Agglomerates of different morphologies are fragmented via random bond removal. The fragmentation algorithm relies on mapping each agglomerate onto an adjacency matrix. The numerically-determined fragment size distributions are U-shaped, clusters break predominantly into two largely dissimilar fragments, becoming more uniform as the fractal dimension decreases. A symmetric beta distribution reproduces the fragment distribution rather accurately. Its exponent depends on the structure (fractal dimension) and number of monomers of the initial agglomerate. A universal fragment distribution, a function only of the initial fractal dimension, is derived by requiring that it satisfy the fragmentation conversation laws and the straight-chain limit. We argue that the fragmentation rate is proportional to the initial agglomerate size.

Introduction. -Fragmentation is important in nature and industry. Polymer degradation [1], fracture of solids and volcanic eruptions [2], network topology and resilience [3], droplet [4] and agglomerate [5] breakup in turbulent flow are just a few examples where fragmentation determines the outcome of a process.
The fragmentation of linear chains has been extensively studied in polymer science since bond degradation or depolymerization breaks up the polymer. In most cases random scission [1,6,7], whereby all bonds break with equal probability, is considered. In colloids, the size distribution of colloidal aggregates in a fluid flow, and thence the suspension rheology, is greatly influenced by hydrodynamic shear-induced stresses that may break or restructure the particles [8].
In aerosol science, which plays a fundamental role in e.g., atmospheric sciences (cloud formation), air pollution (ultimate fate of atmospheric particulate pollutants), and industrial production of pharmaceuticals, the particle (a) E-mail: yannis.drossinos@ec.europa.eu (corresponding author) size distribution is one of the most important parameters that specifies the chemical and physical properties of the aerosol.
The dynamics of the size distribution, as determined by agglomeration, whereby the colliding particles retain their identity (or coagulation whereby colliding particle coalesce) and fragmentation, is usually described by meanfield rate equations. These population balance equations are generalizations of the Smoluchowski coagulation equation [9]. The kinetics of irreversible agglomeration has been extensively analysed, often with emphasis on scaling properties of the size distribution or the characteristic cluster size [10]. As clusters grow their number decreases continuously, their characteristic size increases, and the size distribution shrinks. In contrast, if fragmentation occurs, a likely event as the ever-growing clusters are more likely to fragment, a steady-state distribution may be established in the long-time limit.
Our primary interest in this work lies in the characterization of the size distribution of fragments (and the associated fragmentation kernel) upon random scission of p-1 arXiv:1909.12003v1 [cond-mat.soft] 26 Sep 2019 fractal-like agglomerates, and its dependence on the morphology of the initial agglomerate, as described by its fractal dimension (for fixed fractal prefactor). We consider numerically the linear, binary fragmentation of synthetic, loop-less, fractal-like agglomerates upon random link removal. Agglomerates that satisfy a scaling law over a finite range (fractal-like) are generated in silico. They are composed of equal-size, point-touching spherical monomers.
Discrete fragmentation kernel. -Consider the irreversible fragmentation equation [11], for the discrete particle size distribution n i (t), the number of agglomerates composed of i monomers (of size i) per unit volume at time t where a i is the (herein taken to be time-independent) fragmentation rate of a particle of size i, and p ij is the distribution of fragments of size i resulting from the break-up of a particle of initial size j. The maximum particle size i max is usually taken to be infinity. The fragment size distribution must obey a number of conservation laws [11,12]. The first constraint gives the expected number of fragments upon single-bond removal. For binary fragmentation it reduces to ∀j. (2) Since each particle breaks into two fragments, the fragment size distribution must be symmetric, p ij = p (j−i)j . The second conservation law is mass conservation: the sum of the number of monomers in all the fragments must equal the number of monomers in the initial agglomerate, In Appendix A (Supplementary Material) we show that for binary fragmentation eq. (2) implies eq. (3), rendering only one conservation law independent. In addition, eq. (2) suggests that g ij = p ij /2 defines a (discrete) probability distribution. Then, eq. (3) shows that the average fragment size is j/2 (as expected). The fragmentation rate a j is estimated by assuming that all bonds break with equal probability. Let the lifetime of each bond be characterized by a probability distribution, e.g., a Poisson distribution with average bond lifetime (characteristic time of fragmentation) τ . The number of bonds in an j-mer that would break in time t = τ is the number of monomers times the probability of a bond breaking. Thus, the fragmentation frequency is proportional to the initial size of the agglomerate, a j = j/τ . Continuous fragmentation kernel. -It is occasionally convenient to approximate the discrete particle size distribution by a continuous distribution. This approximation is very accurate for j 1. We also introduce a vector of morphological descriptors s to account for the morphology of the aggregate [13]. Accordingly, the discrete size distribution n i (t) is replaced by the continuous distribution n(x, t; s) where n(x, t; s)dx is the number of particles per unit volume in the size range x and x + dx and at morphological state s. For the cases considered herein, the additional morphological descriptor is the fractal dimension s = d f . The discrete variables i and j are replaced by the continuous variables x (fragment size) and y (initial agglomerate size), respectively. The continuous version of eq. (1) becomes where p(x, y; d f ) is the fragmentation density function to obtain a x-mer from the fragmentation of a y-mer. It is related to its discrete counterpart via p(x, y; d f )dx = p ij . The fragmentation frequency is a(x) = x/τ . The fragmentation density function must satisfy the conservation laws (see, also, Appendix A, Supplementary Material) with the symmetry condition p(x, y; d f ) = p(y − x, y; d f ).
As in the discrete case, only one conservation law is independent. The function p(x, y; d f ) is defined for 1 ≤ x ≤ y − 1. Henceforth, we will present results for the continuous fragment size distribution. The discrete case is summarized in Appendix C (Supplementary Material).
Generation and fragmentation of synthetic agglomerates. -We generated in silico independent, fractal-like agglomerates of specified structure (fixed fractal prefactor k f = 1.3 and varying dimension d f = 1.6, 1.8, 2.1) and mass (number of monomers j). We used the hierarchical tunable, off-lattice, cluster-cluster agglomeration algorithm, initially proposed by Thouy and Jullien [14] and later modified by Filippov et al. [15], to create mostly loop-less, agglomerates composed of equalsized, spherical, point-touching monomers. We only considered loop-less agglomerates. The algorithm does not aim to reproduce a physical agglomeration mechanism in that the motion of the colliding clusters is not explicitly simulated. Instead, it is based on geometric (including steric) considerations. The advantage of such a geometric algorithm is that the desired cluster morphology (d f , k f ) is an input (whereas in explicit cluster simulations it is an output). The synthetic agglomerates satisfy by construction exactly the fractal-like scaling law j = k f (R g /R 1 ) d f where R g is the radius of gyration of the cluster [16] and R 1 is the monomer radius. As summarized in Appendix B (table 1, Supplementary Material) the number of agglomerates varied from 400,000 to 25,000. The number of monomers per cluster varied from approximately 80 to 800.
According to the cluster generation algorithm (Appendix B, Supplementary Material,) clusters are created via monomer-monomer contact. If the resulting cluster does not have loops, the addition of a cluster onto an existing cluster generates only one bond. Thus, the total number of bonds in a loop-less j-cluster is j − 1; the sum of bonds connecting all monomers in a cluster [over both pairs (i, j − i) and (j − i, i)] is 2(j − 1). If f ij is the number of bonds that upon breaking would divide the initial cluster into two fragments of sizes i and j − i, then their sum is , as discussed extensively in Ref. [17]. The corresponding fragment distribution becomes p ij = f ij /(j − 1). Equivalently, these clusters have a (mean) coordination number of c j = 2(j − 1)/j [16,18].
The fragmentation algorithm is based on mapping the initial agglomerate onto a symmetric adjacency matrix. Each monomer configuration (i.e., cluster) has a unique representation as a graph via the adjacency matrix. A symmetric adjacency matrix A ij identifies completely a non-directed graph. The adjacency matrix is used to determine the connected components of a structure, thence the fragments upon removal of a bond. A similar approach was used by Isella and Drossinos [19] to identify clusters generated via monomer Langevin dynamics.
The adjacency matrix is constructed from all monomermonomer Euclidean distances: monomer-monomer links (bonds) are represented by one and their absence by a zero. Two monomers (i, j) are considered bonded if their center-of-mass distance D ij is smaller than the monomer diameter 2R 1 plus a threshold distance D thr , i.e., D ij < 2R 1 + D thr . We took D thr = 10 −3 × 2R 1 . A pictorial representation of a fragmentation event (referred to as "neck breakage") via adjacency matrices is shown in fig. 1. The fragmentation algorithm consists of randomly choosing a non-zero element of the adjacency matrix according to a discrete uniform probability distribution, ξ U (0, N link ), where N link is the number of links in the agglomerate. The selected element is replaced by a zero, thereby eliminating the monomer-monomer bond. A new adjacency matrix is formed, and its connected components are determined. If the initial agglomerate fragments, the sizes of the fragments are stored. If the aggregate does not fragment, it is removed from the list of initial clusters.
Only one link is removed per agglomerate, no multiple fragmentation events are considered. Then a new cluster is chosen, and the same procedure is repeated till the list of clusters has been exhausted. We found that fewer than 0.05% of the clusters did not fragment. By removing them we ensure that we consider only loop-less agglomerates. Figure 2 shows two simulated fragmentation events. On the left subfigure the resulting fragments are dissimilar, on the right they are of equal size. We found that the fragments have a slightly smaller fractal dimension, consequently their structure is more open, whereas they have a slightly higher fractal prefactor: they tend to be more locally compact than the initial agglomerate [16]. The structural parameters of the two fragments were not identical: they showed a slight asymmetry, which disappears in the fragment distribution.  . It is apparent that this fragmentation algorithm leads to largely dissimilar fragment sizes: the resulting fragment size distribution is U-shaped. Kalay and Ben-Naim [20] also found a Ushaped fragment distribution in their study of the fragmentation of random trees by removing a single node (a slightly different fragmentation process). A predominance of non-equal-sized fragments upon bond removal in random clusters on a lattice was also found in ref. [3]. The experimental results of Kusters et al. [21], who found that ultrasonic fragmentation of agglomerated particles in liquids leads to erosion, namely only a limited number of monomers break up from the agglomerate, also suggest a U-shaped fragment size distribution.
Empirical fragment size distribution. -We fitted the empirically determined fragment distribution to a symmetric beta distribution because it reproduces a Ushaped distribution for certain values of its parameter. Of course, our choice is not unique: other distributions could p-3 Y. Drossinos et al. reproduce our data, for example, eq. (11). However, the beta-distribution fit is more parsimonious: it depends on one parameter only [for fixed number of monomers in the agglomerates and fixed aggregate morphology Moreover, a symmetric beta distribution was introduced ad hoc previously to model the fragment size probability density function [22], see also ref. [23]. The functional form of the fragment density function p(x, y), expressed in terms of the relative number of monomers z = x/y, was chosen to be By construction the continuous probability density func- , and it satisfies the conservation laws Note that a fragment distribution proportional to [i −1 + (j − i) −1 ] α , a frequently used choice [24], would lead to a beta-like dependence proportional to [z(1 − z)] −α . For a given initial cluster size y, and a fractal dimension, the fragment distribution may be obtained by minimizing the distance, in the sense of the l 2 -norm, between the numerical results (the histograms shown in fig. 3 appropriately normalized) and the predictions of eqs. (6). The minimization procedure would minimize the distance between b(z)dz and p ij /2 for all i = 1, . . . j −1 to ensure that b(z)dz = p ij /2. For j − 1 fragments dz = 1/j, and the optimization condition becomes b(z) = jp ij /2. This procedure leads to a global optimization, i.e., an optimization over all fragment sizes.
Instead, we opted for a local optimization. It is more important to reproduce the small-fragment behaviour (z → 1/y or z → 1 − 1/y), the size of most fragments, than the complete distribution. Moreover, since the number of equal-size fragments is relatively small, the quality of the simulation data deteriorates as z → 1/2. We determined the probability of occurrence of a monomer fragment for all the clusters we generated for the three fractal dimensions. We found that provides an excellent fit of our simulation data. Equation (8) suggests that p 1y depends only on the fractal dimension and not on the size of the initial agglomerate. The exponent β was obtained by requiring that the probability of monomer occurrence equal the empirically determined monomer occurrence probability. We used eq. (8) to generate monomer occurrence probabilities for j = 10, 11, . . . , 1010, and we fitted the resulting beta-distribution exponent to stant asymptotic value as y → ∞. In that limit the fragment density function p(x, y; d f ) becomes approximately a homogeneous function under the scaling (x, y) → (λx, λy) since lim y→∞ β(d f , y) = β(d f ). The kernel does not become exactly homogeneous because some y-dependence remains in the integration limits. The y-dependent limits are required because the asymptotic exponent is β(d f , ∞) < −1 (for the studied fractal dimensions). The beta distribution is properly defined only for β > −1, the reason why our empirical fit eq. (6) is defined for z (0, 1). This is not a limitation since there is always a finite monomer size [23]. Homogeneous fragmentation kernels have been used extensively in the past, see, e.g., refs. [11,12,25,26]. We found that the coefficients presented in eq. (10b) depend slightly on whether a discrete or continuous distribution was modelled. The differences are small, but noticeable; see fig. 9 in Appendix C (Supplementary Material). Predictions of our empirical fit are compared to numerical simulations of single-bond, random fragmentation of DLCA agglomerates in fig. 5. Simulation results (filled square red symbols) were obtained by normalizing (to unity) the histogram of fragment sizes with j −1 bins. The beta distribution points, converted from the continuous density function to a discrete mass function, are denoted by filled blue pentagrams (note the logarithmic y-scale of the main figure; inset double logarithmic plot).
We also compare the numerical distributions to an expression proposed by Odriozola et al. [17], who analysed the fragment distribution of DLCA clusters generated by a slightly different algorithm. They fitted the fragment size distribution, specifically, f ij which is the number of bonds that upon breaking would divide a cluster into fragments of sizes i and j − i, to a four-parameter function, with p 1 = 0.4391, p 2 = 1.006, p 3 = 1.007, and p 4 = −0.1363. The product of the second term times the third term is reminiscent of the continuum Brownian agglomeration kernel for two fractal-like clusters of dimension d f . Note that p 2 ∼ p 3 ∼ 1, suggesting that eq. (11) approximates a symmetric beta distribution with a fixed exponent, independent of j, and for a fixed fractal dimension (d f = 1.8).
As argued earlier, eq. (11) may by normalized to unity by dividing it by the total number of bonds to obtain .
The predictions of eq. (12) are denoted by filled, green circles in figs. 5, 6, and 7. We note reasonable agreement of simulation data with theoretical predictions. It should be stressed that the normalization condition f ij is approximately equal to 2(j − 1). In fact, the authors of ref. [17] used this condition to estimate the accuracy of their fit. This implies that eq. (12) does not ensure mass conservation, rendering its use in population balance equations problematic. This is in sharp contrast to eqs. (6) that by construction satisfy exactly mass conservation.
The effect of morphology is studied in fig. 6. The two empirical fits are compared graphically to numerically determined fragment distributions arising from the break up of 192-monomer clusters for the three fractal dimensions. As the fractal dimension increases the probability of obtaining small fragments increases, the fragment distribution sharpens. As before the agreement is very good, even though eq. (11) was derived for relatively small DLCA clusters (j < 100).
We also present the probability of occurrence of a monomer fragment as a function of initial agglomerate size parametrized by the fractal dimension (of the initial agglomerate) in fig. 7. Since the exponent β was determined from the monomer probability, our proposed fit reproduces very well the simulation data, and it shows a dependence on the fractal dimension. As the fractal dimension decreases the probability of monomer occurrence decreases, In fact, we expect that as the fractal dimension decreases the fragment distribution would flatten out: the straight-chain limit (d f = 1) is a constant probability distribution independent of fragment size, but dependent on y, b(z; d f = 1) = 1, p(x, y; 1) = 2/y. parameter-free distribution) that depends only on the initial cluster fractal dimension. We constructed it by requiring that it satisfy exactly the conservation laws eqs. (7) and that it have the correct limit for straight chains. The universal fragment size probability distribution is see Eq. (13) It is easy to show that the conservation law is satisfied: function of z only and independent of y, under the scaling of sizes, b uni (z; d f ) remains invariant. Its homogeneity index is zero, leading to a homogeneous fragmentation kernel. The predictions of eq. (13) are compared to those of the beta distribution in fig. 8, top subfigure [for the asymptotic beta-distribution exponent β(d f , 1000)]. It is apparent that the parameter-free distribution does not reproduce the simulation data for d f = 1 as accurately as the beta distribution. However, it does reproduce the straight-chain limit, bottom subfigure. That subfigure shows that the effect of the fractal dimension becomes significant only for fractal dimensions very close to the singular limit d f → 1. The transition from a pronounced U-shaped distribution to the limiting case of a constant (flat) distribution of a straight line becomes noticeable for d f < 1.1 in that at d f = 1.1 the probability of a monomer fragment occurring is still predicted to be almost twice that of a monomer fragment of a straight chain.
Conclusions. -We studied numerically the fragment size distribution upon random scission of in silico fractallike agglomerates as a function of the initial size of the agglomerate and its structural properties, fractal dimension (d f = 1.6, 1.8, 2.1) and fixed prefactor (k f = 1.3). The fragment size distribution is an essential ingredient of the fragmentation kernel. Earlier studies, for example refs. [10,12,24,26], used fragmentation kernels in population balance equations whose analytical form was dictated either by physical arguments or by properties of related homogeneous agglomeration kernels, instead of originating from simulation data or experimental measurements. We followed the latter approach to obtain an analytical expression for the fragment size distribution that reproduces rather accurately our numerical simulations.
Our numerical simulations showed that, for our fragp-6 Random binary fragmentation of in silico fractal-like agglomerates mentation algorithm, the fragment size distributions are U-shaped, namely most clusters fragment into largely dissimilar fragments. We showed that a symmetric beta distribution reproduces rather accurately the empirical fragment size distributions. We found that the betadistribution exponent depends on the initial agglomerate morphology (via the fractal dimension for fixed prefactor) and the number of monomers y in the initial agglomerate, tending to a d f -dependent asymptotic limit for large y. We, also, derived a universal (parameter-free) fragment size distribution, dependent on the fractal dimension of the initial cluster, by requiring that it satisfy the two fragmentation conservation laws and that it reproduce the straight-chain limit. * * * The views expressed are purely those of the authors and may not in any circumstances be regarded as stating an official position of the European Commission. A.D. Melas was partially supported by the Horizon 2020 EU Framework Programme through the SUREAL-23 project (Grant Agreement 724136).

Supplementary Material
Appendix A: Fragmentation conservation laws.
-We first show that the mass conservation law, eq. (3) in the main text, for the discrete fragment size distribution is a direct consequence of the first conservation law [expected number of fragments, eq. (2)]. The derivation is based on splitting the sum into two parts, redefining a summation index, and then using the symmetry property of the fragment size distribution, p ij = p (j−i)j . Specifically, = j A similar decomposition of the first conservation law leads to Substitution of eq. (20) into eq. (19) gives The second sum in eq. (21) is zero because the factor j/2−i changes sign with respect to j/2, whereas the distribution p ij is symmetric. Specifically for an odd number of monomers in the cluster j = 2n + 1 the sum becomes since p (n+1)(2n+1) = p n(2n+1) . For j = 2n even, the relevant fragment sizes are n−(1/2) and n+(1/2). Since they are non-integers the distribution functions p (n−(1/2))(2n) and p (n+(1/2))(2n) are zero as fragments do not have noninteger number of monomers (alternatively, if half-integer fragments are considered, the previous argument for the sum being zero due to the symmetry properties of p ij still holds ensuring that the sum is zero even for an even number of monomers in the initial cluster). Hence, for the discrete fragment distribution, the first conservation law (expected number of fragments) along with the symmetry property of the fragment size distribution imply the second conservation law (mass conservation), The derivation applies mutatis mutandis to the continuous distribution p(x, y; d f ). Specifically, the mass conservation law becomes The integral on the RHS of eq. (23) vanishes because the function (y/2−x)p(x, y) is an odd function with respect to y/2. Therefore, as in the case of the discrete fragment size distribution, only one conservation law is independent.
Appendix B: Cluster generation and description. -The tunable cluster-cluster agglomeration algorithm [1,2] used to generate the synthetic agglomerates is hierarchical in that after the generation of the initial building blocks, the clusters combine pairwise. A new cluster is generated from two pre-existing clusters by first choosing randomly a sticking point and a sticking angle: the two clusters stick at that point with the appropriate orientation. The two random choices ensure that the resulting agglomerate is unique. Then, one of the initial clusters is rotated by randomly choosing the three Euler angles. A distance condition on the rotated cluster, which ensures that the resulting cluster would satisfy the scaling law, is checked. If satisfied, the code checks whether monomers overlap: if they do not, the new agglomerate is accepted [3].
Each pairwise binding defines a generation. We considered two types of initial block units: dimers (N init = 2) and a collection of k-mers randomly chosen to have between six and eight monomers, N init = 6, 7, 8. The number of monomers in a cluster of generation n is N init × 2 n . The random choice of initial building blocks of a varying number of monomers gives clusters containing a range of monomer numbers, centered about a mode (the most frequently occurring number of monomers in a set of clusters) that depends on the number of monomers in the initial building blocks (and the generation number). We tested two slightly different generation algorithms because we noted that random fragmentation of clusters generated only with dimers as the initial blocks tended to produce fragments containing a "magic" number of monomers (usually multiples of 2). Hence, fragmentation is a severe test of the cluster-generation algorithm. All the clusters considered in this work were generated with initial blocks of varying number of monomers.
We generate clusters with fixed prefactor (k f = 1.3) and variable fractal dimension d f = 1.6, 1.8, 2.1. The clusters we analysed belonged to generations n = 4 to n = 7, their number varying from 200,000 (400,000) to 25,000 (50,000) for d f = 1.6, 2.1 (d f = 1.8). At every generation the number of clusters decreases by 2. The number of monomers per cluster varied from approximately 80 to 800, depending on the generation number. The number of clusters considered, their size rage, and mode are shown in table 1.