Gelation in coagulating systems

Physica D 222 (2006) 37–53 www.elsevier.com/locate/physd

Review

Gelation in coagulating systems A.A. Lushnikov ∗ Department of Physical Sciences, University of Helsinki, P.O. Box 64, FIN-00014, Helsinki, Finland Available online 6 September 2006

Abstract A review of recent achievements in the theory of sol–gel transitions is given. This paper outlines the basic ideas of a stochastic approach and describes its possible application as a basis for a theory of gelation in coagulating systems. The kinetics of coagulation is described in terms of the probability W to find in the system a given set of monomers, dimers, etc., at time t. This approach allows for a consideration of finite coagulating systems, i.e., those containing a finite number of monomeric units. The evolution equation for the generating functional of W is formulated and solved exactly for the coagulation kernel proportional to the product of masses of coagulating particles for which the traditional Smoluchowski’s scheme leads to a paradox related to the loss of the total mass of the coagulating particles. This paradox is resolved within the scheme described in this paper. It is shown that a single giant object (a gel) invisible in the thermodynamic limit arises after a critical time and begins to consume the mass of smaller but higher populated fraction (the sol). This approach is then applied for considering the time evolution of random graphs. The distribution of the linked components of a random graph over their orders and numbers of edges is found exactly. The paper is concluded with a short consideration of the sol–gel transition in systems with instantaneous sinks at large particle mass. It is demonstrated how the sol–gel transition manifests itself in such systems. c 2006 Elsevier B.V. All rights reserved.

Keywords: Sol–gel transition; Finite coagulating systems; Exactly solvable model

1. Introduction Although the phenomena related to the appearance of new collective degrees of freedom are quite common for theoretical physics, every time they appear it impels theorists to invent something special for their consideration. It is enough to mention phase transitions, turbulence, superconductivity and superfluidity, formation of surfaces, rotation of atomic nuclei, and many other interesting manifestations of collective effects whose considerations have made the theorists rack their brains in attempts to find suitable routes for treating these phenomena. The sol–gel transition in coagulating systems is also not an exception. At first sight the coagulation process looks rather offenceless: a system of M monomeric objects begins to evolve by pair coalescence of g- and l-mers according to the scheme, (g) + (l) −→ (g + l). ∗ Fax: +358 9 191 3860.

E-mail address: [email protected]. c 2006 Elsevier B.V. All rights reserved. 0167-2789/$ - see front matter doi:10.1016/j.physd.2006.08.002

(1)

And there is not a problem to write down the kinetic equation governing the process; everyone can do it: g−1

∞ X dcg 1X = K (g − l, l)cg−l cl − cg K (g, l)cl . dt 2 l=1 l=1

(2)

This is the famous Smoluchowski equation. Here the coagulation kernel K (g, l) is the transition rate for the process given by Eq. (1). The first term on the right-hand side (RHS) of Eq. (2) describes the gain in the g-mer concentration cg (t) due to coalescence of (g − l)- and l-mers, while the second one is responsible for the losses of g-mers due to their sticking to all other particles. In what follows we will use the dimensionless version of this equation, i.e., all concentrations are measured in units of the initial monomer concentration c1 (0) and the time in units of 1/c1 (0)K (1, 1). More details can be found in review articles [1–3]. Half a century ago it had become clear that there is something wrong with Eq. (2). An attempt [4] to find an exact solution to this equation for the kernel proportional to the

38

A.A. Lushnikov / Physica D 222 (2006) 37–53

and the particle concentrations cg (t) = n¯ g (t)/V.

Fig. 1. The total number and total mass concentrations of sol particles are shown as functions of time (dimensionless units). After the critical time t = tc the mass concentration ceases to conserve, for a massive gel particle forms and begins to consume the mass of the sol. On the other hand, the number concentration does not feel the loss of one (although very big) gel particle. Still the post-critical behavior of the number concentration found from Eq. (77) differs from that predicted by the Smoluchowski equation (n(t) = 1 − t, dash line). From Ref. [19].

masses of coalescing particles K (g, l) ∝ gl

(3)

had led to a strange conclusion that the total mass concentration ceases to conserve after a finite time t = tc (in what follows tc is referred to as the critical time)P and the second moment of the particle mass distribution φ2 = g 2 cg has a singularity, φ2 (t) ∝

1 . tc − t

(4)

Even more strange is the fact that at t = tc nothing wrong happens either to the particle mass spectrum or the particle number concentration. The whole situation is displayed in Fig. 1. Immediately the problems of existence of the solution to Eq. (2) and of its uniqueness were posed and resolved [5–8]. But this did not help to answer the question, what does happen after t = tc ? On the other hand, it is clear that if we consider a finite coagulating system, then at any time we see a number of bigger and bigger particles whose total mass M does not disappear somewhere. Such a system can be characterized by the set {n g } of occupation numbers of g-mers. We can also introduce the probability W ({n g }, t) for the realization of a given set at time t. Now the evolution of the coagulating system can be fully described in terms of W . But a truncated description is also admissible: we can introduce the average occupation numbers X n¯ g (t) = W ({n g }, t)n g (5) {n g }

(6)

Here V is the total volume of the coagulating system. Now we are ready to return to the question, what is going on in our system? The point is that the concentrations appearing in the Smoluchowski equation are defined as a thermodynamic limit of the ratio n¯ g /V , (V −→ ∞, n g −→ ∞, their ratio is finite), i.e., if n¯ g are not large enough (e.g., some of them ∝ V α , α < 1), then these particles are not “seen” in the thermodynamic limit, even if they exist. Still these particles can have large masses, comparable to the mass of entire system, and thus contribute to the mass balance. In coagulating systems such particles can form spontaneously during a finite time. This is the gelation. Two approaches had been applied for considering the sol–gel transition. The first one does not refuse from the Smoluchowski description of gelling systems, i.e., it starts with the Smoluchowski equation (e.g., [9,10]). The mass deficiency appearing after the critical time is attributed to an infinite cluster (a gel) which is introduced “by hand” (its existence does not follow from the Smoluchowski equations) and serves only for restoring the mass conservation. The gel can be assumed to be either passive or active with respect to the sol fraction (defined as the collection of particles whose population numbers are macroscopically large). In the former case the gel grows due to a finite mass flux toward infinite particle sizes. The active gel grows, in addition, by consuming the sol particles. These two mechanisms are thoroughly discussed in Ref. [10], where the reader will find references to earlier works. More accurately (but again, within the Smoluchowski scheme) this process had been considered in Refs. [11,12], where an instantaneous sink of particles with masses exceeding a large one, G, was introduced. Then the kinetics of coagulation can be described by a truncated Smoluchowski equation; no paradox with the total mass concentration appears, for the mass excess is attributed to the deposit: the particles with masses g > G consumed by the sink. Nevertheless, the difference between gelling and non-gelling systems manifests itself in the fact that during the whole pregelation period the mass is almost conserved, and only right after the critical time does a noticeable mass loss appear (the description of this model is given in this paper). An alternative approach applying the Marcus scheme [13] to the gelation problem had appeared earlier in Refs. [14– 17]. The idea of this approach relies upon the consideration of finite coagulating systems. As has been mentioned above, this approach operates with the occupation numbers and the probability for the realization of a given set of occupation numbers. Within this scheme the gel manifests itself as a narrow hump in the distribution of the average particle numbers over their masses. This hump appears after the critical time at macroscopically large mass g ∝ M and behaves itself like the active gel. The master equation governing the time evolution of the probability is extremely complicated, but on replacing it by another one, for the generating functional of the probability

A.A. Lushnikov / Physica D 222 (2006) 37–53

W , it acquires a similarity with the Schr¨odinger equation for interacting quantum Bose fields. Although many features of the solution to the evolution equation for the generating functional had been clear almost three decades ago, only very recently was I able to find the exact solution to this equation in a closed form and to analyze the behavior of the particle mass spectrum in the thermodynamic limit [18–21]. This paper does not review all achievements in the theory of sol–gel transitions. I am focusing on the exact results for the kernel equation (3) found within the Markus scheme. The point is that this approach has not found (yet, I hope) a wide recognition as an effective tool for studying the sol–gel transition, mainly because of its complexity. I will attempt to demonstrate below that the results and the technique itself are so impressive that they are deserving of more attention from the scientific audience involved in the problems related to coagulation. The remainder of this paper is organized as follows. In the next section one finds a sufficiently detailed introduction to the theory of finite coagulating systems and a sketch of the solution to the evolution equation for the generating functional (without tedious details). In Section 3 this approach is applied to the already very old, but still magnificent, problem of time evolution of random graphs. The thermodynamic limit of the solution found in Section 2 is analyzed in Section 4, where the picture of gel formation is described in detail. Section 5 reproduces the solution of the gelation problem found long ago in [11,12] within a truncated version of the Smoluchowski equation. It is shown that this model gives a reasonable picture of this phenomenon. Concluding, Section 6 discusses (qualitatively) possible scenarios of the sol–gel transition. Two Appendices introduce special families of polynomials entering the final expression of the particle mass spectrum. 2. Coagulation in finite systems The description of the coagulation process in terms of occupation numbers (numbers of g-mers considered as random variables) had been first introduced by Marcus [13]. This approach had been then reformulated by me [14–21] in the form strongly reminiscent of the second quantization. Below I outline this approach. 2.1. Illustrative example The idea of this approach is very simple. Let us consider a process in which a pair of identical particles A, on colliding, produces one A-particle (the process A + A −→ A). Let there be M such particles moving chaotically in the volume V . They collide and coalesce. Two particles produce one. Exactly like in a coagulation process. The collision rate (the probability per unit time for a pair of particles to collide) is introduced as κ/V , where κ is the rate constant of the binary reaction A+ A −→ A. The rate equation for the particle number concentration c(t) describes the kinetics of the process, dc = −κc2 . dt

(7)

39

We, however, can choose an alternative route and introduce the probability W (N , t) to find exactly N particles at time t in our system. It is also easy to guess that dW (N , t) κ = [(N + 1)N W (N + 1, t) dt 2V − N (N − 1)W (N , t)].

(8)

The first term on the RHS of this equation gives the positive contribution to the rate dt W because of the coalescence of two particles ((N + 1)N /2 is the number of ways to choose a couple of coalescing particles from N +1 particles). The second term describes the negative contribution to dt W , because the particles continue to coalesce and transfer the system from the state with N particles to the state with N − 1 particles. Hence, a simple nonlinear equation (7) is replaced by a set of linear equations (8). We, however, do not stop at this point and will make a step forward. We introduce the generating function for our probability, X Ψ (z, t) = W (N , t)z N . (9) N

From Eq. (8) we can derive the equation for Ψ , V

∂Ψ κ ∂ 2Ψ = (z − z 2 ) 2 . ∂t 2 ∂z

(10)

This equation should be supplemented with the initial condition, Ψ (z, 0) = ψ0 (z),

(11)

where ψ0 (z) is a reasonable function (it should be analytical at z = 0). For example, if we fix the number of particles N0 at the beginning of the process, then ψ0 (z) = z N0 . Or the function ψ0 (z) = e N0 (z−1) corresponds to the initial Poisson distribution. Two questions immediately arise. (i) Why should we introduce such a complex scheme for describing the kinetics of the reaction, why not use Eq. (7)? and (ii) If the second scheme describes the same process as Eq. (7), then how do we derive Eq. (7) from Eq. (10)? First I answer the second question. Let us expand the RHS of Eq. (10) near z = 1, i.e., we replace z − z 2 ≈ −2ξ , where ξ = z − 1  1. We find from Eq. (10), V

∂Ψ ∂ 2Ψ = −κξ 2 . ∂t ∂ξ

(12)

Now it is easy to solve this equation, Ψ (z, t) = ec(t)V (z−1) .

(13)

On substituting this into Eq. (12) we come to the conclusion that the concentration c(t) meets Eq. (7). Thus the probability has Poisson’s form. This approximation works well at very large N . Now let us return to the first question. If we want to describe a finite system, then eventually we should use Eq. (8) or, better, (10). Because the description of a gel requires us to make a step

40

A.A. Lushnikov / Physica D 222 (2006) 37–53

beyond the scopes of the thermodynamic limit I will use this very approach, although it requires much more serious efforts for operating and understanding the final results.

xl+m

2.2. Basic equation Let there be M monomers in a volume V . The monomers move, collide, and coalesce producing dimers, trimers, etc., along the scheme given in equation (1). Let then Q = {n 1 , n 2 . . . n g . . .}

(14)

be the state of the system given by the set of integers n g , the numbers of particles of mass g, with g being the number of monomeric units in a g-mer. A single coagulation act (the collision of two particles + their coalescence) changes a preceding state Q = {n 1 , . . . nl + 1, . . . n m + 1, . . . n g − 1 . . .} −

(15)

if l 6= m and Q − = {n 1 , . . . nl + 2, . . . . . . n g − 1 . . .}

(16)

if l = m, 2l = g, to the state Q by coalescing the particles with masses l and m to one particle with mass g. In its turn, a next coagulation act transfers the state Q to the state Q + according to the scheme Q − −→ Q −→ Q + ,

The equation for Ψ is readily derived from Eq. (19). Indeed, noticing that

(Q + )− = Q.

(17)

The probability per unit time for two particles to collide and to coalesce is K (l, m)/V , where K (l, m) is the coagulation kernel (the efficiency of the coagulation process). The rate of the process Q − −→ Q is then K (l, m) nl (Q − )(n m (Q − ) − δl,m ). (18) 2V Here δa,b stands for the Kronecker delta. The combinatorial multiplier here is just the number of ways to get a successfully coalescing pair of l- and m-mers. Next, we introduce the probability W (Q, t) to find the system in the state Q at time t. We can write down the master equation for the probability W (Q, t). It is X dW (Q, t) A(Q, Q − )W (Q − , t) = dt − Q X − A(Q + , Q)W (Q, t). (19) A(Q, Q − ) =

Q+

∂2 + X Q = (nl (Q)n m (Q) − δl,m n m (Q))X Q ∂ xl ∂ x m

and xl x m

∂2 X Q = (nl (Q)n m (Q) − δl,m n m (Q))X Q ∂ xl ∂ x m

we find instead of Eq. (19), V

∂Ψ ˆ , = LΨ ∂t

(22)

where the evolution operator Lˆ is defined as ∂2 1X K (l, m)(xl+m − xl xm ) Lˆ = . 2 l,m ∂ xl ∂ x m

(23)

At this step we introduce the operators of occupation numbers, total particle number, and total particle mass. They are: X X ∂ nˆ l = xl , Nˆ = nˆ l , Mˆ = l nˆ l . (24) ∂ xl l l Any average value of interest can be expressed in terms of Ψ . For example, the average mass spectrum n¯ g (t) = P Q n g (Q)W (Q, t) is n¯ g (t) = nˆ g Ψ (X, t) X =1 . (25) It can be readily checked that Mˆ commutes with the evolution operator, ˆ M] ˆ = 0, [L,

(26)

which means total mass conservation. The analogy with the second quantization and theory of interacting Bose fields is now clearly seen: the operator ∂x acts like an annihilation operator and x like a creation operator. Their commutator ∂x x − x∂x = 1. Hence, the first term on the RHS of Eq. (23) replaces two particles with masses l and m by one with the mass equal to l + m. The second term written as nˆ l nˆ m − nˆ l δm,l removes the pair of particles with masses l and m from a given state of the coagulating system. The rate of both these processes is K (l, m). Eq. (22) reminds us of the Schr¨odinger equation. 2.3. Exact solution

Technically it is much more convenient to deal with the generating functional X Ψ (X, t) = W (Q, t)X Q (20)

Considerable simplifications of Eq. (22) arise for the kernels K (g, l) = g f (l) + l f (g), with f being an arbitrary function. In particular, for f (x) = x one finds (see also Refs. [15,19]),

Q

Lˆ =

X

lmxl+m

X ∂2 + l 2 nˆ l − Mˆ 2 . ∂ xl ∂ x m l

(27)

rather than the probability W (Q, t). Here X stands for the set n (Q) n (Q) x1 , x2 . . . and X Q = x1 1 x2 2 . . .. P The normalization of W (Q, t) to unity, i.e., Q W (Q, t) = 1, corresponds to

If we work with the functionals belonging to a given mass M, i.e.,

Ψ (X = 1, t) = 1.

ˆ M = MΨ M , MΨ

(21)

l,m

(28)

41

A.A. Lushnikov / Physica D 222 (2006) 37–53

then the operator Mˆ in Eq. (27) can be replaced by c-number M. The functional Ψ M (X, t) can be now constructed in the form " # ∞ X Ψ M = M!Coefz z −M−1 exp z g ag (t)x g . (29)

It is convenient to introduce the variable τ=

t . V

The solution to Eq. (35) is then readily found as a formal series,

g=1

The operation Coef in Eq. (29) is defined as follows: X Coefz bm z m = b−1 .

(37)

D(z, t) =

∞ g X z g2 τ e . g! g=0

(30)

We thus find from Eq. (38)

The operation Coef is of great use when one deals with formal variables and divergent series. This operation possesses many properties of ordinary residues [22]. The reader can find the rules of the game with the formal divergent series in the review article [23]. In principle, there is nothing special in these rules. No error arises if one performs exact operations and does not try to do numerical calculations with functions of a formal variable. Every theorist is well familiar with such functions. For example, the functions of operator variables are often met in many branches of theoretical physics. It is easy to check that

Coefz (z −M+g−1 eG(z,τ ) ) =

(38)

m

• the functional Ψ M Eq. (29) corresponds to the initially monodisperse sol Ψ M (X ) = x1M if ag (0) = δg,1 , • the functional Ψ M meets Eq. (28), • the functional Ψ M is the solution to Eq. (22) with Lˆ given by Eq. (27) if the coefficients ag (t) are determined from the set of equations, g−1

V

X dag = l(g − l)al ag−l − Mgag + g 2 ag , dt l=1

(31)

• the functional Ψ M satisfies the normalization condition Ψ M (1, t) = 1. Let us introduce the generating function G(z, t) for ag (t), G(z, t) =

∞ X

ag (t)z g .

(32)

1 2 e(g −Mg)τ . (M − g)!

(39)

In order to restore ag (t) we use the Knuth identity [24] (see Eq. (A5)), ln

∞ X

x n(n−1)/2

n=1

∞ X zn zn = (x − 1)n−1 Fn−1 (x) n! n! n=1

(40)

which can be rewritten as ln D(z, τ ) =

∞ g Mgτ X z e (e2τ − 1)g−1 Fg−1 (e2τ ), g! g=1

(41)

where Fg (x) are the Mallows–Riordan polynomials discussed in [26,24,25] (see also Appendix A). Combining this result with Eqs. (32), (34), (36), (39) and (41) we finally restore the exact average particle mass spectrum,   M (g2 −2Mg+g)τ 2τ n¯ g (τ ) = e (e − 1)g−1 Fg−1 (e2τ ). (42) g The polynomials Fg (x) were also discussed in great details in Ref. [19], where one finds the derivation of several recurrence relations, the integral equation for generating function Fg (x), and their asymptotic form at large g and the argument close to unity. The reader can find a much more perfect introduction to the theory of these polynomials at the end of this article in Appendix A.

g=1

Multiplying both sides of Eq. (31) by z g and summing over all g gives the closed equation for G(z, t),   ∂G ∂G 2 ∂G ∂ ∂G V = z − Mz +z z (33) ∂t ∂z ∂z ∂z ∂z

2.4. Coagulating mixtures

with the initial condition G(z.0) = z. The substitution

with the initial condition D(z, 0) = ez . This is the central moment of the solution. According to Eqs. (25) and (29) the average occupation numbers are expressed as

where the notation m, n stands for an (m + n)-mer comprising m and n monomers of the first and the second component respectively. The most evident example is coagulation of a binary mixture, where a gas of M + N monomers of two sorts begins to form mixed clusters containing m and n monomers of each sort. This type of coagulation had been considered almost three decades ago in Ref. [27]. As has been shown in Ref. [29], coagulation in mixtures can also lead to gelation. In Refs. [21,28] an approach very similar to that described in the preceding section has been applied to the coagulating mixture with the kernel

n¯ g (t) = M!ag (t)Coefz −M+g−1 eG(z,t) .

K (k, l; p, q) = κ(kq + lp).

−Mt/V

G(z, t) = ln D(ze

, t)

(34)

allows Eq. (33) to be cast into the linear equation for D(z, t), V

∂D ∂ ∂D =z z ∂t ∂z ∂z

(35)

(36)

A wide variety of phenomena of quite different nature can be described as an aggregation process of the type (k, l) + ( p, q) −→ (m, n) = (k + p, l + q),

(43)

(44)

42

A.A. Lushnikov / Physica D 222 (2006) 37–53

Now we introduce the probability W (Q, t) to find the graph in the state Q at time t and its generating functional, Y n(g,ν|Q) X W (Q, t) x g,ν , (49) Ψ (X, t) =

Here we cite only the final result. It is    Mg Mr n¯ m,n (t) = emnt/V −m Mr t−n Mg t m n × (et/V − 1)m+n−1 Fm−1,n−1 (et/V ), (45) P where n¯ m,n (t) = ν n¯ m,n;ν (t), Mg , Mr are the numbers of monomers of the first and the second kind in the particle, Mg,r = Mg,r /V , and the polynomials Fm,n (x) play the same role as Mallows–Riordan polynomials Fg (x) in Eq. (42). The properties of these polynomials are discussed in details in Appendix B. 3. Time evolution of a random graph As an application of the above approach we consider the famous problem on the time evolution of a random graph [30]. So far the approaches to the solution of this problem used discrete methods based on combinatorial approaches [30–32]. The results presented below have recently appeared in Ref. [33]. Although the link between evolution of a random graph and coagulation in the system with the kernel equation (3) had been known (I cannot say who first mentioned this analogy, but I heard about it for the first time in 1977 from my colleague E.R. Domilovskii. The first solution (not full) similar to that presented below had been given in my Doctoral dissertation as an Appendix [34]). Let there be a graph of order M comprising N linked components. Each linked component of the graph can be characterized with its order g, the number of vertices, and the degree of filling ν (the number of edges in the component). It is clear that g − 1 ≤ ν ≤ g(g − 1)/2. The minimal value of ν corresponds to a tree of order g and the maximal one is the number of edges in the complete graph of order g. A bare vertex is also regarded as a linked component of order 1. Any state of the graph can be given by the set of the population numbers Q = {n g,ν },

(46)

(47)

or to filling a given linked component with one extra-edge, (g, ν) −→ (g, ν + 1).

g,ν

where n(g, ν|Q) stands for the occupation number n g,ν belonging to a given state Q and X = {x g,ν } is the set of independent formal variables. The functional Ψ obeys the equation ∂Ψ = (Lˆ f + Lˆ c )Ψ , (50) ∂t where the multiplier V defines the scale of time. In the theory of coagulation V is the total volume of the system. Here V depends on our good will. The RHS of this equation contains two differential operators, Lˆ f and Lˆ c . The operator Lˆ f is responsible for the evolution of the filling of linked clusters  X  l(l − 1) ∂ − λ + 1 xl,λ Lˆ f = 2 ∂ xl,λ−1 l,λ    ∂ l(l − 1) − λ xl,λ . (51) − 2 ∂ xl,λ V

The coalescence operator has the form of Lˆ in Eq. (27), 1 X ∂2 Lˆ c = . lm(xl+m,λ+µ+1 − xl,λ xm,µ ) 2 l,λ;m,µ ∂ xl,λ ∂ xm,µ

(52)

Eq. (50) replaces the master equation governing the time evolution of the probability W (Q, t). As the order of the initial graph conserves X M= lnl,λ , (53) l,λ

the solution to the evolution equation (50) can be found in the form ! X Ψ = M!Coefz z −(M+1) exp ag,ν (t)x g,ν z g , (54) g,ν

where n g,ν is the number of linked components of order g having exactly ν edges (g, ν-component). Let us consider an initially empty graph of order M (simply M bare vertices) and begin to add to this graph the edges (one edge at a time) linking two valent vertices. This process gives rise to either a coalescence of two linked components, (l, λ) + (m, µ) −→ (l + m, λ + µ + 1)

Q

(48)

The graph thus evolves due to changing the number of linked components, their order, and their degree of filling. The efficiency of the coalescence process is proportional to lm, the number of ways to connect two linked components by an extra-edge. The efficiency for filling a linked component is proportional to g(g − 1)/2 − ν, the number of possible vacant placements of the extra-edge.

where the multiplier M! provides the correct normalization of the generating functional, Ψ (1, t) = 1. Substituting Eq. (54) into Eq. (50) gives the system of equations for ag,ν (t):   dag,ν 1 V = g(g − 1) − ν + 1 ag,ν−1 dt 2   1 − g(g − 1) − ν ag,ν 2 1 X + lmal,λ am,µ δg,l+m δν,λ+µ+1 2 l,λ;m,µ g2 1 (55) − Mgag,ν + ag,ν . 2 2 Here δα,β stands for Kronecker’s delta. The initial condition for Eq. (55) is chosen in the form ag,ν = δg,1 δν,0 .

(56)

43

A.A. Lushnikov / Physica D 222 (2006) 37–53

It is easy to check that this initial condition corresponds to M , i.e., to the initially empty graph comprising Ψ (X, 0) = x1,0 only M bare vertices. Next, all ag,ν (t) = 0, once ν lies beyond the permitted interval (g − 1) ≤ ν ≤ g(g − 1)/2. Now let us try to solve Eq. (55). To this end we introduce the bivariate generating function for ag,ν (t), X G(z, ζ ; t) = z g ζ ν ag,ν (t). (57)

Next, we use the expansion (see Appendix A and Ref. [24]), w g−1 Fg−1 (1 + w) =

This equation reduces to a linear one for the function D(z, ζ ; t)   ζ ∂ ∂D ∂D 1 ∂D ∂D = z z − (ζ − 1) ζ + z , (59) V ∂t 2 ∂z ∂z ∂ζ 2 ∂z where D(z, ζ ; t) = exp[G(ze

Mt/2V

, ζ ; t)].

The initial condition to this equation is D(z, ζ ; 0) = ez .

(60)

Eq. (59) is readily solved by separating variables. Let X D(z, ζ ; t) = bn,κ Tκ (t, )Z n,κ (ζ )z n . n,κ

Then Tκ (t) = eκt and κ Z n,κ

  dZ n,κ n ζ 2 + Z n,κ , = n Z n,κ + (1 − ζ ) ζ 2 dζ 2

(61)

where κ is a separation constant. The solution to this equation is Z n,κ (ζ ) = bn,κ (1 − ζ )n

2 /2−κ

ζ κ−n/2 .

(62)

The function D should contain only integer powers of ζ . Hence, κ − n/2 = s, where s is a nonnegative integer. Next, the coefficients bn,κ should be chosen  from the initial condition Eq. (60). It is easy to see that bn,κ = to the result, D(z, ζ ; t) =

(n 2 −n)/2 s

∞ n X z nt/2V 2 e [ζ et/V + (1 − ζ )](n −n)/2 . n! n=0

Eq. (40) allows us to restore A g (ζ, t) = A g (ζ, t) =

. We then come

P

ν

(63)

ag,ν (t)ζ ν ,

(65)

where C g,ν is the number of labelled linked graphs of order g having ν edges. Eq. (65) is readily applied for restoring ag,ν (t). The result is ag,ν (t) =

1 −g(M−1)t/2V t/V e (e − 1)ν C g,ν . g!

(66)

Now we are ready to find the average number of linked components of order g having exactly ν edges. From Eq. (54) we have n¯ g,ν (t) = M!ag,ν (t)Coefz z −M+g−1 D(z, 1, t).

(67)

At ζ = 1, Eq. (63) gives D(z, 1; t) =

∞ n X z n 2 t/2V e . n! n=0

We thus come to the result,   M (g2 −2Mg+g)t/2V t/V n¯ g,ν (t) = e (e − 1)ν C g,ν . g

(68)

(69)

Because the function D(z, 1; t) coincides with the D-function corresponding to the spectrum of coagulating particles in the system with the kernel K (g, l) = gl considered above (Section 3), we can easily derive the expression for the spectrum of linked P components over their orders (numbers of vertices), n¯ g = ν n¯ g,ν . The result has the same form as the average particle mass spectrum in the coagulating system with the coagulation kernel K (g, l) = gl,   M (g2 −2Mg+g)t/2V t/V n¯ g (t) = e (e − 1)g−1 Fg−1 (et/V ). (70) g Similar results can be derived for random bipartite graphs. In this case two sorts of vertices are randomly connected with edges. The edges connecting the vertices of the same sort are forbidden. In this case the number of ways for two linked components containing m 1 , n 1 and m 2 , n 2 vertices of the first and the second sort to coalesce is proportional to m 1 n 2 + m 2 n 1 . The number of vacant places for edges in the m, n-component is mn − ν. The minimal value of ν is ν = m + n − 1. Using the results of [21], it is possible to write down the expressions for the average numbers of m, n; ν-components,    Mg Mr n¯ m,n;ν (t) = emnt/V −m Mr t−n Mg t m n × (et/V − 1)ν Cm,n;ν

(71)

and for the number of m, n components,    Mg Mr n¯ m,n (t) = emnt/V −m Mr t−n Mg t m n

1 −g(M−1)t/2V t/V e (e − 1)g−1 ζ g−1 g! × Fg−1 (ζ et/V + 1 − ζ ).

C g,ν w ν ,

ν=g−1

g,ν

The summation on the RHS of the above equation goes over all g and ν. The equation for G(z, ζ ; t) immediately follows from Eq. (55), " #  ∂G 2 ∂ ∂G ∂G ζ z +z z V = ∂t 2 ∂z ∂z ∂z   ∂G M ∂G 1 ∂G − (ζ − 1) ζ − z + z . (58) ∂ζ 2 ∂z 2 ∂z

g(g−1)/2 X

(64)

× (et/V − 1)m+n−1 Fm−1,n−1 (et/V ).

(72)

44

A.A. Lushnikov / Physica D 222 (2006) 37–53

P Here n¯ m,n (t) = ν n¯ m,n;ν (t), Mg , Mr are the numbers of vertices of the first and the second sort in the graph, Mg,r = Mg,r /V , Cm,n:ν is the number of m, n; ν-linked components, and the polynomials Fm,n (x) introduced in [21] play the same role as Mallows–Riordan polynomials Fg (x) in the expression for the number of linked component of order g in Eq. (70). In particular, the polynomials w m+n−1 Fm−1,n−1 (1 + w) generate Cm,n;ν , w m+n−1 Fm−1,n−1 (1 + w) =

mn X

Cm,n;ν w ν .

(73)

For n(t) we have dt n = −1 + µ2c (t).

On the other hand, we can derive the equation for N (t) from Eqs. (22), (25) and (27). On applying the operator Nˆ to both sides of Eq. (22), one finds dN = −M 2 + ϕ2 , (79) dt where ϕ2 is the second moment of the particle mass spectrum, V

ν=m+n−1

The polynomials themselves satisfy the identity (see Appendix B) ! ∞ X ξ m ηn X ξ m ηn mn ln x = (x − 1)m+n−1 Fm−1,n−1 (x) m!n! m!n! m,n m,n=1 + ξ + η.

(74)

This is the extension of the Knuth identity equation (40). 4. Thermodynamic limit 4.1. The Smoluchowski spectrum The sol–gel transition happens at finite t. Therefore we explore the asymptotic behavior of the mass spectrum given by Eq. (42) at finite t and V, M −→ ∞, M/V = m < ∞ (the thermodynamic limit). At large V the argument of Fg−1 in Eq. (42) approaches unity. So we try to analyze Eq. (42) at τ ∝ M −1  1. We begin the analysis by considering the limit of finite g and M, V −→ ∞. Here and below we put m = M/V = 1. The replacement t −→ mt restores the dependence of the results on m. As is known [26,19], Fg (1) = (g + 1)g−1 . We thus can write down the mass spectrum at g  M, n¯ g (t) ≈ n¯ (s) g (t) = M

g g−2 −2gt e (2t)g−1 . g!

(75)

(s)

The notation n¯ g (t) stands for the exact solution to the Smoluchowski equation (2). No traces of the catastrophe at (s) t = tc are yet seen. All functions n¯ g (t) are well defined at (s) all t. However, the total mass of the spectrum n¯ g (t) does not conserve at t > tc . It is easy to show [15,19] that the deficit of the mass concentration after the critical time tc is 2t =

1 1 ln µc (t) 1 − µc (t)

or

µc = 1 − e−2µc t .

(76)

This equation has only one root µc (t) = 0 at t < tc and two roots at t > tc . It is clear why should we choose the positive nonzero root after the critical time. The point is that the spectrum equation (72) shrinks after the critical time, so its mass cannot remain constant. The total particle number concentration can also be found, n(t) = m(t) − tm 2 (t), where m(t) = 1 − µc (t).

(77)

(78)

ϕ2 (t) =

∞ X

g 2 n¯ g (t).

g=1

Comparing this with Eq. (78) yields ϕ2 = µ2c (t). (80) V2 Eq. (80) can hold only if at t > tc the mass spectrum n¯ g (t) has a narrow peak at g = Mµc (t). This result evidences in favor of the following facts: lim

• A giant particle with mass g(t) = µc (t)M forms after the critical time tc . • In the limit M −→ ∞ the peak in the particle mass distribution has the zeroth width (the squared gel mass is equal to the second moment of the particle mass distribution). 4.2. Particle mass spectrum At large (but finite) M and finite t (t/M  1) it is more convenient to replace the polynomials Fg (x) by Pg (δ) = Fg (1 + δ) (see Appendix A). As follows from Eq. (42), δ = 2t/M. In order to investigate the particle mass spectrum we need to know the asymptotical behavior of Pg (δ) as g −→ ∞, δ −→ 0, and gδ < ∞. It is (see Eq. (A18) and [18,19]) Pg (δ) ∝

sinhg (gδ/2) . (δ/2)g

(81)

Let us introduce µ = g/M and exponentiate the exact mass spectrum equation (42), n¯ g (t) = e M Φ (µ,t) .

(82)

The function Φ, an analogue of the free energy in statistical mechanics, has the form Φ(µ, t) = −(1 − µ) ln(1 − µ) − µ ln µ + 2(µ2 − µ)t + µ ln(1 − e−2µt ).

(83)

In the vicinity of the transition point we expand Φ(µ, t) (Eq. (83)) in powers of µ and  = tc − t. The result is µ(µ − 4)2 . (84) 8 Next, we will follow the route adopted from the theory of phase transitions, with the functions Φ(µ, ) playing the role of free energy. It is easy to see that Φ(µ, t) has a minimum at Φ(µ, ) = −

45

A.A. Lushnikov / Physica D 222 (2006) 37–53

µ+ = 4/3 and a maximum at µ− = 4. It is important to notice that Φ1 (µ− , t) = 0. The mass distribution in the variables g,  has the form (see also Refs. [18,19])   g3 g2 2 n¯ g (t) = C(g, ) exp − +  − 2g . (85) M 8M 2 Unfortunately, our asymptotic analysis does not for restoring the normalization factor C(g, ). Still conclusions on its form can be retrieved from the conservation, √ M θ () C(g, ) = p , +√ 2π M 2πg 5

allow some mass

(86)

with θ () being the Heaviside step-function. Indeed, below the transition point the total mass conserves and the asymptotic mass spectrum is known. Eqs. (85) and (86) reproduce the latter at g  M. Above the transition point the second term normalizes the peak appearing at g = µ− M to unity. Now it becomes possible to describe what is going on. Below the transition point (at  < 0) the mass spectrum exponentially decreases with increasing g. The terms containing M in the denominators (see Eq. (85)) play a role only at g ∝ M. At these masses the particle concentrations are exponentially small. In short, in the thermodynamic limit and at  < 0 the first two terms in the exponent on the RHS of Eq. (85) can be ignored. The spectrum is thus given by the equation n¯ g (t) = p

M 2πg 5

2

e−2g .

(87)

At the critical point (t = tc or  = 0) the spectrum acquires the form n¯ g (t) = p

M 2πg 5

e−g

3 /8M 2

.

(88)

Although the expression in the exponent contains M in the denominator we have no right to ignore it, for this exponential factor provides the convergence of the integral for the second P moment φ2 = M −1 g 2 n¯ g in the limit M −→ ∞. We thus have Z M −g3 /8M 2 dg 1 e 1 φ2 (tc ) = √ ≈ √ Γ (1/6)M 1/3 . (89) √ g 3 π 2π 0 Here Γ (x) is the Euler gamma-function. It is seen that φ2 (tc ) does not diverge anymore. It remains finite, but now it contains the factor M 1/3 (compare with Eq. (4)). Above the transition point the situation drastically changes. Fig. 2 clearly demonstrates what is going on. Right after the critical time the particle mass distribution splits into two parts: the thermodynamically populated one whose behavior is described by the Smoluchowski equation at g  M and a narrow peak with the mass exactly equal to the difference Mµc (t) = M − Mµs (t), where µs (t) is the mass of the thermodynamically populated fraction and µc (t) is given by

Fig. 2. The particle mass spectrum at  = 0, 0.05, 0.1. It is seen how the gel appears from nothing. (From [19]).

Eq. (76). This peak can be now referred to as gel. The gel peak has a Gaussian form. Indeed, at µ ∝ µ− Eq. (85) gives √  −M(µ−µ− )2 /2 n¯ g = √ e . (90) 2π M √ The width of the gel peak is thus Γ ∝ 1/ M. This peak is well separated from the rest of the mass spectrum when it width much exceeds gc . This happens after a very short time  ∝ 1/M. Thus born gel actively “eats” smaller particles and grows. 5. Truncated models The analytical solutions to the Smoluchowski equation are known for three coagulation kernels: K 1 (g, l) = 1,

K 2 (g, l) = g + l,

K 3 (g, l) = gl.

(91)

Another three exact solution had been added in Refs. [11,12] for the kernels: K (g, l) = K i (g, l)ΘG (g)ΘG (l).

(92)

Here K i is any of the kernels given by Eq. (91) and ΘG (l) = 1 ¯G = at l ≤ G and ΘG (l) = 0 otherwise. We also introduce Θ 1 − ΘG . The meaning of the truncated kernels equation (92) is rather clear. They describe the coagulation process in the systems with the instantaneous sink which converts the active particles g ≤ G (the sol) to a passive deposit (the gel, in what follows) whose spectrum stretches from g = G + 1 to g = 2G. The gel interacting with the sol is introduced by another type of truncation, K (g, l) = K i (g, l)[ΘG (g) + ΘG (l)].

(93)

The particles with g > G cannot interact with each other and continue to grow only by eating the sol particles.

46

A.A. Lushnikov / Physica D 222 (2006) 37–53

In what follows we introduce the notation cg− (t) = ¯ G (g) for gel cg (t)ΘG (g) for sol particles and cg+ (t) = cg (t)Θ particles.

and find zT 0 − T = zT T 0 or

5.1. Passive gel

T e−T = z.

For the kernel equation (92) the equations governing the time evolution of sol and gel have the form

This is the well-familiar tree function, so the coefficients r g are readily restored,

dcg− dt

=

g−1 G X 1X − K (g − l, l)cg−l cl− − cg− K (g, l)cl− , 2 l=1 l=1

(94)

g−1

dcg+

1X − = K (g − l, l)cg−l cl− . dt 2 l=1

(95)

g g−2 . g!

(103)

Let us connect c1 and h. Using Eqs. (97), (98) and (101) gives dτ ln(c1 / h) = −1/ h

These equations are subject to the initial condition: cg (0) = δg,1 ,

rg =

or, if we use the equality dt τ = c1 , (96)

dt (c1 / h) = −(c1 / h)2 .

i.e., the particles are initially monodisperse.

The solution to this equation is

5.1.1. Exact results for K = gl Let us begin with the case of passive gel. We first solve Eq. (94). To this end we introduce the new unknown functions νg (τ ) and the variable τ [35], Z t νg (τ ) = cg /c1 , τ = c1 (t 0 )dt 0 . (97)

c1 =

h . t The spectrum is thus cg− =

g g−2 h g g! t

cg+ =

1X (g − l)lr g−l rl 2 l=1

g−1

0

Eq. (94) then gives G X dc1 = −c1 lνl− (τ ) dτ l=1

dνg− dτ

=

(98)

g−1 G X 1X − (g − l)lνg−l νl− − (g − 1)νg− lνl− . 2 l=1 l=1

(99)

For the gel spectrum we have dcg+ dτ

=

g−1 X 1 − ¯ G (g)ΘG (g − l)ΘG (l). c1 (τ ) (g − l)lνg−l νl− Θ 2 l=1

(100) Eq. (99) is readily solved by separating variables, i.e. we seek the solution in the form νg− (τ ) = r g h g−1 (τ ), g > 1. Substituting this into Eq. (99) gives G dh X + lrl h l = 1, dτ l=1

h(0) = 0

(101)

and g−1

(g − 1)r g =

1X (g − l)lr g−l rl . 2 l=1

We introduce the generating function T (z) =

∞ X g=0

gr g z g

(102)

(104)

(105) t

Z 0

h g (t 0 )dt 0 . (t 0 )2

(106)

The problem is solved, for the equation for h is trivial (however, not analytically solvable for arbitrary G). The t-dependence of τ is immediately restored by combining Eqs. (97) and (104), dτ h(τ ) = . dt t The solution to Eq. (102) can be written down as Z h dh 0 τ= . (107) 0 0 1 − TG (h ) PG l Here TG (x) = l=1 lrl x is the truncated tree function. In contrast to T (x) the function TG (x) is a polynomial with positive coefficients and thus grows with increasing x. The root h m of the equation 1 = TG (h m )

(108)

is the maximal possible value of the function h(τ ) which is reached as τ −→ ∞. The difference between the system with no truncation and the truncated system is now clearly seen. In the former τ = 1 − e−t is restricted, τ (t) ≤ 1, and the function T (h) grows and then decreases with increasing h. At finite truncations τ grows unlimitedly as t −→ ∞ and the function TG (h) also grows on increasing its argument. It is also evident that h m exceeds the convergence radius of the expansion defining the tree function. Eq. (105) thus allows us to conclude that at sufficiently large t all concentrations cg− (t) decrease with time as 1/t.

A.A. Lushnikov / Physica D 222 (2006) 37–53

47

5.1.2. Large G Our further plan is as follows. We first find h m as a function of G, then we calculate the gel spectrum, and then we will P analyze the behavior of mass M − = g gcg− (t). From Eq. (108) we find, at large G,

If K = gl then Eq. (113) gives the same result as the Smoluchowski equation with the fixed mass concentration m(y) = 1. This scenario of gelation is discussed in [10]). The truncated model gives a chance to find the mass spectrum of the gel, but nobody has ever yet tried to do it.

h m = e−1 exp[ξ0 /G],

6. Conclusion

(109)

where ξ0 = 0.854 is the root of the equation Z 1 (eξ0 x − 1)x −3/2 dx = 2.

(110)

0

The derivation is not very complicated [12]. It requires the asymptotic form of the coefficients r g ≈ (2π )−1/2 g −5/2 and some preliminary tricks before replacing the sum by an integral in the definition of TG (x). At t < tc we can replace h(t) by its value found from the Smoluchowski equation with no truncation, h(t) = te−t . Then the sol spectrum is cg− =

g g−2 g−1 −gt t e . g!

(111)

If we calculate the mass of this spectrum at t = √ tc we find that the difference from the total mass is small (∝ 1/ G). Indeed, G X

gcg− =

g=1

∞ X

gcg− −

g=1

∞ X

gcg− ≈ 1 − √

g=G+1

2 πG

.

So we see that the mass almost conserves until t = tc and then begins to decrease as t −1 . The gel spectrum can also be calculated. The result is   1−λ 1 exp[ξ0 (1 + λ)] + 1− cg = . (112) √ t π G2 (1 + λ)2 λ Here λ is introduced as g = G(1 + λ), 0 ≤ λ ≤ 1. 5.2. Active gel It is of interest also to consider the model in which the gel particles can “eat” the particles of the active spectrum (the kernel is given by Eq. (93). In this case Eqs. (94) and (95) are replaced with another two equations, dcg− dt

=

g−1 G X 1X − K (g − l, l)cl− K (g − l, l)cg−l cl− − cg− 2 l=1 l=1 ∞ X + − − cg K (g, l)cl (113) l=G+1

dcg+ dt

=

G 1 X − K (g − l, l)cg−l cl− 2 l=g−G

+

G X

+ K (g − l, l)cg−l cl−

l=g−G−1 G X − cg+ K (g, l)cl− . l=1

The initial conditions are again given by Eq. (96).

(114)

In this paper the “pathological” coagulating systems have been considered, i.e., the systems whose development in time leads to the formation of an object that is not provided for by the initial theoretical assumptions. In our case it is the gel whose appearance breaks the hypothesis that the kinetics of coagulation can be described in terms of the particle number concentrations defined as the thermodynamic limit of the ratio (occupation numbers)/volume. The coagulating system with the kernel proportional to the product of the masses of two colliding particles has been the central object of the present study. Although the main decisive step in understanding the nature of sol–gel transition in finite systems with K ∝ gl had been done long ago [14–17], only recently was I able to find the exact solution of this salient problem [18,19]. The central goal of this paper was to introduce the readers to the main ideas of the approach which I so adore. Here I have demonstrated that this approach is well applicable to the solution of other problems, like time evolution of random graphs or gelation in coagulating mixtures. Ate first sight, the coagulation process cannot lead to something wrong. Indeed, let us consider a finite system of M monomers in the volume V . If the monomers move, collide, and coalesce on colliding, the coagulation process, after all, forms one giant particle of mass M. The concentration of this M-mer is small, c M ∝ 1/M. It is better to say, it is zero in the thermodynamic limit V, M −→ ∞, M/V = m < ∞. In other words, no particles exist in coagulating systems after a sufficiently long time. But still something unexpected goes on in gelling systems after a finite interval of time. The gel forms. Two scenarios of gelation in coagulating systems have been considered in this paper. The first one considers the coagulation process in a system of finite number M of monomers enclosed in a finite volume V . In this case any losses of mass are excluded “by definition”. The gel appears as a single giant particle of mass g comparable to the total mass M of the whole system. What happens then in the system with K (g, l) ∝ gl in the thermodynamic limit? The answer is simple, although in no way apparent. In contrast to “normal” systems, where the time of formation of a large object grows with M, a giant object with the mass of the order of M forms during a finite (independent of V and M) time tc . After t = tc this giant particle (gel) actively begins to “eat” the smaller particles. Although the probability for any two particles to meet is generally small (∝ K (g, l)/V ), in the case of g ∝ M this smallness is compensated by the large value of the coagulation kernel proportional to the particle mass M, which is, in turn, proportional to V . Hence, the gel whose concentration is zero in the thermodynamic limit can play a considerable role in the evolution of the whole system.

48

A.A. Lushnikov / Physica D 222 (2006) 37–53

The structure of the kernel is also the reason why only one gel particle can form. The point is that the time for the process (l) + (m) −→ (l + m) is short for l, m ∝ M: τ ∝ V /K (l, m) ∝ V /M 2 ∝ 1/V −→ 0 in the thermodynamic limit. Of course, the Smoluchowski equation is not able to detect particles with zero concentration. The second (and the most widespread) scenario assumes that after the critical time the coagulation process instantly transfers large particles to a gel state, the latter being defined as an infinite cluster. This gel can be either passive (it does not interact with the coagulating particles) or active (coagulating particles can stick to the gel). In the latter case the gel should be taken into account in the mass balance and no paradox with the loss of the total mass arises (see Ref. [10]). Still neither this definition nor the post-gel solutions to the Smoluchowski equation give a clear answer to the question, what is this, the gel? The situation becomes more clear on considering a class of so-called truncated models (Section 5). In these models a cutoff particle mass G is introduced. The truncation is treated as an instant sink removing very heavy particles with the masses g > G from the system. So we sacrifice the mass conservation from the very beginning. The particles whose mass exceeds G form a deposit (gel) and do not contribute to the mass balance. Of course, the total mass of the active particles plus deposit conserves. The time evolution of the spectrum of active particles (with masses g < G) is described by the Smoluchowski equation as before, with the limit ∞ in the loss term being replaced with the cutoff mass G. The set of kinetic equations then becomes finite and no catastrophe is expected to arise. We have shown that, indeed, nothing wrong happens even for the coagulation kernel K ∝ gl. The total mass concentration of active particles decreases with time, as should be the case, because the largest particles settle out to deposit. But as G −→ ∞ the total mass concentration of active particles almost conserves at t < tc and only after the critical time (tc − t ∝ G −1/2 ) does the deposit begin to form and the mass to decrease with time. The question immediately arises, what kernels are pathological? In 1973 in [35] I had tried to answer this question. A primitive analysis (see also [36]) shows that for homogeneous kernels K (ag, al) = a λ K (g, l) the self-preserving asymptotic solution to the Smoluchowski equation should have the form cg (t) ≈ t −2/(1−λ) ψ(gt −1/(1−λ) ).

(115)

This asymptotic formula shows that something wrong should happen at λ > 1. So the kernels with λ > 1 occurred under suspicious. This opinion has survived until now. People continue to attack this problem, but the problem remains too hard. What is known by now? For the kernels K α (g, l) = g α l α it had been proved that the sol–gel transition exists [37,38]. At α > 1 a gel appears already from the very beginning of the coagulation process (tc = 0). So we know very little. The model K ∝ gl considered in detail here admits the exact solution. As has been shown, this solution is not so simple, especially if one tries to consider a finite system. More general models are less pleasant in this respect. Still the approaches

described above can help to answer quantitatively what is going on in more general systems. For example, the truncated models can be analyzed numerically. We can find the time behavior of particle mass concentration, to detect the gelation point (there is the chance that the mass ceases to conserve very near the critical point) and to try to look for the solution decreasing as t −1 (post-critical behavior). Moreover, such attempts have been reported (see Ref. [39]). Acknowledgements This short review of my recent results appeared, first of all, owing to Wolfgang Wagner who practically returned me to the “coagulation life” which I left many years ago. He invited me first to Berlin and then to Oberwolfach, where I met a nice, although narrow, community of scientists involved in the problems of nonlinear processes. At that time Philippe Chassaing prompted me that my very old work can be completed if I use the Mallows–Riordan polynomials. I was happy, but the nut appeared to be too hard, and I cracked it in three years (not completely), right before the second workshop at Cambridge (2003). And at last, this nice workshop in Edinburgh. I am so happy that I met many people who want to answer the same questions as me. I wish to express my deepest gratitude to the organizers of the International Workshop on Coagulation and Fragmentation 2005 for the financial support, without which my participation at the Workshop would have been impossible. Appendix A. Polynomials Pg (δ) Here I give a short introduction to theory of Mallows–Riordan polynomials appearing in the expressions for the mass spectra of coagulating particles. By definition the polynomials Pg−1 (δ) are introduced as the generating functions for the numbers C g,ν of linked labelled graph of order g,

δ g−1 Pg−1 (δ) =

g(g−1)/2 X

C g,ν δ ν .

(A1)

ν=g−1

The limits of summation in this equation are just a minimal number of edges in a tree of order g (lower limit) and the maximal number of edges g(g − 1)/2 in the complete graph (upper limit). The bivariate exponential generating function for the number of linked graphs with given numbers of vertices and edges is defined as follows: w(ξ, δ) =

∞ X ξg g=1

=

g!

∞ X ξg g=1

g!

g(g−1)/2 X

C g,ν δ ν

ν=g−1

δ g−1 Pg−1 (δ).

(A2)

According to the Riddel theorem [40], 1 + w(ξ, δ) = ln W(ξ, δ),

(A3)

with W(ξ, δ) being the exponential generating function of all labelled graphs. The latter is readily found. Indeed, the number

49

A.A. Lushnikov / Physica D 222 (2006) 37–53



 of ways to connect g vertices by ν edges is g(g−1)/2 . Hence, ν the polynomial (1 + δ)g(g−1)/2 generates the numbers of graphs of order g having exactly ν edges. The exponential generating function for these graphs is thus W(ξ, δ) =

∞ X ξg g=0

g!

(1 + δ)g(g−1)/2 .

(A4)

Hence, ln

∞ X ξg g=0

g!

(1 + δ)g(g−1)/2 =

∞ X ξg g=1

g!

δ g−1 Pg−1 (δ).

(A5)

This identity is normally formulated in terms of polynomials Fg (x) = Pg (x − 1) (see Eq. (40)). It is clear, that ∂ξ w(ξ, δ) = y(ξ δ, δ),

(A6)

where y(ξ, δ) =

∞ X ξg g=1

g!

Pg (δ)

(A7)

ln y(ξ δ, δ) = w[(1 + δ)ξ, δ] − w(ξ, δ) Z (1+δ)ξ y(xδ, δ)dx. =

(A8)

ξ

Replacing ξ δ = ζ and x = ζ (1 + uδ)/δ reduces Eq. (A8) to Z 1 ln y(ζ, δ) = ζ y[ζ (1 + uδ), δ]du. (A9) 0

Eqs. (A3) and (A6) allow for a derivation of a linear set of recurrence relations for the polynomials Pg (δ). The derivation applies the theorem: if ∞ ∞ X X ln ak x k = bk x k k=0

kak =

ln µ − ln(1 − µ) − µαΨ 0 (µα) − Ψ (µα) − α(µ − 1) = 0 (A14) We introduce the variable x = µα and two unknown functions, µ(x) and Ψ (x). We obtain two equations for determining these functions: x (A15) µ ln µ + (1 − µ) ln(1 − µ) − µΨ − (µ − 2) = 0 2 ln µ − ln(1 − µ) − (xΨ )0 − x(1 − 1/µ) = 0.

(A16)

It is easy to check that µ(x) = 1 − e−x ,

and

Ψ (x) = ln[2 sinh(x/2)]

(A17)

is the solution to Eqs. (A14) and (A15). It is also important to notice that if we return to the variables µ, t, then there are two solutions for the mass, µ = 0 for all t and µ = 1 − e−2µt for t > tc . From the second Eq. (A17) we finally find that at large g and small δ Pg (δ) ∝

sinhg (gδ/2) . (δ/2)g

(A18)

The first four F-polynomials are F0 (x) = F1 (x) = 1, and F2 (x) = x + 2,

F3 (x) = x 3 + 3x 2 + 6x + 6,

F4 (x) = x 6 + 4x 5 + 10x 4 + 20x 3 + 30x 2 + 36x + 24, and

then k−1 X

and Φµ0 = 0,

and

is the exponential generating function for Pg (δ). This fact helps us to derive the equation for y(ξ, δ). Let us differentiate both sides of Eq. (A3) over ξ and notice that ∂ξ W(ξ, δ) = W [(1 + δ)ξ, δ]. The latter equality immediately follows from Eq. (A4). The result of these operations is

k=0

Eq. (A11) allows us to derive the asymptotic formula for Pg (δ) in the limit of large g and small δ. The idea is very simple. The maximal value of each term in the sum on the RHS of Eq. (A12) should be of the order of unity. If we represent each term in the exponential form e−gΦ and introduce µ = m/g, α = gδ, and Ψ = ln[δ m Pm (δ)] we find that α Φ = µ ln µ + (1 − µ) ln(1 − µ) − µΨ (µα) − µ(µ − 2) 2 =0 (A13)

P2 (δ) = 3 + δ, (k − s)as bk−s .

(A10)

P3 (δ) = 16 + 15δ + 6δ 2 + δ 3

P4 (δ) = 125 + 222δ + 205δ 2 + 120δ 3 + 45δ 4 + 10δ 5 + δ 6 .

s=1

Applying this theorem to Eqs. (A3) and (A4) yields, after a simple algebra, g   X g m=0

m

δ Pm (δ)(1 + δ) m

(m+1)(m−2g)/2

= 1.

(A11)

Another recurrence relation arises if we apply this theorem to Eq. (A8),  g  1 X g−1 Pg (δ) = [(1 + δ)m − 1]Pm−1 (δ)Pg−m (δ). δ m=1 m − 1 (A12)

Appendix B. Polynomials Pm,n (δ) The polynomials Pm−1,n−1 (δ) are introduced similarly. By definition δ m+n−1 Pm+n−1 (δ) is the generating function for the numbers Cm,n;ν of a linked labelled bipartite graphs of order m, n, δ m+n−1 Pm−1,n−1 (δ) =

mn X

Cm,n;ν δ ν .

(B1)

ν=m+n−1

The limits of summation in this equation are just the minimal number of edges in a tree of order m, n (lower limit) and the

50

A.A. Lushnikov / Physica D 222 (2006) 37–53

maximal number of edges mn in the complete graph (upper limit). The exponential generating function for the number of linked bipartite graphs is defined as follows: ∞ X ξ m ηn m+n−1 w(ξ, η; δ) = δ Pm−1.n−1 (δ). m!n! m,n=1

(B2)

Again, according to the Riddel theorem [40], w(ξ, η; δ) = ln W(ξ, η; δ),

X ξ m ηn m,n

m!n!

(1 + δ)mn .

(B4)

∞ X ξ m ηn m+n−1 δ Pm−1,n−1 (δ). m!n! m,n=1

(B5)

Let us now introduce x = 1 + δ. Then Eq. (B4) takes the form W(ξ, η; x) =

X ηn n

n!

e

ξ xn

=e

ξ

X ηn n

n!

e

ξ(x n −1)

.

(B6)

Hence, the function # " k X ηn ξ(x n −1) w(ξ, η) = ξ + ln e n! n=0

W∂η w = ∂η W.

(B9)

Substituting here expansions (B4) and (B5) gives m,n     X m n x − pn−qm+ pq c p+1,q (x) p q p,q=0 =

m,n     X m n x − pn−qm+ pq c p,q+1 (x) = 1. p q p,q=0

(B10)

These identities are of use in deriving the expression for the composition spectrum in the thermodynamic limit. Now let us derive the set of recurrence relations for cm,n (x) Eq. (B15), it is easy to see that ∂ξ W(ξ, η; x) = W(ξ, xη; x)

∂η W(ξ, η; x) = W(xξ, η; x).

(B7)

Hence ∂ξ w(ξ, η, x) = ew(ξ,xη;x)−w(ξ,η;x), and ∂η w(ξ, η, x) = ew(xξ,η;x)−w(ξ,η;x).

Combining Eqs. (B12) and (B4) yields " # X X ξ m ηn n ξ m ηn = exp cm,n (x) (x − 1) . cm+1,n (x) m!n! m!n! m,n m,n Next, we apply Eq. (B12) for finding ∂ 2 w/∂ξ ∂η,

=

For example, ∂η w(ξ, η) η=0 = e(x−1)ξ

cα,β+1

X ξ α ηβ β+1 ξ γ ηδ (x − 1) cγ +1,δ . α!β! γ !δ! γ ,δ

(B14)

The summations in the two above equations go over all nonnegative integers. Eq. (B14) is equivalent to the set of recurrence relations   m,n  X m n cm+1,n+1 = cm−q, p+1 m−q n−p p,q=0

and thus cm,1 = (x Next, 2 2 = e(x −1)ξ − 2e2(x−1)ξ. ∂η,η w(ξ, η) − 1)m .

η=0

This function generates all cm,2 : cm,2 = (x − 1)m [(x + 1)m − 2m ].

× cq+1,n− p (x p+1 − 1).

It is not so difficult to continue in this spirit: m

m

cm+1,n+1 (x) = (x − 1)m+n+1 Fm,n (x),

In particular, we find c23 = (x − 1)4 (x 2 + 4x + 7),

c33 = (x − 1)5 (x 4 + 5x 3 + 15x 2 + 29x + 31).

(B15)

At this stage we introduce the polynomials Fm,n (x). They are defined as

cm,3 = (x − 1) [(x + x + 1) − 3(x + 2) + 2 · 3 ]. c2,2 = (x − 1)3 (x + 3)

X α,β

m

(B12)

X ∂ 2w ξ m ηn = cm+1,n+1 ∂ξ ∂η m!n! m,n

cm,n (x) = (x − 1)m+n−1 Pm−1,n−1 (x − 1).

2

(B11)

(B13)

generates all cm,n (x) with n ≤ k, where

m

and

and

Here the summation goes over all nonnegative integers m, n except m = n = 0. Hence, X ξ m ηn (1 + δ)mn ln m!n! m,n =

W∂ξ w = ∂ξ W

(B3)

with W(ξ, η; δ) being the exponential generating function of all labelled bipartite graphs. The latter is readily found. Indeed, the
number of ways to connect m, n vertices by ν edges is mn mn generates the numbers ν . Hence, the polynomial (1 + δ) of graphs of order g having exactly ν edges. The exponential generating function for these graphs is thus W(ξ, η; δ) =

Eq. (B3) allows for deriving two sum rules for cm,n (x). It is evident that

or (B8)

Pm,n (δ) = Fm,n (1 + δ).

(B16)

51

A.A. Lushnikov / Physica D 222 (2006) 37–53

Substituting Eq. (B16) into Eq. (B15) gives the recurrence for the polynomials Fm,n (x),   m,n  X m+1 n+1 Fm+1,n+1 (x) = m+1−q n+1− p p,q=1 x p+1 − 1 . × Fm−q, p (x)Fq,n− p (x) x −1

(B17)

The lowest order polynomials are F0,0 = 0, F0,1 = F1,0 = 1, and F1,1 (x) = x + 3,

F1,2 (x) = F2,1 (x) = x + 4x + 7, 2

3

∞ X

Fm−1.n−1 (x)

(B18)

m,n=1

ξ m ηn + ξ + η. m!n!

(B19)

On the other hand, using Eq. (B16), we find Φ[(x − 1)ξ, (x − 1)η; x] = (x − 1)w(ξ, η; x).

(B20)

Eq. (B12) allows us to derive the set of equations for Φ(ξ, η; x). We introduce (x − 1)ξ = a and (x − 1)η = b. Then Eq. (B20) can be rewritten as   a b , ;x . (B21) Φ(a, b; x) = (x − 1)w x −1 x −1 Differentiating both sides of this equation with respect to a and b and using Eq. (B12) yields the set of functional equations for Φ,   ∂Φ(a, b; x) Φ(a, xb; x) − Φ(a, b; x) = exp , ∂a x −1   ∂Φ(a, b; x) Φ(xa, b; x) − Φ(a, b; x) = exp . (B22) ∂b x −1 Eq. (B22) can be rewritten in a more symmetric form, Z 1 ln X δ (a, b) = b Yδ [a, (1 + uδ)b]du, ln Yδ (a, b) = a

Z

Pm,n (0) = (m + 1)n (n + 1)m .

0

Indeed, let us introduce

u = ηev ,

and

v(ξ, η) = ln Y (ξ, η; 0).

v = ξ eu .

(B26)

Then we can find Pm,n (0) from the following chain of equalities: X (ξ, η; 0) ξ m+1 ηn+2 eu(m+2)+v(n+2) −(u+v) = m!(n + 1)!Coefu,v e (1 − uv). v m+1 u n+2

Pm,n (0) = m!(n + 1)!Coefξ,η

Elementary algebraic transformations prove Eq. (B25). The multiplier e−(u+v) (1 − uv) in the last term is just the Jacobian appearing in replacing the variables ξ, η −→ u, v. It is easy to find the first several polynomials. From Eq. (B18), we have P1,1 (δ) = δ + 4,

P1,2 (δ) = P2,1 (δ) = δ 2 + 6δ + 12,

P2,2 (δ) = δ 4 + 9δ 3 + 36δ 2 + 78δ + 81, . . . . The generating function Φ(ξ, η; 1) can be expressed in terms of X 0 (ξ, η) and Y0 (ξ, η) as follows: Φ(ξ, η; 1) = ln X 0 + ln Y0 − ln X 0 ln Y0 .

(B23)

Yδ (a, b) = ∂b Φ(a, b; 1 + δ).

It is clear that X δ and Yδ are the bivariate exponential generating functions for the polynomials Pm,n−1 (δ) = Fm,n−1 (1 + δ) and

(B27)

In order to prove this identity it is enough to differentiate Eq. (B27) with respect to η. The left-hand side gives simply Φη0 (ξ, η; 1) = Y0 (ξ, η). We will show that the right-hand side of this equation reproduces the same result. Indeed, as follows from Eq. (B23), X 0 and Y0 obey the set of transcendent equations, ln X 0 = ηY0 ,

ln Y0 = ξ X 0 .

Y0 ∂ X0 = , ∂η 1 − ln X 0 ln Y0

Here δ = x − 1 and X δ (a, b) = ∂a Φ(a, b; 1 + δ),

(B24)

(B25)

These equations allow us to find the derivatives,

0 1

X δ [(1 + uδ)a, b]du.

Pm,n−1 (δ)

It is easy to show that

The exponential generating function for the polynomials Fm,n (x) is introduced as Φ(ξ, η; x) =

X

From Eq. (B23) we find

We finally come to the central identity allowing for deriving the exact mass spectrum (45), ! ∞ X ξ m ηn X ξ m ηn mn ln x = (x − 1)m+n−1 Fm−1,n−1 (x) m!n! m!n! m,n m,n=1 + ξ + η.

ξ m ηn , m!n! X ξ m ηn Yδ (ξ, η) = Pm−1,n (δ) . m!n! X δ (ξ, η) =

u(ξ, η) = ln X (ξ, η; 0)

F2,2 (x) = x + 5x + 15x 2 + 29x + 31, . . . . 4

Pm−1,n (δ) = Fm−1,n (1 + δ),

∂Y0 Y0 ln Y0 = . ∂η 1 − ln X 0 ln Y0

Now we differentiate the right-hand side of Eq. (B27), ln0η X 0 + ln0η Y0 − ln0η X 0 ln Y0 − ln X 0 ln0η Y0 = ln0η X 0 (1 − ln ηY0 ) + ln0η Y0 (1 − ln ηX 0 ) =

Y0 (1 − ln Y0 ) + Y0 ln Y0 (1 − ln X 0 ) = Y0 . 1 − ln X 0 ln Y0

This result proves Eq. (B27).

52

A.A. Lushnikov / Physica D 222 (2006) 37–53

Now we use Eq. (B10) for deriving the asymptotic formula for the polynomials Pm,n (δ) in the limit m, n −→ ∞, δ −→ 0, mδ, nδ < ∞. The answer can be easily guessed:

It is easy to check that

Pm,n (δ) ∝ m n n m h m (nδ)h n (mδ),

meet Eq. (B38). At small ζ , we find

(B28)

where sinh x/2 h(x) = . x/2

(B29)

In order to prove Eq. (B28) we rewrite the first Eq. (B10) in the exponential form, X 0 0 eΩ (m,n, p ,q ;δ) = 1, (B30) p,q

where p 0 = p/m, q 0 = q/n, and Ω (m, n, p 0 , q 0 ; δ) = −m[(1 − p 0 ) ln(1 − p 0 ) + p 0 ln p 0 ] − n[(1 − q 0 ) ln(1 − q 0 ) + q 0 ln q 0 ] − mn( p 0 + q 0 − p 0 q 0 )δ + mp 0 H (nq 0 δ) + nq 0 H (mp 0 δ).

(B31)

The idea is that the terms of the order of unity contribute to the right-hand side of Eq. (B30). Next, these terms correspond to the maximum of eΩ . We assumed also that C p,q ≈ e p H (qδ)+q H ( pδ) , in accordance with Eq. (B28). Let us analyze the condition for the maximum of Ω with respect to p 0 . It gives −m[− ln(1 − p 0 ) + ln p 0 ] − mnδ(1 − q 0 ) + m H (nq 0 δ) + mnq 0 δ H 0 (mp 0 δ) = 0.

(B32)

Now we introduce the variables x = mp 0 δ and y = nq 0 δ and rewrite Eq. (B32) in terms of these variables, ln(1 − p 0 ) − ln p 0 −

y + y + H (y) + y H 0 (x) = 0. q0

(B33)

A similar equation holds for q 0 : ln(1 − q 0 ) − ln q 0 −

x + x + H (x) + x H 0 (y) = 0. p0

(B34)

Together with the condition Ω (m, n, p 0 , q 0 ; δ) = 0,

(B35)

we have three equations for determining p 0 , q 0 and H as functions of x and y. Eqs. (B33) and (B34) can be resolved under the condition that p 0 (x, y) = f (y) and q 0 (x, y) = f (x). In this case a separation of variables is possible and both Eqs. (B33) and (B34) give one and the same couple of equations for f (ζ ) and H (ζ ), ln(1 − f ) − ln f + ζ + H (ζ ) 1 = − H 0 (ζ ) = a. ζ f

(B36)

Here a is a separation constant. Its value determined from the condition (B35) is a = 1/2. Then we find, from Eq. (B36), f (ζ ) = 1 − e−ζ

and

H (ζ ) = ln(2 sinh(ζ /2)).

(B37)

p 0 (x, y) = 1 − e−nq (x,y)δ , q 0 (x, y) = 1 − e−mp (x,y)δ . (B38) 0

H (ζ ) ≈ ln ζ +

ζ2 . 24

0

(B39)

References [1] R.L. Drake, in: G.M. Hidy, J.R. Brock (Eds.), Topics in Current Aerosol Researches (Part II), Pergamon, New York, 1972, pp. 201–376. [2] D.J. Aldous, Deterministic and stochastic models for coalescence (aggregation, coagulation; review of the mean-field theory for probabilists), Bernoulli 5 (1999) 3–122. [3] F. Leyvraz, Scaling theory and exactly solved models in the kinetics of irreversible aggregation, Phys. Rep. 383 (2003) 95–212. [4] Z.A. Melzak, The effects of coalescence in certain collision processes, Q. J. Appl. Math. 11 (1953) 231–236. [5] J.B. McLeod, On an infinite set of nonlinear differential equations, Q. J. Math. 2 (1962) 119–124. [6] J.B. McLeod, On an infinite set of nonlinear differential equations, Q. J. Math. 2 (1962) 192–196. [7] W.H. White, A global existence theorem for Smoluchowski’s coagulation equation, Proc. Amer. Math. Soc. 80 (1980) 273–281. [8] P.G.J. van Dongen, M.H. Ernst, On the occurrence of a gelation transition in Smoluchowski’s coagulation equation, J. Statist. Phys. 44 (1986) 785–792. [9] E.M. Hendrics, M.H. Ernst, R.M. Ziff, Coagulation equation with gelation, J. Statist. Phys. 31 (1983) 519–563. [10] R.M. Ziff, G. Stell, Kinetics of polymer gelation, J. Chem. Phys. 73 (1980) 3492–3499. [11] A.A. Lushnikov, V.N. Piskunov, Three new exactly solvable models in theory of coagulation, Dokl. Akad. Nauk SSSR 267 (1982) 127–130. [12] A.A. Lushnikov, V.N. Piskunov, Analytical solutions in the theory of coagulating systems with sinks, Appl. Math. Mech. 47 (1983) 931–939. [13] A.H. Marcus, Stochastic coalescence, Technometrics 10 (1968) 133–143. [14] A.A. Lushnikov, Coagulation in systems with distributed initial conditions, Dokl. Akad. Nauk SSSR 236 (1977) 673–676. [15] A.A. Lushnikov, Some exactly solvable models of stochastic coagulation, Dokl. Akad. Nauk SSSR 237 (1977) 1122–1125. [16] A.A. Lushnikov, Coagulation in finite systems, J. Colloid Interface Sci. 65 (1978) 276–285. [17] A.A. Lushnikov, Some new aspects of coagulation theory, Izv. AN SSRR, Ser. Phys. Atmos. Ocean. 14 (1978) 1046–1054. [18] A.A. Lushnikov, From sol to gel exactly, Phys. Rev. Lett. 93 (2004) 198302. [19] A.A. Lushnikov, Exact kinetics of the sol–gel transition, Phys. Rev. E 71 (2005) 046129. [20] A.A. Lushnikov, Exact particle mass spectrum in a gelling system, J. Phys. A 38 (2005) L35–L39. [21] A.A. Lushnikov, Sol–gel transition in a coagulating mixture, J. Phys. A 38 (2005) L383–L387. [22] G.P. Egorychev, Integral Representation and Counting the Combinatorial Sums, Nauka, Moscow, 1977 (in Russian). [23] I. Niven, Formal power series, Amer. Math. Monthly 76 (1969) 871–889. [24] D.E. Knuth, Linear probing and graphs, Algorithmica 22 (1998) 561–568. [25] P. Flajolet, P. Poblete, A. Viola, On the analysis of linear probing hushing, Algorithmica 22 (1998) 490–515. [26] G. Kreweras, Une famille de polynomes ayant pluseurs proprietes enumeratives, Period. Math. Hungar. 11 (1980) 309–320. [27] A.A. Lushnikov, Evolution of coagulating systems. III. Coagulating mixtures, J. Colloid Interface Sci. 54 (1976) 94–103.

A.A. Lushnikov / Physica D 222 (2006) 37–53 [28] A.A. Lushnikov, Exact kinetics of sol–gel transition in a coagulating mixture, Phys. Rev. E 73 (2006) 036111. [29] E.R. Domilovskii, A.A. Lushnikov, V.N. Piskunov, A new exactly solvable model of a coagulating mixture, Dokl. Akad. Nauk SSSR 243 (1978) 407–410. [30] P. Erd¨os, A. R´enyi, On evolution of random graphs, Publ. Math. Inst. Hung. Acad. Sci. 5 (1960) 17–61. [31] B. Bollobas, Random Graphs, Academic, London, 1985. [32] R. Albert, A.-L. Barab´asi, Statistical mechanics of complex networks, Rev. Modern Phys. 74 (2002) 47–97. [33] A.A. Lushnikov, Time evolution of a random graph, J. Phys. A 38 (2005) L777–L782. [34] A.A. Lushnikov, Evolution of coagulating systems, Doct. Dissertation, Karpov Institute of Physical Chemistry, Moscow, 1978.

53

[35] A.A. Lushnikov, Evolution of coagulating systems, J. Colloid Interface Sci. 45 (1973) 549–559. [36] A.A. Lushnikov, M. Kulmala, Singular self-preserving regimes of coagulation–condensation process, Phys. Rev. E 65 (2002) 041604. [37] I. Jeon, Existence of gelling solutions for Coagulation–fragmentation equations, Comm. Math. Phys. 194 (1998) 541–567. [38] F. Leyvraz, Existence and properties of post-gel solutions of the equations of coagulation, J. Phys. A 16 (1983) 2861–2871. [39] M.H. Lee, A survey of numerical solutions to the coagulation equation, J. Phys. A 34 (2001) 10219–10241. [40] R.J. Riddel, G.E. Ulenbeck, On the theory of virial development of the equation of state of monoatomic gases, J. Chem. Phys. 21 (1953) 2056–2064.