Instantaneous Gelation in Smoluchowski’s Coagulation Equation Revisited

# Instantaneous Gelation in Smoluchowski’s Coagulation Equation Revisited

Robin C. Ball Centre for Complexity Science, University of Warwick, Gibbet Hill Road, Coventry CV4 7AL, UK Department of Physics, University of Warwick, Gibbet Hill Road, Coventry CV4 7AL, UK    Colm Connaughton Centre for Complexity Science, University of Warwick, Gibbet Hill Road, Coventry CV4 7AL, UK Mathematics Institute, University of Warwick, Gibbet Hill Road, Coventry CV4 7AL, UK    Thorwald H.M. Stein Department of Meteorology, University of Reading, Earley Gate, Reading, RG6 6BB, UK    Oleg Zaboronski Mathematics Institute, University of Warwick, Gibbet Hill Road, Coventry CV4 7AL, UK
July 7, 2019
###### Abstract

We study the solutions of the Smoluchowski coagulation equation with a regularisation term which removes clusters from the system when their mass exceeds a specified cut-off size, . We focus primarily on collision kernels which would exhibit an instantaneous gelation transition in the absence of any regularisation. Numerical simulations demonstrate that for such kernels with monodisperse initial data, the regularised gelation time decreases as increases, consistent with the expectation that the gelation time is zero in the unregularised system. This decrease appears to be a logarithmically slow function of , indicating that instantaneously gelling kernels may still be justifiable as physical models despite the fact that they are highly singular in the absence of a cut-off. We also study the case when a source of monomers is introduced in the regularised system. In this case a stationary state is reached. We present a complete analytic description of this regularised stationary state for the model kernel, , which gels instantaneously when if . The stationary cluster size distribution decays as a stretched exponential for small cluster sizes and crosses over to a power law decay with exponent for large cluster sizes. The total particle density in the stationary state slowly vanishes as when . The approach to the stationary state is non-trivial : oscillations about the stationary state emerge from the interplay between the monomer injection and the cut-off, , which decay very slowly when is large. A quantitative analysis of these oscillations is provided for the addition model which describes the situation in which clusters can only grow by absorbing monomers.

83.80.Jx

## I Introduction

The determination of the statistical time evolution of an ensemble of particles undergoing irreversible binary coagulation is a problem which arises in many different contexts in the physical, chemical and biological sciences. See ERN1986 () for a list of applications. Of pre-eminent interest is the cluster size distribution, , which specifies the average density of particles of mass, , at a given time . Often one aims to derive from a given model of the microscopic dynamics governing cluster coagulation. This problem has been extensively studied for almost a century starting with the seminal work of Smoluchowski SMO1917 () who showed that, at the mean field level, the cluster size distribution of a statistically homogeneous system evolves according to the Smoluchowski coagulation equation:

 ∂tNm = ∫m0dm1K(m1,m−m1)Nm1Nm−m1 − 2Nm∫∞0dm1K(m,m1)Nm1+Jm0 δ(m−m0).

Here the coagulation kernel, , is proportional to the probability rate for a cluster of mass merging with a cluster of mass . At this level of description, it encodes all relevant information about the underlying micro-physics. is the rate of injection of monomers, having mass , into the system. may be zero depending on the application. In this paper, we take Eq. (I) as our departure point. It is expected to apply when the clusters remain well-mixed (despite aggregation), but the elucidation of the exact conditions under which Eq. (I) is rigorously obtained as the mean field limit of an underlying stochastic process can be a subtle question. Mathematically inclined readers are referred to the review by Aldous ALD1999 () for some common mathematical perspectives on this issue. A physical example of a situation in which the applicability of Eq. (I) fails due to the generation of spatial correlations between particles by diffusive fluctuations is discussed in detail in CRZ2006 ().

In many applications the microphysics is scale invariant, at least over some range of cluster sizes. That is to say, the coagulation kernel is a homogeneous function of its arguments, the degree of homogeneity of which we shall denote by :

 K(cm1,cm2)=cλK(m1,m2). (2)

In many cases, such kernels result in solutions of Eq. (I) which exhibit self-similarity. This means that the cluster size distribution asymptotically takes the form where is a characteristic cluster size which grows in time, is a dynamical scaling exponent and denotes the scaling limit: and with fixed. Furthermore, it has been proven for the constant, sum and product kernels that such self-similar solutions are attracting for any initial data with finite moment MP2004 (). Much work in the physics literature on the theory of coagulation, following the work of Van Dongen and Ernst EVD1988 (), has focused on the problem of determining the properties of the scaling function and the exponent from the scaling properties of . An almost complete scaling theory of Eq. (I) is now known. See the review by Leyvraz LEY2003 () for a modern discussion.

A key feature of this scaling theory, and one which is of considerable importance for what follows, is the fact that the scaling properties of Eq. (I) are often very sensitive to the rate of coagulation of clusters of widely different masses. In order to parameterise this dependence in a general way, and following the notation of LEY2003 (), we introduce the scaling exponents and which specify the asymptotic behaviour of the coagulation kernel in the limit where the mass of one cluster greatly exceeds that of the other:

 K(m1,m2)∼mμ1mν2  m1≪m2. (3)

Clearly from Eq. (2), we must have . We shall use the notation to denote the -moment of the cluster size distribution:

 Mα(t)=∫∞0mαNm(t)dm. (4)

The first moment, the total mass density, , is governed by conservation of material and is thus formally conserved by Eq. (I) when . If we simply have . This accords with the intuition imparted by the mass-conserving character of the individual coagulation events. This intuition is challenged, however, by one of the more interesting phenomena to emerge from the scaling theory of the Smoluchowski equation: the so-called gelation transition. Originally conjectured by Lushnikov LUS1977 () and Ziff ZIF80 () and put on a firmer theoretical footing by Van Dongen and Ernst VDE1986 (), gelation refers to the fact that, for kernels having , mass-conservation can be spontaneously broken in finite time. That is to say, there exists a time , known as the gelation time such that

 M1(t)<∫∞0mNm(0)dm  t>t∗.

The characteristic cluster size, , typically diverges at the gelation time. The “missing” mass can be interpreted as going into a cluster of infinite mass or “gel”. The generation of arbitrarily large clusters in finite time might seem surprising but can be perfectly physical in some cases. One example is polymer gelation in which clusters merge by the formation of crosslinks and therefore do not need to move in order to coalesce. For other examples of kernels having , the solution of Eq. (I) describes only the intermediate asymptotics of the underlying physical system over some range of cluster sizes and for times less than the gelation time. Once this intermediate asymptotic range has been identified, there is no conceptual problem with the loss of mass conservation in the Smoluchowski or the generation of infinite clusters. Indeed, the gelation transition has even been observed experimentally, for example in polymer aggregation LNP1990 (), and found to exhibit dynamics in reasonable agreement with the predictions of the Smoluchowski equation.

It turns out, however, that this is not the full story. Soon after the discovery of the gelation phenomenon, it was realised by Hendriks, Ernst and Ziff HEZ1983 () that the kernel produced a series expansion for the second moment of the cluster size distribution which seemed to have zero radius of convergence in time. This led the authors to suggest the possibilty of gelation occuring at time . Detailed study of the scaling properties of the Smoluchowski equation subsequently led to formal arguments VDO1987 (), later made rigourous CC1992 (), which show that for coagulation kernels having exponent , the gelation time is actually zero. This surprising phenomenon is referred to as instantaneous gelation. Even more surprisingly, the gelation process can be complete in the sense that for . See JEO1999 () for a mathematical discussion of this point. A situation in which all mass vanishes from the system in time clearly cannot describe even the intermediate asymptotics of any physical coagulation problem. For applications such as the coagulation of polymers or colloidal aggregates, the fact that the available surface area of an aggregate cannot grow faster than its volume means that the exponent cannot be greater than 1. If the only applications of the Smoluchowski equation came from polymer science, the phenomenon of instantaneous gelation would be regarded as a mathematical pathology which need not concern physicists. Yet there are models for which it can plausibly be argued that the exponent is greater than 1. One such example is the astrophysical phenomenon of gravitational clustering SW1978 (); KON01 () which is thought to play a role in determining the large scale matter distribution of the universe. A second important example is that of differential sedimentation of water droplets falling at their terminal velocity KC2002 (); HNS2008 (), one of the processes responsible for the observed droplet size distribution in clouds FFS2002 (); FSV2006 (). Furthermore, there are even proposed heuristic solutions of Eq. (I) in the literature for such models KON01 (); HNS2008 () which seem reasonably supported by numerical simulations. The question of how this is possible, given the known mathematical results on instantaneous gelation discussed above, is the principal topic of this paper. It should be clear from the outset that the presence of an instantaneous gelation transition for a particular kernel indicates that the underlying physical model omits some process which becomes important for short timescales or for large masses. For example, in the case of droplet coalescence in clouds, the Stokes assumption for the terminal velocity of a droplet is ultimately responsible for the exponent exceeding 1. This assumption ceases to be valid for sufficiently large droplets (see KC2002 () and the references therein).

## Ii Instantaneous Gelation in the Regularised System

The gelation phenomenon in particle systems is best understood by considering a regularisation of the system and studying the behaviour as this regularisation is removed. Two natural regularisations can be found in the literature. One approach is to consider the stochastic dynamics of a finite number of particles as has been done by Lushnikov LUS2005 () for the product kernel. Another approach is to introduce a mass cut-off, , into the Smoluchowski equation. This has been done by Filbet and Laurençot FL2004 (). These are different regularisations but both show singular behaviour as the regularisation is removed when .

In this paper, we regularise by the latter method. As in FL2004 (), the cut-off is introduced in such a way that clusters having mass are removed from the system:

 ∂tNm = 12∫m1dm1K(m1,m−m1)Nm1Nm−m1 − Nm∫M−m1dm1K(m,m1)Nm1+J δ(m−1) − DM[Nm(t)]

where

 DM[Nm(t)]=Nm∫MM−mdm1K(m,m1)Nm1 (6)

describes the removal of clusters having . This regularisation explicitly breaks mass conservation. Furthermore, we shall implicitly measure all masses in terms of the monomer mass from now on so that the lower cut-off, , is set equal to 1. Clearly this is not the only way in which the problem can be regularised. In particular, by omitting from Eq. (II) one would obtain an explicitly conservative regularisation. One would expect the behaviour to be completely different in this case. The differences between conservative and non-conservative regularisations for a related aggregation–fragmentation model arising in the kinetics of waves are explored in CON2009 (). The fact that the gelation transition involves a loss of mass from the system suggests that the non-conservative regularisation is more natural. While we implement this regularisation as a hard cut-off, , we would expect the behaviour described below to be qualitatively the same for smoother regularisations provided that the non-conservative property is retained.

Figs. 1 and 2 show the mass contained in this regularised system as a function of time for a sequence of values of for the non-gelling kernel (Fig. 1) and for the gelling kernel (Fig. 2). These were obtained by numerical solution of Eq. (II) with monodisperse initial data.

All numerical results in this paper were obtained using the algorithm developed by Lee in LEE2000 (); LEE2001 (). This algorithm involves coarse-graining the masses into bins whose width increases exponentially with mass and then using the Smoluchowski equation to approximate the mass transferred per timestep as a result of the aggregation of clusters in each pair of bins. The approximation explicitly enforces mass conservation. Full details of the coarse-graining and formulae for the computation of mass transfer are provided in LEE2000 (). In a slight modification of Lee’s algorithm, the regularisation described by Eq. (6) is implemented by introducing an infinite width bin, , which accumulates any mass transferred to clusters larger than . We also modified the timestepping procedure, replacing Lee’s original explicit timestepping method with an adaptive implicit trapezoidal rule. This is essential because the Smoluchowski equation becomes increasingly stiff for larger values of to the extent that the results presented here could not be obtained with explicit timestepping. The implicit trapezoidal rule requires the solution of a set of nonlinear equations at each timestep which was done using the GSL GSL () implementation of the Rosenbrock algorithm ROS1960 (). Explicit formulae and further discussion of these modifications of Lee’s algorithm can be found in CON2009 ().

We see that for the non-gelling system, mass conservation is restored in the limit whereas for the gelling system, it is not. This latter situation will be recognisable to readers familiar with the theory of turbulence where it is widely believed that energy conservation is not restored when the limit of zero viscosity is taken in the Navier-Stokes equations FRI1995 (), a phenomenon referred to as the dissipative anomaly. A similar phenomenon is observed in the kinetics of wave turbulence CON2009 (). There is no physical contradiction in the fact that mass conservation is broken in the Smoluchowski equation in the gelling regime. It simply means that for times larger than the gelation time, the underlying conservative coagulation dynamics must be modified for the largest clusters.

Let us now consider what happens when an instantaneously gelling kernel is inserted into the regularised Smoluchowski equation. An archetypal instantaneously gelling kernel, which we shall use extensively in what follows, is the generalised sum kernel:

 Kϵ(m1,m2)=m1+ϵ1+m1+ϵ2. (7)

For this kernel and . According to the classification of Van Dongen and Ernst, it is nongelling for and instantaneously gelling for . In the marginal case, , it is the simple sum kernel which is exactly solvable, at least in the absence of a source of mononers LEY2003 (), and turns out to be non-gelling. Fig. 3 shows the sol mass in the regularised system with the generalised sum kernel with for a sequence of increasing values of the regularisation mass, . In the presence of the cut-off, the regularized gelation time, , is clearly identifiable. This regularised gelation time, although finite, decreases as increases. Extrapolating the behaviour seen in Fig. 3, it is plausible that the gelation time vanishes as consistent with the expectation that the unregularised system exhibits complete instantaneous gelation. Instantaneous gelation has not, to the best of our knowledge, been successfully numerically demonstrated in the literature previously. In the most extensive numerical study of the Smoluchowski equation to date, that of Lee LEE2001 (), the numerical difficulties posed by kernels like Eq. (7) were explored and it was concluded that there are no self-consistent solutions of Eq. (I) for such kernels. Our results demonstrate that the non-conservative regularisation, Eq. (II), provides one way around these difficulties, at least numerically, although we expect that the regularisation could also be useful for rigorous mathematical studies.

From a physical point of view, the most important observation about the results presented in Fig. 3 is that the regularised gelation time decreases extremely slowly as the cut-off is increased. decreases by a factor of less than 2 as is increased by 8 orders of magnitude. This very weak dependence means that, in practice, kernels which would exhibit instantaneous gelation, even complete instantateous gelation, for , can still have smooth, physically reasonable behaviour for finite . The regularised gelation time, , depends sufficiently weakly on the actual value of that such regularised systems may still be useful in modelling. On the basis of our numerics, we conjecture that as , the regularised gelation time decreases as for some . This would complement heuristic arguments put forward by Ben-Naim and Krapivsky BNK2003 () for gelation in finite systems of particles undergoing exchange-driven growth for which it is argued that the gelation time decreases as a power of the logarithm of the initial number of particles when the aggregation rate increases sufficiently quickly as a function of the mass of the larger cluster and thus becomes instantaneous in the limit of an infinite number of initial particles. An important question, which is left open, is to develop an analytic approach allowing the determination of the functional dependence of on and the value of if the dependence conjectured above in indeed correct. We feel it is unlikely that the numerics can be extended to sufficiently large values of to determine this dependence unambiguously although on the basis of the numerics we have available and taking into account the analytic work on the corresponding problem with a source of monomers reported below, we would not be surprised if .

We devote much of the remainder of the paper to studing what happens when a source of monomers is introduced into a system with an instantaneously gelling kernel. Such a situation has been partially analysed by Kontorovich KON01 () in the context of gravitational clustering and by Horvai et al. HNS2008 () in the context of differential sedimentation. Both studies concluded that the system should reach a stationary state at large times in which injection of monomers is balanced by the aggregation of smaller clusters into larger ones (a mass “cascade”). In both cases, the cascade is non-local in the sense that the transfer of mass from small clusters to large is dominated by the interaction of very large and very small clusters. For a detailed discussion of the criteria for locality of mass transfer in cascade solutions of the Smoluchowski equation see CRZ2004 ().

## Iii The addition model: a simplified description of runaway absorption of small clusters by large ones

Instantaneous gelation is driven by the runaway absorption of small clusters by large ones. This fact is most easily seen from the analytically tractable (but non-gelling) marginal kernel, , with source of monomers. For , the kernel in Eq.(7) is the standard sum kernel. In this case, integrating Eq.(I) with respect to (we have set ) allows us to obtain the following equation for the total number of particles in the system:

 dNdt=−2M1(t)N(t)+J. (8)

There is no gelation transition at finite time in this case so and we get a closed equation for :

 dNdt=−2JtN(t)+J. (9)

The solution is

 N(t)=12√Jπ Erfi(√Jt) e−gJt2, (10)

where is the “imaginary error function”. This formula shows an interesting feature which is surprising at first sight: while the total number of particles in the system initially increases as we add particles, it subsequently reaches a maximum and starts to decrease as shown in Fig. 4. It tends to zero () as time gets large. Meanwhile the total mass in the system increases linearly. This tells us that the long time behaviour of the system is dominated by an ever decreasing number of increasingly massive particles which immediately eat all the monomers injected into the system. Thus we already see the essential feature of non-local interactions at the level of the sum kernel. Fig. 4 demonstrates this behaviour. The numerical solution follows the analytic prediction, Eq. (10), for longer and longer times as the mass cut-off is increased. This case is marginal in the sense that there is no finite time gelation even though all the mass gets concentrated in larger and larger clusters. If , and the exponent , then big clusters are so “hungry” that they consume all smaller clusters at a rate which diverges with the cut-off, .

The dynamics of a coagulation system dominated by the absorption of small clusters by large ones is modelled in extremis by the so-called “addition model” (Brilliantov & Krapivsky BK1991 ()) in which clusters are only permitted to grow by reacting with monomers. Many of the interesting phenomena which we observe in the full coagulation problem have qualitative analogues for the addition model so we shall devote some time to studying it. The addition model corresponds to the kernel

 K(m1,m2)=12(κ(m1)δ(m2−1)+κ(m2)δ(m1−1)), (11)

where the function is typically taken to be a homogeneous function of degree . This considerably simplifies the coagulation dynamics. Although the addition model can be analyzed for a fairly general set of gelling reaction rates , we will use the following example for illustrative purposes:

 κ(m)={γm=1m1+ϵ,m=1,2,…,M, ϵ>0. (12)

Here we have changed notation slightly and written so that the instantanenous gelation regime corresponds to . The reason for introducing a separate parameter, , to describe the monomer reaction rate will become clear in Sec. VI where we show that by choosing appropriately, it is possible to make the connection between the addition model and the Smoluchowski equation more quantitative. We shall also impose the initial condition

 Nm(t)=N0δm,1, (13)

excluding the possibility of zero initial momomer concentration (this is to regularize a change of time to be defined below).

Insertion of Eqs. (11) and (12) into Eq. (I) yields the discrete set of equations

 dN1dt(t) = −2γN21−M∑m=2κ(m)N1Nm+J, (14) dNmdt(t) = κ(m−1)N1Nm−1−κ(m)N1Nm (15) 2≤m≤M,

where explicit dependence of the on has been suppressed for brevity and a cut-off, , has been introduced. Instantaneous gelation in the addition model without a cut-off for in the absence of a source was conjectured in BK1991 () and subsequently proven by Laurencot LAU1999 () providing support for our intuition that the process of aggregation of monomers by large particles is a runaway process in the full aggregation equation when . Before continuing, we would like to emphasise again that we only expect the addition model to be similar to the Smoluchowski dynamics in the instantaneously gelling regime when the dynamics is dominated by the aggregation of large clusters and monomers. It is known that the addition model has very different dynamics to the Smoluchowski equation in the non-gelling regime, . In particular it does not exhibit scaling in the absence of a source BK1991 ().

Eqs. (14, 15) simplify after a non-linear change of physical time to “monomer” time defined as follows:

 dτ(t) = N1(t)dt, (16) τ(0) = 0. (17)

Once the solution parameterized by the monomer time has been found, the inverse map to the physical time is given by

 t(τ)=∫τ0dτ′N1(τ′) (18)

This integral is proper at due to the initial condition, Eq. (13). In terms of the monomer time, Eqs. (14, 15) take the following form:

 N′1(τ) = −2γN1(τ)−M∑m=2κ(m)Nm+JN1, (19) N′m(τ) = κ(m−1)Nm−1−κ(m)Nm, (20) 2≤m≤M,

where denotes the monomer time derivative . Note that Eq. (20) is an inhomogeneous linear system of ODE’s with the inhomogeneity being a nonlinear function of the monomer concentration . It can be solved with respect to the concentrations of polymers using the Laplace transform in the monomer time. The solution can be written in the form

 Nm(τ) = ∫∞0dτ′Km(τ−τ′)N1(τ′), (21) m=1,2,…,M,

where

 Km(τ)=γκ(m)−1∫Γdω2πeωτm∏k=211+ωκ(k), (22)

the integration contour , and is a real constant: . It is easy to see from Eq. (22) that

 K(τ)=0, τ<0,

which is an expression of causality in the addition model. Substituting (21) in (19) we arrive at a nonlinear integral equation for the monomer concentration:

 N′1(τ) = −2γN1(τ)−κ(1)∫τ0dτ′QM(τ−τ′)N1(τ′) (23) +JN1(τ),

where

 QM(τ)=M∑m=2∫Γdω2πeωτm∏k=211+ωκ(k) (24)

It is clear from Eq. (23) that stays positive at all times. The transformation to monomer time is therefore well defined. To find the steady state solution of Eq. (23) we note that all the poles in the integrand of Eq. (24) have negative real parts. Thus,

 limτ→∞∫τdτ′QM(τ−τ′)N1(τ′)=N1(∞)∫∞0dτQM(τ).

Substituting this answer into Eq. (23) and setting we find that

 N1(∞) = √Jγ(2+∫∞0QM(τ)dτ) (25) = √Jγ(M+1),

where the last equality results from direct integration of the kernel, Eq. (24), over time using Cauchy’s theorem. For the specific case of the kernel, Eq. (12), the corresponding steady state answer for the polymer concentrations is

 Nm(∞) = γN1(∞)κ(m)−1 (26) = √γ√JM+1m−1−ϵ.

We acknowledge that it would be easier to find the steady state directly from Eqs. (14, 15), but Eq. (23) will be better suited to study the time evolution in what follows.

## Iv The full coagulation problem with ν>1 in the presence of a source of monomers

In the presence of a source of monomers, Eq. (I) has a formal stationary solution which scales for large masses as which describes a cascade of mass from small masses to large CRZ2004 (). This solution is only valid if the collision integral is convergent. For the case of the kernel Eq. (7), this convergence criterion fails for CRZ2004 () presumably reflecting the tendency for the system to gel instantaneously in this regime. When a cut-off is introduced, a stationary state may be reached if a source of monomers is present although this stationary state must involve the cut-off rather than being of the cascade type. In HNS2008 () such a stationary state was found for the case of differential sedimentation () which scaled for large masses as .

In this section we perform a systematic analysis of the stationary state for the model kernel

 K(m1,m2)=max(m1,m2)ν, (27)

which captures the essential features of instantaneously gelling systems while permitting some convenient simplifications of the equations. If we assume that the system reaches a stationary state and define

 ΓM=M∑m=1mνNm, (28)

then some manipulations of the discrete analogue of Eq. (I) show that the stationary state can be expressed via a recursion relation:

 Nm=12∑m−1p=1max(p,m−p)νNpNm−pΓM+∑m−1p=1(mν−pν)Np. (29)

which can be iterated to find once one observes that the monomer concentration is fixed in the stationary state as:

 N1=JΓM (30)

For each value of , this iteration procedure can, in principle, be used to self-consistently determine . In this paper, we adopt a different approach which is based on the assumption that the mass transfer is nonlocal in mass space. Following the approach of HNS2008 (), Eq. (I) can be approximated by:

 ∂Nm∂t=−M1∂∂m(mνNm)−MνNm (31)

where the moments and are now to be computed with the cut-off retained:

 Mα=∫M1mαNmdm. (32)

The source and sink terms have been omitted for now. In obtaining Eq. (31) it has been assumed that the integrals and are dominated by their lower and upper limits of integration respectively. The consistency of this assumption will be determined a-posteriori. Since and do not depend on , Eq. (31) has the stationary solution:

 Nm=Cexp[βm1−νν−1]m−ν, (33)

where is an arbitrary constant and we have, for convenience, introduced the parameter,

 β=MνM1. (34)

Although we do not know the value of a-priori, it can now be determined self-consistently due to the fact that, when is given by Eq.(33), the integrals and can be expressed explicitly in terms of incomplete gamma functions:

 M1(β) = Cν−1(1−νβ)ν−2ν−1[Γ(ν−2ν−1,−βM1−νν−1)−Γ(ν−2ν−1,−βν−1)] (35) Mν(β) = Cν−1(1−νβ)11−ν[Γ(11−ν,−βM1−νν−1)−Γ(11−ν,−βν−1)]. (36)

After some algebra, Eq. (34) reduces to the following consistency condition for the value of :

 Γ(ν−2ν−1,−βM1−νν−1)−Γ(ν−2ν−1,−βν−1)=−1ν−1Γ(11−ν,−βM1−νν−1)+1ν−1Γ(11−ν,−βν−1) (37)

Given , this can be solved numerically for any value of . Such a numerical investigation indicates that is a slowly increasing function of . Furthermore, since we expect the integrals and to be dominated by their respective lower and upper limits of integration as , we would expect Eq. (37) to satisfy the following asymptotic balance for large :

 −Γ(ν−2ν−1,−βν−1)∼−1ν−1Γ(11−ν,−βM1−νν−1). (38)

Since numerics indicate that increases (but only slowly) as grows and , the argument of the left-hand gamma function, , should increase in the limit of interest. Likewise, the argument of the right-hand gamma function, , should decrease in the limit of interest. The relevant leading order asymptotics are then:

 Γ(a,z) ∼za−1e−z as z→∞ (39) Γ(a,z) ∼−1aza as z→0. (40)

Substituting these into Eq. (38), one finds that the leading order terms balance as provided

 β∼log(Mν−1). (41)

Whilst the index could be absorbed within the implied front factor in , we keep it in the following expressions since it captures how the dependence on the cut-off, , goes away when approaches unity. Knowing , Eqs. (40) and (39) can now be used to detemine the asymptotic behaviour of and from Eqs. (35) and (36). After some work one finds:

 M1∼CMlog(Mν−1) as M→∞ (42) Mν∼CM as M→∞. (43)

It now remains to determine the constant in Eq. (33) in order to close the argument and allow us to check for consistency. This can be done using global mass balance. Multiplying Eq. (II) (with the kernel Eq. (27)) by and integrating from 1 to yields the global mass balance condition:

 J=∫M1dm mNm∫MM−mdm1 mν1Nm1. (44)

Our assumption of nonlocality allows to replace the inner integral by in the limit of large since it is dominated by its upper limit. Consequently the global mass balance condition in the limit of large is

 J=M1Mν=C2M2log(Mν−1), (45)

where we have used Eqs. (42) and (43). This gives

 C=√Jlog(Mν−1)M. (46)

Putting this all together with Eq. (33) the asymptotic solution of Eq. (31) is

 Nm=√Jlog(Mν−1)MMm1−νm−ν. (47)

Everything becomes self-consistent: becomes independent of and diverges as when gets large thereby justifying the use of Eq. (31). This theory does a good job of capturing the general features of the stationary state: stretched exponential decay at small cluster sizes followed by a cross-over to a power law decay with exponent for large cluster sizes but with an amplitude which decreases with the cut-off, . Fig. 5 compares Eq. (47) with the stationary state obtained from numerical simulations of Eq. (II) with for a range of values of . The agreement is excellent given that there are no adjustable parameters in Eq. (47). Note that the large mass behaviour, , is in agreement with the stationary solution of the addition model, Eq. (26), obtained in Sec. III.

Finally, we can now quantify the rate at which the system becomes singular as the cut-off, , is increased. Using Eq. (47) to calculate the total particle density in the stationary state we get

 N=√J(M−MM1−ν)M√log(Mν−1)∼√Jlog(Mν−1)  as M→∞. (48)

The density of particles in the stationary state thus vanishes as the cut-off is removed which is the signature of instantaneous gelation. The approach to zero is logarithmically slow however. The asymptotic estimate, Eq. (48), is compared against the total density measured from a sequence of numerical simulations with and different values of the cut-off in the inset of Fig. 5. The theory again produces convincing agreement with numerics without any adjustable parameters.

## V Approach to the Stationary State

Numerical simulations indicate that the approach to the stationary state is interesting. Fig. 6 shows a plot of the total density in the system as a function of time for the generalised sum kernel, Eq. (7), for a range of values of with cut-off . In the nonlocal regime, , the approach to the stationary state is characterised by long-lived transient oscillations of the density in the system. This is surprising since there is nothing of an oscillatory or cyclic character in the underlying coagulation dynamics. We note, however, that long-lived oscillations of the total droplet density have been observed in experiments on the phase separation of binary mixtures with nucleation and coagulation of droplets VAV2007 (); BV2010 () which may be related to the effect noted here. From Fig. 6, it is not obvious that these oscillations are decaying at all for . From our numerical investigations we believe that these oscillations are indeed transient although the rate of decay becomes very slow as the cut-off is increased. Although the oscillations are quite nonlinear in character, we defined a quantity analogous to the Q-factor of a linear oscillator by measuring the ratio of subsequent maxima of the signal and defining the Q-factor to be the estimated numerical limit of this ratio as (for all of our numerical experiments with , we found that this ratio does indeed seem to converge to a finite value). Thus a Q-factor of 1 would correspond to persistent oscillations.

A summary of these numerical experiments is presented in Fig. 7 which shows numerically estimated Q-factors for the oscillations observed for a range of values of and . The most evident trend from the data is the fact that the oscillations come closer and closer to a Q-factor of 1 as is increased. These measurements indicate that the oscillations are decaying in all cases but only very slowly when the cut-off becomes large. The behaviour as a function of is less clear with some possible evidence for a weak maximum.

We offer the following heuristic explanation for how oscillations can be generated in this system. For , the case for which instantaneous gelation should occur in the absence of a cut-off, the indication of Fig. 3 is that the particle density drops close to zero in finite time (but this time is not zero). When monomers are injected into the system, until the cut-off is felt, the mass in the system just increases linearly. However, as soon as large particles are generated, their absorption of the monomers is a runaway process which very rapidly converts the monomers which have accumulated in the system into large particles which are immediately removed by the cut-off. This then resets the system close to its initial state in which there are almost no particles in the system. The dynamics then repeats.

## Vi Approach to the stationary state in the addition model

We do not, as yet, have an analytic characterisation of the oscillations described in the previous section. Instead, we conclude this article with an analysis of the analogous problem for the simpler case of the addition model discussed in Sec. III for which the oscillatory approach to the stationary state, Eq. (26), can be derived.

The addition model aims to model the most important contribution to the dynamics of the full aggregation model with non-local kernels by concentrating on interactions between “light” and “heavy” particles only. At the level of modeling, there is no reason to expect that the parameters of the addition model are the same as the analogous parameters in the full model. Before proceeding with calculations, we first present a semi-heuristic argument for the most appropriate choice of parameters. Let us first clarify what is meant by “light” and “heavy” in terms of the steady state of the full model, Eq. (47). Heavy particles have masses in the polynomial tail of Eq. (47), . Light particles belong to the stretched exponential region of Eq. (47), . Therefore, neglecting logarithmic corrections, we can treat all light particles as effective monomers. Hence, the parameters and for in the addition model are the same as in the full model, but the parameter should be chosen as an effective interaction rate of all light particles in the full model with themselves. The easiest way to choose this rate is to require that the steady state of the addition model gives the correct description of the mass distribution of the full model at large masses. Comparing Eq. (26) with the large- behaviour of Eq. (47) we see that we should choose:

 γ∼1M. (49)

For large times let us look for a solution of Eq. (23) in the form

 N1(τ)=√Jγ(M+1)(1+Aeζτ), (50)

where . The limit of large times formally corresponds to the limit of small amplitude . The requirement that Eq. (50) solves Eq. (23) up to terms of order leads to the following non-linear eigenvalue problem:

 ζ+γ(M+3)+γM∑m=2e−∑mk=2log(1+ζκ(k))=0. (51)

Notice that for gelling kernels,

 limm→∞m∑k=2log(1+ζκ(k))<∞

Therefore, the sum in the left hand side of Eq. (51) diverges as the first power of :

 limM→∞1MM∑m=2e−∑mk=2log(1+ζκ(k))=e−∑∞k=2log(1+ζκ(k))

As we are interested in the limit of large cut-off mass , let us solve the eigenvalue problem, Eq.  (51), in the vicinity of . Bearing in mind that , the leading term of the large -expansion of the left hand side of Eq. (51) gives

 e−log(1+ζ)−∑∞k=2log(1+ζκ(k))=−1, (52)

or equivalently,

 log(1+ζ)+∞∑k=2log(1+ζκ(k))=πi(1+2q) (53)

where . We need to find the solution to the above equation with the largest real part. It turns out that an analytical solution is possible for the set of reaction rates Eq. (12) in the limit Let us Taylor expand the logarithms in Eq. (53) and switch the orders of summation in anticipation of seeking an asymptotic expansion for . We obtain:

 ∞∑n=1(−1)n+1Snζn=πi(1+2q)

where

 Sn=1n(1+∞∑k=2κ(k)−n).

For the reaction rates given by Eq. (12) with , we have the following behaviour for the :

 S1= ∼ 1ϵ, Sn= ∼ 1, n=2,3,….

The fact that for small values of allows one to find the solution to Eq. (53) with the largest real part ( or ) as an asymptotic expansion in :

 ζ = ±i(π1S1+π3S3S41−2π3S22S51)+(−π2S2S31+π4S4S51) (54) +O(ϵ6).

We observe that and . Therefore the monomer concentration slowly decays to the stationary value in an oscillatory manner with many periods of oscillations per inverse decay rate. Accordingly the oscillations’ quality factor is close to for :

 Q=e−2π2S2S21+O(ϵ3). (55)

The independence of the quality factor of the cut-off mass (up to possible logarithmic corrections) is confirmed. It is also straightforward to check that the physical period of oscillations is -independent: notice that . Therefore, at large times the transformation Eq. (16) reads:

 τ(t)=τ0+N1(∞)t+O(eReζt).

Substituting this expression into Eq. (50) we find that in the limit ,

 N1(τ)=√J(1+Aeζ√Jt), (56)

where is given by Eq. (54). Therefore, the period of oscillations does not depend on the cut-off. It scales as the inverse square root of the flux of monomers:

 T=2√JS1(1+O(ϵ3)). (57)

The period diverges in the limit , which implies that oscillations disappear as we increase the locality of the aggregation kernel. These conclusions are in agreement with our numerical measurements on the full coagulation equation presented in Sec. V.

## Vii Conclusions

To conclude, the fact that certain aspects of the regularized Smoluchowski equation, such as the total density, depend so weakly on the value of the cut-off, , used to regularise the system means that instantaneously gelling kernels are still potentially reasonable as physical models. Their use to describe gravitational clustering or differential sedimentation is not in contradiction with mathematical results which indicate that the density vanishes for all positive time in the unregularised system provided one is willing to accept non-universal dependencies on a cut-off as part of the model. We have presented an essentially complete analysis of the stationary state of the regularised system in the presence of a source of monomers when the regularisation is done be removing all clusters having size larger than the cut-off. In this case, a power law scaling is observed for large cluster sizes whose exponent is universal (depending only on the scaling properties of the kernel) but whose prefactor is strongly cut-off dependent. The scale-to-scale mass balance is of a nonlocal character rather than being described by a mass “cascade” as is the case for regular gelling and non-gelling kernels. Finally, we found that this stationary state is approached in an interesting way with long-lived oscillatory transient resulting from interaction between the cut-off and the source. A quantitative analysis of the approach to the steady state was provided in the simplified case of the addition model which expresses the most extreme case of nonlocal mass transfer in which clusters can only grow by interaction with monomers.

## Acknowledgements

C.C. would like to acknowledge useful conversations with E. Ben-Naim, P. Krapivsky and S. Nazarenko.

## References

• (1) M. Ernst, in Fractals in Physics, edited by L. Pietronero and E. Tosatti (North Holland, Amsterdam, 1986) p. 289
• (2) M. V. Smoluchowski, Z. Phys. Chem. 91, 129 (1917)
• (3) D. J. Aldous, Bernoulli 5, 3 (Feb. 1999)
• (4) C. Connaughton, R. Rajesh, and O. Zaboronski, Physica D 222, 97 (2006)
• (5) G. Menon and R. L. Pego, Comm. Pure Appl. Math. 57, 1197 (2004)
• (6) M. H. Ernst and P. G. J. van Dongen, J.Stat. Phys. 1/2, 295 (1988)
• (7) F. Leyvraz, Phys. Reports 383, 95 (Aug. 2003)
• (8) A. A. Lushnikov, Dokl. Akad. Nauk SSSR 236, 673 (1977)
• (9) R. Ziff, J. Stat. Phys. 23, 241 (1980)
• (10) P. G. J. van Dongen and M. H. Ernst, J. Stat. Phys. 44, 785 (1986)
• (11) A. Lushnikov, A. Negin, and A. Pakhomov, Chem. Phys. Lett. 175, 138 (1990)
• (12) E. M. Hendriks, M. H. Ernst, and R. M. Ziff, J. Stat. Phys. 31, 519 (1983)
• (13) P. van Dongen, J. Phys. A: Math. Gen. 20, 1889 (1987)
• (14) J. Carr and F. P. da Costa, Zeit. Math, Phys. 43, 974 (1992)
• (15) I. Jeon, J. Stat. Phys. 96, 1049 (1999)
• (16) J. Silk and S. D. White, Astrophys. J. 223, L59 (1978)
• (17) V. Kontorovich, Physica D 152–153, 676 (2001)
• (18) V. I. Khvorostyanov and J. A. Curry, J. Atm. Sci. 59, 1872 (2002)
• (19) P. Horvai, S. V. Nazarenko, and T. H. M. Stein, J. Stat. Phys. 130, 1177 (2008)
• (20) G. Falkovich, A. Fouxon, and M. G. Stepanov, Nature 419, 151 (2002)
• (21) G. Falkovich, M. G. Stepanov, and M. Vucelja, J. Appl. Met. Clim. 45, 591 (2006)
• (22) A. A. Lushnikov, Phys. Rev. E 71, 046129 (2005)
• (23) F. Filbet and P. Laurençot, Archiv der Mathematik 83, 558 (2004)
• (24) C. Connaughton, Physica D 238, 2282 (2009)
• (25) M. Lee, Icarus 143, 74 (2000)
• (26) M. Lee, J. Phys. A: Math. Gen. 34, 10219 (2001)
• (27) M. Galassi, J. Davies, J. Theiler, B. Gough, G. Jungman, P. Alken, M. Booth, and F. Rossi, GNU Scientific Library Reference Manual - 3rd Edition (Network Theory Ltd., 2009) ISBN 0-9546120-7-8
• (28) H. H. Rosenbrock, Comp. J. 3, 175 (1960)
• (29) U. Frisch, Turbulence: the legacy of A. N. Kolmogorov (Cambridge University Press, Cambridge, 1995)
• (30) E. Ben-Naim and P. Krapivsky, Phys. Rev. E 68, 031104 (2003)
• (31) C. Connaughton, R. Rajesh, and O. Zaboronski, Phys. Rev. E 69, 061114 (2004)
• (32) N. V. Brilliantov and P. L. Krapivsky, J. Phys. A Math. Gen. 24, 4789 (Oct. 1991)
• (33) P. Laurençot, Nonlinearity 12, 229 (Mar. 1999)
• (34) J. Vollmer, G. K. Auernhammer, and D. Vollmer, Phys. Rev. Lett. 98, 115701 (2007)
• (35) I. J. Benczik and J. Vollmer, Europhys. Lett. 91, 36003 (2010)
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters