Anomalous dimensions of the Smoluchowski coagulation equation

The coagulation (or aggregation) equation was introduced by Smoluchowski in 1916 to describe the clumping together of colloidal particles through diffusion, but has been used in many different contexts as diverse as physical chemistry, chemical engineering, atmospheric physics, planetary science, and economics. The effectiveness of clumping is described by a kernel $K(x,y)$, which depends on the sizes of the colliding particles $x,y$. We consider kernels $K = (xy)^{\gamma}$, but any homogeneous function can be treated using our methods. For sufficiently effective clumping $1 \ge \gamma>1/2$, the coagulation equation produces an infinitely large cluster in finite time (a process known as the gel transition). Using a combination of analytical methods and numerics, we calculate the anomalous scaling dimensions of the main cluster growth, calling into question results much used in the literature. Apart from the solution branch which originates from the exactly solvable case $\gamma = 1$, we find a new branch of solutions near $\gamma = 1/2$, which violates scaling relations widely believed to hold universal.

The coagulation (or aggregation) equation was introduced by Smoluchowski in 1916 to describe the clumping together of colloidal particles through diffusion, but has been used in many different contexts as diverse as physical chemistry, chemical engineering, atmospheric physics, planetary science, and economics.The effectiveness of clumping is described by a kernel K(x, y), which depends on the sizes of the colliding particles x, y.We consider kernels K = (xy) γ , but any homogeneous function can be treated using our methods.For sufficiently effective clumping 1 ≥ γ > 1/2, the coagulation equation produces an infinitely large cluster in finite time (a process known as the gel transition).Using a combination of analytical methods and numerics, we calculate the anomalous scaling dimensions of the main cluster growth, calling into question results much used in the literature.Apart from the solution branch which originates from the exactly solvable case γ = 1, we find a new branch of solutions near γ = 1/2, which violates scaling relations widely believed to hold universal.
Smoluchowski's equation [1] has its origin in physical chemistry, but more generally furnishes a fundamental description of the formation of larger objects by the aggregation of smaller entities.It appears in many physical problems such as planetesimal accumulation, mergers in dense clusters of stars, aerosol coalescence in atmospheric physics, and polymerization and gelation (see [2][3][4][5][6][7]), but also in chemical engineering [7], and the social sciences [8,9] It describes the evolution of the density c(x, t) of particles of size x at time t, taking into account the formation of new clusters of size x by the aggregation of pairs of size x − y and y respectively, as well as the disappearance of clusters of size x forming a larger one: Here the function K(x, y) (known as the coagulation kernel) describes the probability for two particles of sizes x and y to stick together.The behavior of solutions to (1) depends crucially on the degree of homogeneity of K. To explore this, here we restrict ourselves to the class of models described by K(x, y) = (xy) γ , for which the degree is 2γ.This kernel applies to branched polymers with surface interactions [4,10], and to fractal clusters more generally [11,12], but stands for a much broader class of models whose asymptotic behavior for large cluster sizes scales with an exponent 2γ.One of the fundamental problems in the field is to relate, by solving (1), γ to the scaling exponent β determining the typical size of clusters, and the gel exponent τ , giving the power-law size distribution of clusters [4].Thus by measuring β or τ , one is then able to infer fundamental mechanisms of aggregation, in phenomena as diverse as planatesimal formation [7], aerosol dynamics [13], or pipeline fouling caused by asphaltenes [14].
Only for γ = 1 can (1) be solved explicitly [15][16][17][18], for more general kernels studies have relied on discrete particle simulations and ad-hoc scaling arguments (see [19], [20] and references therein).It is therefore of enormous importance to develop mathematical methods able to provide novel information on the behavior of solutions to (1).
For 1/2 < γ ≤ 1, (1) develops singularities in finite time, such that, starting from an initial particle size distribution c(x, 0) with all its moments M i = ∞ 0 x i c(x, 0)dx bounded, there is a certain time t 0 such that all moments M i for i ≥ i 0 diverge (see [21] and references therein).This phenomenon, which has the character of a phase transition [4], is called finite time gelation (at a gelation time t 0 ), and indicates the aggregation of particles in a single cluster of infinite mass.In practice, of course, the singularity will be cut off by the finite size of the total number of particles available.On the other hand, if γ ≤ 1/2, solutions exist globally in time [22].
As in many other physical problems involving diverging quantities (cf.[23]), we assume that the approach to the singularity is selfsimilar, of the form c(x, t) = t α ψ(xt β ); here t = t 0 − t is the time distance to the singularity.However, this selfsimilar structure has so far only been established for γ = 1, while for γ < 1 selfsimilar solutions have not been determined explicitly.Discrete numerical simulations (cf.[24], [20]) appear to show selfsimilarity of the second kind [25], for which similarity exponents cannot be determined from dimensional considerations, or from symmetry arguments.This fact was proven in [26] for 1 − γ sufficiently small, without calculating the exponents explicitly.
Inserting the similarity form into (1) and balancing powers of t , one finds α − 1 = 2α − (2γ + 1)β, so that which will form the basis of our analysis; s(t) = constt −β is known as the typical cluster size [19].In addition, one obtains an integral equation for the similarity profile ψ(ξ).The scaling form ( 2) is also known as the "selfpreservation hypothesis" [3]; for example, using a rescaling analogous to (2), in Fig. 7.11 of [3] the distribution of aerosol particles, taken from experiment, is collapsed onto a single profile ψ(ξ).
Imposing that the mass M 1 = ∞ 0 xc(x, t)dx is conserved by selfsimilar solutions, one would obtain β = 1/(2γ − 1).However, this need not be the case; rather, mass only has to be conserved by the full solution to (1), and not necessarily by an asymptotic solution of the form (2). Below we will consider such asymptotic solutions which do not conserve mass and for which the cluster size diverges, and hence β cannot be deduced from mass conservation.
By considering the Laplace transform of ψ [26], defined as the equation for Φ becomes an equation similar to those of non-local transport equations treated by us [27], and for which we have developed efficient numerical treatments.The behavior of Φ(η) for large arguments represents the distribution of small clusters; if ψ(ξ) ∝ ξ −τ [19], then Φ(η) ∝ η τ −2 , where τ is known as the (pre)gel exponent [28].This implies that as t 0 is approached, the cluster size distribution function approaches a power law c(x, t 0 ) ∝ x −τ , which will once more link to the kernel's exponent γ.Similar power law distributions occur in many other fields, such as Zipf's law [8], and have been proposed to understand phenomena such as bank mergers [9].
To derive the integral equation for Φ(η) to be used in the following, we start from (2.6) of [26], which, adopting the notation of the present paper, reads where ν = (1 − (2γ − 1)β)/β.Defining and inserting (3) into (5), one obtains, interchanging the order of integration, and performing the integral over ζ, that F equals the square bracket in (4).This shows that with F defined by (5).
The left hand side of ( 6) corresponds to the time derivative of the Smoluchowski equation, which is expected to vanish for large η in order for the solution to match to the "background" of the distribution of small clusters.This matching condition [29] then implies that Φ(η) ∝ η ν for large η, which leads to the scaling relation τ = 2+ν = 1+2γ −σ [19], which relates the gel exponent τ with the exponent σ ≡ β −1 , determining the typical cluster size.We will see below that this is true only for the "lower" branch of solutions, which grows out of the classical case γ = 1, but fails for the "upper" branch, first reported here.
We are looking for solutions to (6) with Φ(η) regular at the origin and Φ(η) ≈ Aη ν for η → ∞, where A is a constant to be found as part of the solution, along with ν.In the exactly solvable case γ = 1, we have F = Φ(η), so (6) becomes νΦ + ηΦ = ΦΦ , the same as for the kinematic wave equation [23,26].This equation has an infinite sequence of regular solutions η = (1 + j)Φ/j + BΦ 1+j , where j = 1, 2, . . ., and B is an arbitrary constant.For the "ground" state j = 1, the cluster size exponent is β = 2, and ν = 1/2.A sequence of non-trivial solution branches, exhibiting anomalous scaling exponents, emanate from each of these exact solutions; we will focus on the ground state branch, which is expected to be attracting, while all other branches are unstable.

pansion Φ ≈ Aη
for large η.As described in more detail in the supplementary material, ( 6), (7) are solved using a Newton method, continuing from the ground state solution at γ = 1.
The resulting values of ν, which make up the "lower branch", are shown in Fig. 1 as the lower solid line, emanating from ν = 1/2 for γ = 1.In a number of widely cited papers [10,19,28,30,31], it was proposed on the basis of ad-hoc conditions on the behavior of the similarity equation for small clusters, that σ = (2γ − 1)/2 and ν = γ − 1/2.The latter is shown as the dotted line in Fig. 1, clearly in strong disagreement with the actual solution of ( 6), as anticipated in [26].Indeed, in (6) the behaviors for small and for large clusters are in fact coupled, which leads to anomalous scaling exponents [25,29], invalidating a simple linear scaling.Our results also agree very well with the numerics of [20] (circles), obtained using a particle-based description.For convenience, in the supplementary material we also give an interpolation formula, which describes the solution branch to three decimal places.
We have also calculated solution branches which emanate from the higher-order solutions at γ = 1, which are known to be unstable [29].It is therefore likely that the entire higher-order branches are unstable.Indeed, we have also solved the time-dependent evolution equations in Laplace space for a particular value γ = 0.7942, and found the solution to converge onto the stable ground state solution, shown in Fig. 1.
As γ decreases toward 1/2, the correction exponent 1− 2γ +2ν in (7) becomes ever closer to ν, so a larger domain is needed to correctly describe the asymptotics for large η, as seen in Fig. 2. The two exponents become identical for ν = 2γ − 1 ≡ δ, which suggests the appearance of a new branch of solutions, for which Φ(η) ∼ η δ , and which we call the "upper branch", also shown in Fig. 1.
To understand the transition between the two branches, we write a formal solution of the similarity equation ( 6): If the integral in ( 8) is convergent, then the second term scales like η ν , and the first term in ( 8) is seen to be subdominant.We anticipate a non-uniform convergence of the lower branch toward the upper branch as γ → γ c .Let us assume that as suggested by Fig. 2, for 1 , where the exponent is positive on the lower branch, so that as η c → ∞, the prefactor A diverges.Indeed, a best fit to the numerical data yields A ≈ A 0 /(γ − 0.521) α , with α = 0.348 and A 0 = 0.44, a blowup at γ very close to the extrapolated value of γ = γ c = 0.5214.
If on the other hand the integral in ( 8) is divergent, both terms scale in the same way, and using the known asymptotic behavior of F , balancing both sides yields to leading order as η → ∞.
To find what we call the upper branch, we enforce (10) instead of (7) for large η; to account for the scale invariance of Φ, we impose Φ (0) = 1.A typical profile on the upper branch is shown in Fig. 3, using a logarithmic scale, except near the origin.The dashed line is the expected asymptotics (10), and γ is the same as in Fig. 2, showing the lower branch.This demonstrates that for a range of γ values above γ c , there are multiple solutions.This is also clear from the phase diagram in Fig. 1, where both branches are shown.The lower branch ends at γ c , where it meets the upper branch, as seen in the inset.We show in the supplementary material that for γ = 1/2, the scaling function Φ(η) behaves asymptotically like a logarithm.Thus for γ 1/2 one needs an ever larger computational domain to describe the crossover between logarithm and power law, and we are not able to continue the upper branch all the way to γ = 1/2.A more fundamental concern is that on the upper branch, η ν ≡ η α/β does not equal the asymptotic behavior Φ ∝ η δ , represented by circles in Fig. 3.As a result, the Laplace transform of the cluster size distribution behaves for large arguments like ω(µ, t) ≈ β Āµ δ t −1 , which for t → 0 diverges, and hence does not match the expected static distribution at large arguments.However, the outer solution ω(µ, t) = t λ ϕ(ξ), ξ = µt −(1+λ)/δ , (11) succeeds in bridging this time dependence with a static outer distribution; λ is another anomalous exponent to be determined.The similarity equation becomes −λϕ+(1+λ)ξϕ ξ /δ = F F , which for ξ → 0 has solutions of the form ϕ ≈ β Āξ δ + Gξ α .Linearizing around the small perturbation Gξ α , one finds For small δ, λ, this simplifies to whose dominant (smallest) solution is α = α 0 = 1.30737 . . . .The exponent λ plays the role of a nonlinear eigenvalue, which has to be found such that the asymptotics of ϕ at infinity are satisfied.On the other hand, for large ξ we have to require that ϕ(ξ) → Φ 0 ξ λδ/(1+λ) .This will ensure that for small cluster sizes, ω matches onto a static cluster size distribution.There are many directions in which to extend the present research.First, it would be interesting to consider (1) for times after the singularity, and establish postgel solutions, including the scaling relations they satisfy.Second, we have not yet explored the range of exponents 0 ≤ γ ≤ 1/2, for which an infinitely large cluster is formed in the limit t → ∞ only.For example, it is not known whether in this case anomalous dimensions once more appear, or whether the explicit expressions for σ and τ in terms of γ, established in the literature [19], hold.
Third, many studies have looked at other types of homogeneous kernels [4], such as K = x µ y ν or K = x λ +y λ .For example, an important question is whether exponents only depend on the degree of homogeneity λ = µ + ν, or whether they are more sensitive to the structure of the kernel.Indeed, our methodology extends to much more general kernels, for example of the form K(x, y) = (x µ y ν + x ν y µ ) /2 (with µ + ν = 2γ).This would lead to the integro-differential equation ( 6), with the right hand side replaced by 1 2 , where F µ and F ν are defined as in (5), with γ replaced by µ and ν, respectively.In wave turbulence [32][33][34][35], similar integral equations arise, which have not been solved explicitly, as we do here.Instead, the theory rests on scaling assumptions similar to those which in coagulation theory were found to be invalid.Here a stationary turbulent spectrum would correspond to postgel solutions, which evolve out of the initial singularity [35].It is therefore possible that a more careful treatment of the integral equations of wave turbulence yields anomalous dimensions, as has been conjectured [34], which would change the scaling exponents of (say) the velocity field of the turbulence.
In conclusion, integro-differential equations represent an area where some of today's most challenging unsolved problems in statistical mechanics and in fluid dynamics come together.Using self-similar solutions to the Smoluchowski equation, we showed that such integral equations have many unexpected properties, which challenge longheld beliefs.

FIG. 1 .
FIG.1.The exponent ν for γ between 1/2 and 1.One branch (the lower branch) starts from γ = 1, ν = 1/2 on the right, and ends by intersecting another branch (the upper branch) at γc ≈ 0.5214.The scaling ν = γ − 1/2[19] is shown as the dotted line.The inset shows a detail of the bifurcation, where the two branches meet at γc.The circles mark the particle-based simulations of Lee[19].

FIG. 2 .
FIG.2.The local exponent of Φ(η) over a large domain,(7) being used for η > η k = 10 70 , for γ = 0.521829.The solid line is the full solution on the lower branch, the dashed line is(7).The solid horizontal line on the left is δ = 2γ − 1 = 0.04366, the dashed horizontal line on the right is ν = 0.02808.