Stationary mass distribution and nonlocality in models of coalescence and shattering

Many-body systems controlled by coalescence arise in many branches of science. The microscopic particles constituting such systems have a tendency to merge irreversibly upon collision or contact. Some examples include hydrogels for biomedical applications Berger et al. (2004), supramolecular polymer gels Sangeetha and Maitra (2005), aerosol formation Friedlander (2000), cloud formation Falkovich et al. (2002); Horvai et al. (2008), ductile fracture Pineau et al. (2016), and charged biopolymers Tom et al. (2016, 2017). By understanding the kinetics of coalescence, the macroscopic properties of such systems can be related to the microphysics of the fundamental collision and merging processes. For more applications and known results see the reviews Leyvraz (2003); Connaughton et al. (2010). In some applications, colliding particles may also fragment or shatter into smaller particles. Whether the collision between particles will result in coagulation or fragmentation of the constituent particles depends on the energy of the colliding particles Wada et al. (2009); Brilliantov et al. (2009). Typically particles with higher kinetic energy fragment, while slow moving ones coalesce or rebound. The size distribution of the fragmented particles is typically a power law distribution Nakamura and Fujiwara (1991); Giblin et al. (1998); Arakawa (1999); Kun et al. (2006), reproducible in simple tractable models Spahn et al. (2014); Dhar (2015). Such fragmentation processes find application in geophysics Grady and Lipkin (1980), astrophysics Michel et al. (2003); Nakamura et al. (2008), glacier modeling Åström et al. (2014), etc.

In this paper we are interested in situations where both coalescence and collisional fragmentation are simultaneously relevant. Examples include the fluctuations of phase coherent domains in high-temperature superconductors Johnson et al. (2011), the statistical properties of insurgent conflicts Bohorquez et al. (2009), the dynamics of herding behavior in financial markets Teh and Cheong (2016), and the formation and stability of planetary rings Davis et al. (1984); Longaretti (1989) where a coalescence–fragmentation model has recently been proposed to explain the particle size distribution of Saturn’s rings over several orders of magnitude Brilliantov et al. (2009, 2015). When coalescence and fragmentation occur together, one might expect the system to reach a nonequilibrium steady state in which the depletion of smaller particles due to coalescence is balanced by the depletion of larger particles due to fragmentation. These nonequilibrium states are expected to be insensitive to fine details of aggregation-fragmentation processes provided that the mass scales at which fragmentation acts as an effective source of light particles and the sink of heavy particles are widely separated. The simplest model of fragmentation for which the described scale separation occurs naturally is such that all fragmented particles are of the size of the smallest possible particle Brilliantov et al. (2009, 2015). We refer to such extreme fragmentation as shattering and use $m_{0}$ to denote the mass of the smallest particles generated. The expected universality of coalescence-shattering models explains the diversity of their applications and motivates a parametric study of how the particle size distribution depends on the form of the collision kernel, $K(m_{1},m_{2})$. The kernel depends on the nature of particle motion and interaction and gives the dependence of the microscopic collision rate on the masses, $m_{1}$ and $m_{2}$ of the colliding particles.

A particularly important class of collision kernels are homogeneous functions describing collisions that do not have a characteristic mass scale. This class includes many scientifically relevant cases, including the previously mentioned examples. We therefore restrict our attention to homogeneous kernels and denote the degree of homogeneity by $\beta$. Let us now consider how the stationary mass distribution should scale for a general homogeneous kernel. In the limit of small shattering probability $p$, there is a divergent mass scale $M(p)$ beyond which there are very few heavy particles due to a high cumulative probability of shattering. For a large range of masses $m_{0}\ll m\ll M(p)$, we expect $N(m)$ to be such that the flux of mass due to coalescence $J(m)$ is $m$ independent due to local mass conservation. In the well-mixed limit, the mass scaling of the flux can be determined from the following mean field scaling, (see Connaughton et al. (2004) for details):

Constant aggregation flux, $J$, implies that

By analogy with wave turbulence, we refer to the scaling exponent $(3+\beta)/2$ as the Kolmogorov-Zakharov or constant flux exponent. It depends only on the kernel homogeneity, $\beta$, and therefore possesses a high degree of universality. In particular, it will not change if we perturb the kernel while preserving the value of $\beta$ or change the nature of of sinks and sources (e. g., by removing heavy particles from the system once they become heavier than a fixed mass $M$ as in Connaughton et al. (2004) or removing colliding particles at a certain rate as in Connaughton et al. (2017)).

It is natural to ask when these constant flux solution is realized? This question can be answered by substituting the scaling solution as in Eq. (1) into the analytic integral expression for the flux $J$ and checking that it remains finite in the limit of small shattering probability $p$ (see Refs. Connaughton et al. (2004, 2008) for details of the derivation). It turns out that the realizability of constant flux solution depends on the locality of the collision mechanism. To characterize locality carefully, let us further reduce the class of kernels we are studying by assuming that in addition to having homogeneity degree $\beta$, $K(m_{1},m_{2})$ has the following asymptotic scaling when one of the colliding particles is much heavier than the other:

Clearly, $\beta=\mu+\nu$. This second reduction of generality is as natural for scale-free kernels as the homogeneity. It turns out that the realizabiliity of the Kolmogorov-Zakharov scaling depends, not on $\beta$, but on the difference

which we call the locality exponent. In Connaughton et al. (2004), we showed that for pure coalescence, the constant flux scaling (1) is realized if the locality exponent $\theta<1$. Physically, kernels with $\theta<1$ lead to the flux of mass due to coalescence being dominated by collisions between similar mass particles. Hence the term “locality.” For pure coalescence, the scaling of the particle size distribution for local kernels is strongly universal: when the characteristic scales of the source and sink are widely separated, it becomes independent of their details and depends only on the degree of homogeneity, $\beta$, of the collision kernel. We therefore expect intuitively that for the coalescence-fragmentation case with $\theta<1$, the scaling of the particle size distribution in the limit of small shattering probability $p$ is given by the constant flux expression (1).

This paper confirms this intuition and addresses the nonlocal case, $\theta>1$. We will use the simplest family of collision kernels labeled by two exponents $\nu$ and $\mu$ (or equivalently by $\beta=\nu+\mu$ and $\theta=\nu-\mu$):

Without loss of generality, we assume that $\nu\geq\mu$. This family has been widely used in general studies of aggregation (see Refs. Leyvraz (2003); Connaughton et al. (2010) for a review). Furthermore, we will always assume that the mean field approximation is applicable, which will allow us to calculate the mass distribution $N(m)$ by deriving and solving the corresponding Smoluchowski equation Smoluchowski (1917).

For the family of kernels (2) we will show that if $\theta>1$, the scaling exponent of the mass distribution in the limit of small shattering probability is both $\beta$- and $\theta$ dependent. Moreover, this dependence is different depending on whether $1<\theta<2$ (weak non-locality) or $\theta>2$ (strong non-locality). The amplitude of the mass distribution in the non-local regime is non-universal in the sense that it depends on the effective shattering scale and the mass of dust particles. It is interesting to note that these results for the aggregation-shattering model which conserves mass, parallel the answers for a non-conserved system with coalescence, input of small particles, and collision-dependent evaporation studied in  Connaughton et al. (2017). It seems that fine details of the mechanisms leading to effective sources and sinks of particles are irrelevant for a large class of coalescent models even in the non-local regime. An important conclusion from our analysis is that the mass distribution for the model kernel (2) for $\theta>1$ is different from the mass distribution for the local kernel $(m_{1}m_{2})^{\beta/2}$ with the same degree of homogeneity.

The paper is organized as follows. In Sec. II we define the model precisely and state our main quantitative results. In Sec. III we discuss the numerical algorithm that we use to solve for the steady state mass distribution. It is an iterative procedure that we show to reproduce known exact solutions. In Sec. IV we solve the model exactly for two special cases: first when $\mu=\nu$ ($\theta=0$) and second the addition model in which two particles coalesce only when at least one particle is a monomer ($\theta=\infty$). These exact solutions help us to benchmark the numerical algorithm. It is possible to obtain exact results for all integer $\theta$’s. This is discussed in Sec. V, where the presence of logarithmic corrections is established for some values of $\theta$. In Sec. VI, the small mass behavior of the mass distribution is studied using the exact relations between different moments. This enables us to determine the exponents when the kernel is local, and relations between the exponents when the kernel is non-local. In Sec. VII, we analyze the large mass behavior of the mass distribution by studying the singularities of the generating functions. By stitching together the small and large mass behavior, we are able to determine both the small and large mass asymptotic behavior of the mass distribution. In Sec. VIII we discuss the implications of our findings for the specific case of planetary rings since it is an interesting example where both local and nonlocal cases may be relevant. Finally, we conclude with a overview of results and directions of future research in Sec. IX.

Consider a collection of particles, each characterized by a single scalar parameter, mass. The mass of particle $i$ will be denoted by $m_{i}$, $i=1,2,,\ldots$, and will be measured in terms of the smallest possible mass in the system $m_{0}$, corresponding to the smallest possible dust particle, such that $m_{i}$ is an integer. Given a certain initial configuration, the system evolves in time via coagulation and collision-dependent fragmentation. Two particles of masses $m_{1}$ and $m_{2}$ collide at rate $(1+\lambda)K(m_{1},m_{2})$, where $K(m_{1},m_{2})$ is the collision kernel. On collision, with probability $1/(1+\lambda)$, the two particles coalesce to form a particle of mass $m_{1}+m_{2}$, and with probability $\lambda/(1+\lambda)$, fragment into $(m_{1}+m_{2})$ particles of mass $1$. Note that both the dynamical processes conserve mass, so that total mass is a constant of motion. We would be interested in the limiting case when the fragmentation rate tends to zero, i.e., $\lambda\to 0$. Also, we will be considering the well-mixed mean field limit when the spatial correlations between the particles may be neglected.

Let $N(m,t)$ denote the number of particles or mass $m$ per unit volume at time $t$. In the well-mixed dilute limit, the time evolution of $N(m,t)$ is described by the Smoluchowski equation:

The first term in the right hand side of Eq. (3) is a gain term that accounts for the number of ways a particle of mass $m$ may be created through a coalescence event. The second term is a loss term that accounts for the number of ways in which $N(m,t)$ decreases due to coalescence or fragmentation. The last term describes the creation of particles of mass $1$ due to fragmentation events. It is easy to check that the mean mass density is conserved. In this paper, we will be interested in the steady state solution of Eq. (3) obtained by setting the time derivative to $0$. We will denote the steady state solution by $N(m)$. $N(m)$ satisfies the equation

We consider the general class of kernels given by

The kernel may also be classified using two other exponents. The first is the homogeneity exponent $\beta$ defined through $K(hm_{1},hm_{2})=h^{\beta}K(m_{1},m_{2})$:

The second is the nonlocality exponent $\theta$ defined as

When $\beta>1$, the kernel is referred to as a gelling kernel and non-gelling otherwise. We will refer to kernels with $\theta<1$ as local kernels and non-local otherwise.

We also consider another kernel that corresponds to the so called addition model Hendriks and Ernst (1984); Brilliantov and Krapivsky (1991); Laurençot (1999); Ball et al. (2011); Blackman and Marshall (1994); Chávez et al. (1997) . Here collision events are allowed only if at least one of the particles has mass $1$. The kernel for the addition model is

which is characterized by a single exponent $\nu$. This kernel turns out to be exactly solvable (see Sec. IV.2).

In this paper, we will determine the asymptotic behavior of $N(m)$ through analysis of the moments as well as the singularities. Moments and generating function are defined as:

Clearly $F_{\alpha}(1)=\mathcal{M}_{\alpha}$. Multiplying Eq. (4) by $x^{m}$ and summing over all $m$, we obtain a relation between moments and generating functions,

We also define the exponents that characterize the mass distribution $N(m)$. We assume that the only relevant mass scale in the problem is the cutoff mass $M$ and hence $N(m)$ has the scaling form:

where $\tau$ is an exponent and $f(x)$ is a scaling function. $M$ denotes the cutoff scale below and above which $N(m)$ behaves differently. There are two cutoff mass scales in the problem. One is the total mass in the system and the other is the scale introduced by fragmentation. We will be working in the limit when total mass is infinite, but mean density is finite, leaving only one cutoff scale. The divergence of the cutoff mass scale as the fragmentation rate $\lambda\to 0$ is captured by

where the exponent $y$ will depend on the kernel. To characterize the scaling behavior for small and large masses, we introduce four new exponents $\tau_{s}$, $\eta_{s}$, $\tau_{\ell}$ and $\eta_{\ell}$ which are defined as:

The exponential decay for large mass is a conjecture. For $m\gg M$, the exponential decay with mass will be supported by exact solutions for special cases and numerical observation for more general cases. Further justification for arbitrary kernels follows from the additivity principle using which it has been argued that, for generic conserved mass models, the mass distribution has an exponential decay Chatterjee et al. (2014); Das et al. (2016). The four exponents are not independent. It is straightforward to obtain from Eq. (12) that

The results obtained in this paper for the different exponent as a function of the exponents $\theta$ and $\beta$ are summarized in Table. 1.

In this section, we describe the numerical scheme for obtaining the steady state mass distribution $N(m)$. Solving Eq. (4) in the steady state for $m=1$, we obtain

For $m\geq 2$, $N(m)$ may be determined from Eq. (4) in the steady state, provided ${\mathcal{M}}_{\mu}$, ${\mathcal{M}}_{\nu}$, and all $N(k)$ for $k<m$ are known:

Thus, ${\mathcal{M}}_{\mu}$ and ${\mathcal{M}}_{\nu}$ determine $N(m)$ for $m\geq 1$.

Consider scaled variables $\widetilde{N}_{m}=N(m)/N(1)$ and $\widetilde{{\mathcal{M}}}_{\alpha}={\mathcal{M}}_{\alpha}/N(1)$. In terms of these variables, Eqs. (17) and (18) reduce to

$\widetilde{{\mathcal{M}}}_{\mu}$ and $\widetilde{{\mathcal{M}}}_{\nu}$ determine $\widetilde{N}(m)$ for $m\geq 1$. The two unknowns, $\widetilde{{\mathcal{M}}}_{\mu}$ and $\widetilde{{\mathcal{M}}}_{\nu}$ are not independent and related to each other through Eq. (19).

To determine $\widetilde{{\mathcal{M}}}_{\mu}$, we follow an iterative procedure as summarized in the flowchart shown in Fig. 1. We start by assigning a numerical value (close to $1.0$) for $\widetilde{{\mathcal{M}}}_{\mu}$. $\widetilde{{\mathcal{M}}}_{\nu}$ is determined from Eq. (19). We then solve for $\widetilde{N}(m)$ up to a value of $m$ for which $\widetilde{N}(m)$ is larger than a predetermined value ($10^{-16}$ in our analysis), setting $N(m)=0$ for larger values of $m$. We then check for self consistency, i.e whether $\sum m^{\mu}\widetilde{N}(m)$ is equal to the preassigned value of $\widetilde{{\mathcal{M}}}_{\mu}$. We increment $\widetilde{{\mathcal{M}}}_{\mu}$ in small steps until the self-consistency condition is satisfied to the required precision. In our numerical analysis, we demand that the difference between the assumed and calculated values of $\widetilde{{\mathcal{M}}}_{\mu}$ should be smaller than $10^{-10}$.

To determine the unscaled variables $N(m)$, we use the fact that mass is conserved: $\sum_{1}^{\infty}mN(m)=\rho$, where $\rho$ will be treated as a parameter. We then scale all $\widetilde{N}(m)$ by the same factor so that the desired mass density $\rho$ is achieved, thereby determining $N(m)$. In our numerical measurements, we set $\rho=1$. There is no proof that the algorithm will result in the convergence of $N(m)$ to its correct value. However, we verify the convergence for special solvable kernels (see Sec. IV), leading us to expect that the mass distribution converges to its correct value for more general kernels.

From the numerically computed $N(m)$, we observe that, for all values of $\mu$ and $\nu$ that we have studied, $N(m)$ decays exponentially at large masses [as in Eq. (15)]. The exponential cutoff mass $M$ is determined by solving for the three parameters ($\tau_{\ell}$, $M$, and $a_{\ell}/M^{\eta_{\ell}}$) in Eq. (15) using $N(m)$ for three consecutive $m$’s and extrapolating to large $m$. Once $M$ is determined, the compensated mass distribution $N(m)e^{m/M}$ is a power law with exponent $\tau_{\ell}$, allowing us to verify the theoretical results for the exponents at large mass.

The steady state mass distribution $N(m)$ may be determined exactly for two cases: (1) when $\mu=\nu=\beta/2$ and (2) the addition model (defined in Sec. IV.2).

It is possible to obtain exact results for the case when the locality exponent $\theta$ is an integer, namely $\theta=n$, $n$ is an integer. The starting point is Eq. (II). If $\nu=\mu+n$, where $n=1,2,\ldots$, Eq. (II) reduces to a closed differential equation for $F_{\mu}$:

We expect the singularities to occur at the points in the complex $x$-plane where the coefficient in front of the highest order term is zero. Therefore, at the singular point $x_{c}$, $F_{\mu}$ satisfies

Introduce new variables $f(x)$ as follows:

where $f(x_{c})=0$ and $t_{c}=\ln(x_{c})>0$. Then Eq. (V) may be rewritten as

Equation (39) becomes more tractable under the following transformations:

such that we obtain

where

We note that $g(\tau_{c})=0$. We now analyze Eq. (42) for specific integer values of $\theta$.

An exact solution is possible only when $\theta=0$. In this section, we use moment analysis to determine some of the exponents characterizing $N(m)$ for general $\theta$ and $\beta$. In particular, we study the small mass behavior of the mass distribution $N(m)$ as described by Eq. (14). Our aim is to determine the exponents $\tau_{s}$, $\eta_{s}$, and $y$ as a function of $\beta$ and $\theta$. For this, we will require the equations satisfied by the different moments of $m$. These may be obtained by differentiating Eq. (II) with respect to $x$ or by multiplying Eq. (4) by $m^{n}$ and summing over $m$ from $1$ to $\infty$. Doing so gives