The Becker-Döring equations with exponentially size-dependent rate coefficients

This paper is concerned with an analysis of the Becker-Döring equations which lie at the heart of a number of descriptions of non-equilibrium phase transitions and related complex dynamical processes. The Becker-Döring theory describes growth and fragmentation in terms of stepwise addition or removal of single particles to or from clusters of similar particles and has been applied to a wide range of problems of physicochemical and biological interest within recent years. Here we consider the case where the aggregation and fragmentation rates depend exponentially on cluster size. These choices of rate coefficients at least qualitatively correspond to physically realistic molecular clustering scenarios such as occur in, for example, simulations of simple fluids. New similarity solutions for the constant monomer Becker-Döring system are identified, and shown to be generic in the case of aggregation and fragmentation rates that depend exponentially on cluster size. PACS numbers: 64.60.-i, 82.60.Nh 05.70.Fh


Introduction
The Becker-Döring equations were originally developed to describe the kinetics of firstorder phase transitions [4] and as such have become a central element in the description of the kinetics of classical nucleation theory.The Becker-Döring equations describe the stepwise growth and fragmentation of clusters (say of condensing liquid droplets within a vapour comprised of monomer molecules) in terms of the rates of the individual processes wherein single monomer molecules join or leave each cluster.It should be noted that although the Becker-Döring model provides a detailed description of such cluster growth and fragmentation processes, it is not fully 'microscopic' in that it provides information only on the average dynamical behaviour of clusters, but does not say anything about fluctuations.In this sense, the Becker-Döring model provides us with a mean-field theory.However, the general subject of nucleation and growth processes is one that has received comparatively little attention-of either an experimental or a theoretical nature-in the past few decades, largely owing to the twin difficulties of experimental reproducibility and reliable measurement methods on the one hand, and immense theoretical complexity on the other.
As will be described below, the Becker-Döring model corresponds mathematically to an infinite set of coupled ordinary differential equations describing the temporal evolution of the cluster size distribution.For each cluster size, the time evolution of its concentration depends on two parameters, the aggregation ('forward') and fragmentation ('backward') rate coefficients so that in general one needs to have experimental data on all of them to describe the dynamical behaviour completely.Such information is unlikely ever to be available.
Faced with this impasse, we have developed two methods of analysis.In certain cases where the form of the aggregation and fragmentation rates are known and have a relatively simple form, standard asymptotic techniques can be applied to deduce the kinetics at large times and large cluster sizes.Results for the case of rates which depend on size in a power-law fashion (a case of direct relevance to classical nucleation theory) have been derived by King & Wattis [13,26].Coveney, Wattis and Bolton have shown that by eliminating a proportion of the concentration variables, a coarse-grained mesoscopic description can be obtained [10,24,25,6].The resulting model also has the form of a Becker-Döring system, with modified aggregation and fragmentation rate coefficients; the coarse-graining procedure can be reapplied and thus has the form of a renormalisation transformation, with each application filtering out more of the detailed dynamics of the cluster size distribution.
The purpose of the present paper is to introduce new sets of rate coefficients which depend exponentially on cluster size, and to show that the large-time and large cluster size asymptotic analysis is tractable.The exponential dependence of rate coefficients on size within Becker-Döring theory is far from being simply a technical curiosity; rather, it describes cluster-formation processes of wide-ranging and considerable practical interest.Preliminary work of Wonczak & Strey [27] shows using molecular dynamics that the decay rates of argon clusters scale exponentially with cluster size.
Several other authors have studied the Becker-Döring cluster kinetic equations in a variety of settings, from the general theory of cluster formation as studied by Brilliantov & Krapivsky [7], by da Costa [9] and Bolton & Wattis [5], to the specific case of late time scaling scaling laws, where the Becker-Döring model should match into the Lifshitz-Slyozov-Wagner (LSW) theory of cluster growth [15,21].Penrose [18] analysed the late time asymptotics of the Becker-Döring equations to show a connection with the LSW theory of coarsening.Velazquez [20] used a combination of asymptotics and renormalisation ideas to analyse the constant mass formulation of the Becker-Döring equations and draw connections with Lifshitz-Slyozov scaling theory.Several authors have used asymptotic methods to study the Becker-Döring equations.Penrose & Lebowitz [19] aimed to understand the separation of fast and slow dynamics in the Becker-Döring equations in order to estimate nucleation rates.
In Section 2 we summarise the various formulations of the Becker-Döring system of equations.Section 3 describes the simplest form of Becker-Döring model whose rate coefficients depend exponentially on cluster size; such models are shown to give rise both to systems with and without a nucleation barrier.We also describe the equilibrium solution, the steady-state solution and the kinetics expected to be observed in such systems, identifying the novel self-similar form present in the large-time asymptotics.
In section 4 we illustrate this solution via numerical solutions of the Becker-Döring rate equations and the appendix shows the results of numerical simulations of the coarsegrained approximation of the Becker-Döring system.The paper concludes in Section 5 with a general discussion of the theory we have developed here and its consequences for the modelling of complex aggregation-fragmentation processes.

The Becker-Döring equations
The model underlying the theory presented here is that of a constant mass system in which monomer is added to or removed from the system by way of a precursor species P which undergoes a reversible reaction to create monomer C 1 .This is described by the simple reaction P C 1 with forward rate k f (p, c 1 ) and reverse rate k b (p, c 1 ).Denoting the concentration of P by p(t) we then have the system of rate equations where the concentration of clusters C r of size r (i.e.aggregation number r) is denoted by c r .For isothermal systems, the aggregation rates a r and fragmentation rates b r+1 are assumed to depend only on cluster size.The quantity J r is the net flux of material from cluster size r to size (r + 1).In this model, clusters grow and fragment by the addition or removal of only a single particle at a time.This system has a conserved quantity, tot = p + ∞ r=1 rc r ; note however, that the amount of mass in the form of clusters ({c r } ∞ r=1 , namely BD (t) = ∞ r=1 rc r ) is not constant.There are two simplifications of this system which merit a more detailed study.The first, and the form which we study in this paper, is the constant monomer concentration formulation ) We typically take c 1 = 1 to simplify the algebra.This simplification assumes that the precursor chemical decays at a rate which maintains the monomomer concentration constant, that is, It is also valid in the case where there is an overabundance of monomer in the system, so that further increases (or even moderate decreases) in the concentration of monomer do not affect the rates of reaction (i.e. the monomer concentration is so high that any change in its concentration is negligible, the so-called 'pool chemical approximation').
The system (2.5) has various properties which make it interesting and applicable in a variety of physical and physicochemical situations.Firstly there is the total mass of the system (t) = ∞ r=1 rc r (t), (2.6) which varies due to the input and removal of monomeric material, ˙ = J 1 + ∞ r=1 J r .Provided the quantity is bounded below, it qualifies as a Lyapunov function and can be used to prove the convergence of all solutions to the equilibrium solution c r = Q r c r 1 .If V is the volume of the system then k B T VV ({c r }) is the free energy of the system (2.5).Here Q r is the partition function generated from the rate coefficients by Finally there is the identity valid for sequences {g r } ∞ r=1 which do not grow too rapidly as r → ∞ (for details see Ball et al .[3]).In cases where the equilibrium solution is divergent, the system will approach a steady-state solution, given by J r = J independent of cluster size r.The flux through the system is then given by (2.9) The second simplification of (2.1)-(2.4)corresponds to the later stages of the process, in situations when the precursor species is exhausted and no further monomer is introduced to the clustering part of the system.We formally derive this by assuming that the decay of precursor is irreversible, that is k b = 0 and p = 0.This yields the constant mass formulation of the Becker-Döring equations together with (2.3)-(2.4).This system has mass (or density) as a conserved quantity, a similar set of identities and a similar Lyapunov function 12) It has the same form of equilibrium solution as the constant monomer system, though given a set of initial data, determination of the equilibrium solution is non-trivial since the monomer concentration must be chosen in such a way that the equilibrium concentrations c r have total mass (2.6) which matches that of the initial data.Since the monomer concentration is then one of the variables to be determined, this formulation is somewhat harder to analyse than the constant monomer case.The kinetics of the Becker-Döring equations can be studied using a variety of methods.This paper will be concerned with asymptotic analysis and accompanying numerical simulations Asymptotic analysis is most suited to the constant monomer concentration formulation of the Becker-Döring equations.Where this method is less useful, namely the case of constant mass, more reliance has to be placed on numerical simulations, for example Carr, et al [8] and Duncan & Soheili [11].In these papers the authors have analysed extremely large systems of such equations, and use a variety of techniques to overcome the attendant problems.The numerical solution of such systems requires a delicate balancing of errors, since the systems respond on a wide variety of timescales.

Systems with exponential rate coefficients
In this section we investigate the physical situations which may be modelled by rate coefficients of the form We study the form of equilibrium and steady-state solutions so generated.In the next section we analyse the kinetic behaviour exhibited by models of this form.Note that for the existence theorem for the constant mass formulation of the Becker-Döring equations of Ball et al [3] to be valid we require a r (g r+1 −g r ) = O(g r ) for all r with g r+1 −g r ≥ δ > 0 for some constant δ > 0. The formula gives a r (g r+1 − g r ) = Cg r for all r.For γ < 0 this also gives a sequence satisfying g r+1 − g r ≥ δ > 0, but for γ > 0, g r → g ∞ as r → ∞ so the condition g r+1 − g r ≥ δ > 0 fails.Theorem 2.7 of [3] thus guarantees the existence of solutions for γ < 0 but shows that no solution to the constant mass formulation can have γ > 0 and η ≤ γ.We, however, will be considering the kinetics of the constant monomer system, and will restrict the parameter regime to γ ≤ 0, so that we are can be sure of the existence of an unique solution.

Equilibrium behaviour
In order to see why such a functional dependence of rates on aggregation size is physically interesting, let us examine some of the more immediately obvious properties of a system in which the rate coefficients are governed by equations (3.1).The statistical mechanical partition function is which is rapidly decaying in the large aggregation number limit provided η > γ.All possible qualitatively different partition functions Q r are plotted in Figure 1, using the parameter ) , and Q r = 1 at r = 1 and r = r c .Cases where Q r is monotonic in r ≥ 1 correspond to r c < 1.At large aggregation numbers, the equilibrium concentrations have the asymptotic form where θ = ac 1 /b.
In Figures 1(a)-(b), η < γ and b a so that Q r diverges at large cluster sizes, after first decaying from r = 1 to some critical cluster size r = r c .Thus in these cases Q r has the form of the partition function used in classical nucleation theory with a nucleation barrier; see, for example, Lewis [14].The shapes of the partition function displayed in Figures 1(c)-(d) have the familiar single-humped form encountered in the cluster distributions of real, that is polydisperse, colloidal and surfactant systems including micelles and vesicles.The polydispersity of the distribution is a function of the parameters a, b, γ, η.The dispersity of such systems has been measured in experiments, for example by Walde et al [22].Marques [16] analysed the slow evolution of vesicle size distributions over periods up to 3 years.Thus, given an experimentally determined cluster size distribution for such a system, values for a, b, γ, η can be determined by fitting the data to our partition function.A more detailed theoretically-based model for the polydispersity of sodium dodecyl sulfate micelles has been provided by Hall [12], whose data is consistent with rates of the more general form exp(Cr 1/3 ).
The remaining graphs in Figure 1 are, while perhaps less interesting, still physically relevant.Figures 1(e)-(f) represent the case of classical nucleation without a critical cluster size while Figures 1(g)-(h) illustrate undersaturated systems in which only very small clusters form.Whilst these observations are based on qualitative analysis, Wonczak [27] has observed from quantitative molecular dynamics simulations of argon clusters that fragmentation rates have an exponential size dependence.

Steady-state solutions
In cases where γ > η the equilibrium solution is divergent at large aggregation numbers, so we consider the broader family of steady-state solutions as potential attractors for the large-time limit of the system.Steady-state solutions have the form J r = J independent of size and thus in which the steady-state flux J is given by The large cluster size asymptotics for the steady-state solution (3.6) yield In the case of 0 < γ − η 1 and θe γ = 1 + o(1), the above summations (3.6)-(3.7)can be approximated by their continuum limits, yielding with steady-state flux J given by where erfc is the complementary error function (see Abramowitz & Stegun [1] for details).Instead of assuming that the system will converge to the equilibrium solution, let us now assume that the system converges to a more general steady-state, that is where J r = J for some given value J.We then use (3.7) to determine which value of J is approached.The equilibrium solution is a possible solution, since it corresponds to J = 0, but this will only be appropriate in cases where the sum on the right-hand side of (3.7) diverges; this occurs when η ≥ γ, since in the case η = γ we also get J = 0 since γ < 0.
There are three generic cases depending, in the first case, on the relative size of γ to η and if η < γ then on the relative size of γ to log(a/b): (i) if η > γ then the equilibrium solution decays to zero in the large-r limit, with (3.5) showing that the decay is extremely rapid.Thus the equilibrium solution will be approached in the large-time limit.
(ii) if η < γ then the equilibrium solution diverges in large r limit, and grows more rapidly than any exponential; whereas the steady-state solution (3.7)only grows exponentally.Thus in the large-time limit it is the steady-state solution which is approached.
(iii) η = γ.In this case the equilibrium solution behaves like an exponential in the large r limit (c eq r = θ r−1 c 1 ), as does the steady-state solution (c sss r = (J/ac 1 )e −γr ).Thus the large-time behaviour depends on which solution decays the most rapidly at large r.We have two subcases: (iii.a) η = γ and θ ≤ e −γ .In this case the steady-state solution grows faster than the equilibrium solution at large r, thus the system approaches the equilibrium solution in the large-time limit.
(iii.b) η = γ and θ > e −γ .In this case the equilibrium solution grows faster at large r than the steady-state, so the system will approach the steady-state solution.
These cases are summarised in Figure 2.

Example application: the kinetics of self-reproducing micelles
As an example where normal size distributions such as those in Figure 1(d) may be applicable we consider experiments performed by Luisi's group [2] which have possible relevance to the origin on life on Earth.In a well-stirred vessel, water-immiscible ethyl caprylate is hydrolysed under alkaline conditions to form the water soluble salt of the caprylate surfactant monomer.As the concentration of monomers increases, a well defined point in time is reached after which the total caprylate anion concentration increases sharply.This phenomenon is caused by the crosscatalytic effect of micelle formation: the ester still present is solubilised inside the micelles dispersed in the aqueous phase, wherein it is more rapidly hydrolysed by aqueous hydroxide ions.
Denoting the ethyl caprylate by E, and micelles of r caprylate monomers by C r , we thus have the reaction where c r (t) = [C r ] is the concentration of micelles of size r, ε is the small background rate of ethyl caprylate hydrolysis and the rates k r denote the faster rate due to the cross-catalytic rôle of micelles.The formation of micelles occurs by the standard Becker-Döring cluster formation mechanism, C r + C 1 C r+1 through which clusters grow and fragment by the addition or loss of a single monomer at a time.
The microscopic model is constructed by writing e(t) for the concentration of ethyl caprylate at time t to obtain the kinetic equations ) This model was studied in an earlier paper [10] where a coarse-grained reduction procedure for the Becker-Döring system of equations was outlined.The system corresponds closely to the extended Becker-Döring system described at the start of Section 2, in that a precursor chemical decays to form the monomer.In the case of ethyl caprylate the decay rate depends on the amount of micellar material present in the system and so couples the cluster-formation part of the system to the rate of production of monomer, leading to complex kinetics.In order for the equilibrium solution to have the correct shape, that is, (similar to Figure 1(d).but with the vast majority of the body of the bell-curve in the region r ≥ 1), we choose the rate coefficients to have the form (3.1).For example, we choose a r = a, independent of size, which corresponds to γ = 0, and then determine θ = ac 1 /b and η by requiring the peak of the distribution to occur at r c = 62 -the size of the most probable micelle.This is related to the parameters η, θ = ac 1 /b by If c 1 c rc then the mass of the system is given by and together with r c this determines η and θ.
If we assume that the cross-catalytic rate of hydrolysis of ethyl caprylate by the micelles is a surface phenomenon and then the rate k r is proportional to the surface area of micellar material present in the system and we have k r = kr.The constants k and ε have to be fitted to kinetic data.The rate of loss of ethyl caprylate is then ε + k ∞ r=2 rc r which can be approximated by ε + k( − e).This leads to the equation ė = −e(ε + k − ke) which, in the limit ε 1 has the solution This has the form of a smoothed step function taking the value e(t) = for t < t c := (−1/k ) log(ε/k ) 1, and e(t) = 0 for t > t c , with the width of the transition region being O(1) t c .For t > t c , the system is hence well-approximated by the mass-conserving form of the Becker-Döring equations (2.5) with (2.10).

Large-time solutions: asymptotic and numerical results
We now return to the constant monomer concentration formulation of the Becker-Döring equations.The convergence to the equilibrium or steady-state solution at large aggregation numbers will depend on the relative sizes of γ and η in a way in which we shall now derive.We retain the notation of θ = ac 1 /b for the ratio of aggregation to fragmentation rates.The solution which is approached in the large-time asymptotic regime has been defined above; it is either the equilibrium solution c eq r or the steadystate c sss r , and we determine the kinetics by making the substitution c r (t) = c eq r ψ(r, t) or c r (t) = c sss r ψ(r, t) and developing a partial differential equation for ψ(r, t).As γ varies over γ ≤ 0 and η over the whole of R, there are six main regimes which require separate analysis, summarised in Figure 2. The first three are cases which have been analysed previously, for these we simply quote the large-time asymptotic results from [26]; the other three are new and exhibit novel similarity solutions which we now derive.
Regions of (γ, η) space for which large-time asymptotic solutions are given in Section 4.

Case I: γ = 0 = η
The equilibrium solution is then c eq r = θ r−1 c 1 which is divergent at large r for θ > 1; thus in this case the steady-state solution c sss r = c 1 is the attractor in the large time limit with J = c 1 (ac 1 − b).For θ < 1, the large time limit, from [26], is ).For θ > 1, the steady-state solution is approached by In this case, at large aggregation numbers, the fragmentation rate dominates the aggregation rate, and grows so fast with increasing r that the convergence to the equilibrium state is uniform as Case F(v) of [13] (a r = ar p , b r+1 = br q with q > p).

Case IV: η = γ < 0
As with Case I, the partition function (3.3) simplifies to Q r = (a/b) r−1 and the equilibrium solution is c eq r = θ r−1 c 1 ; however, the range of validity (θ ≤ 1 in Case I) is not the same here.The steady-state solution has the form Thus if θ r decays more quickly than e −γr (that is, θ ≤ e −γ ), the equilibrium solution c eq r = θ r−1 c 1 will be approached (corresponding to J = 0), but if e −γr is the more rapidly decaying function (that is when θ > e −γ ), then a steady-state will be approached, and this has flux so that is appropriate when θ ≥ e −γ .When θ = e −γ , c eq r and c sss r coincide, and so there is a continuous transition from equilibrium to steady-state.We now consider the large time kinetics.

4.4.1.
Case IV(a): η = γ < 0, θ < e −γ .The system converges to the equilibrium solution and the subtitution ψ(r, t) = c r (t)/θ r−1 c 1 yields for which, following the transformation y = e −γr , the similarity form ψ = g(y/t) becomes apparent.The ordinary differential equation generated by this similarity reduction is Table 1.The first few functions of the form which has the solution Let n be the power of x within the above integrands.Then the integral can be evaluated in terms of elementary functions for the special cases of 2n ∈ Z.This corresponds to θ being related to γ by θ = e −γ (2 + γ − 2n)/(2 − γ + 2n).If n is a positive integer (or zero) then the integral is a polynomial of similar degree multiplied by exp(2x/bγ(θ + e −γ )).If the power of x is of the form m − 1 2 for some m ∈ N then the integral has the form of an erfc function summed with the product of an exponential and a polynomial in √ χ.This pattern is summarised in Table 1, and illustrated in Figure 3 for the first few values of m and n.For the integer values of n the term in brackets is the first few terms of the Taylor series of e −cχ expanded about zero; for the intermediate values of n, the term in brackets is the Taylor expansion of 1 2 e −cχ erf( √ −cχ) √ (−π/x) about zero.
The large-time asymptotic solution is then Figure 4 illustrates ψ in the case θ = e −γ (1 + γ)/(1 − γ) where we have This appears to have the form of a travelling wave similar to Case I(a) where the form of the solution is a wavefront with position s(t) and width which scales with w(t), so that its shape is determined by the variable z = (r − s(t))/w(t).However, the solution (4.12) has a different form and is actually described by the single similarity variable χ = e −γr /t.The position of the midpoint (ψ = 1 2 ) is given by r = − 1 γ log(χ c t) where χ c is the critical value for which g(χ c ) = , the earlier ones having more rapid decay at small χ. Figure 5 shows the results of a numerical simulation of the system using Matlab 6 [17], illustrating the convergence of the solution to the similarity solution.The special solvers ode23s and ode15s were used to solve the stiff systems of ordinary differential equations and reach the large times displayed.The dashed lines in Figure 5 were obtained from a numerical solution of (2.5) with rates given by (3.1) with η = γ < 0 and θ < e −γ .4.4.2.Case IV(b): η = γ < 0, θ = e −γ .In this case (4.6) implies J = 0 so the equilibrium solution is approached in the large-time limit.Equation (4.8) simplifies to ψ t = ac 1 e γr ψ rr which again has a similarity solution with variable χ = y/t and y = e −γr .The function g(χ) defined in Case IV(a) diverges at small χ in this case.To accommodate this, we require a modified form of similarity solution; we assume ψ(r, t) = h(t)g(χ), with h(t) = 1/(log t) p with p > 0 and χ = y/t(log t) q ).Balancing the leading-order terms as t → ∞, we find q = 0 and the ordinary differential equation for g(χ), being given by ac 1 γ 2 (χg (χ) + g (χ)) + χg (χ) = 0, (4.13) is independent of p.This equation is solved by g(χ) = AE 1 (χ/ac 1 γ 2 ), where E 1 (•) is the exponential integral [1], this solution satisfies the boundary condition ψ → 0 as r → ∞, which implies g(χ) → 0 as χ → ∞.Imposing the condition ψ(1, t) = 1, implies h(t)g(e −γ /t) → 1 as t → ∞, thus we have p = 1 and A = 1.The large-time asymptotic solution is then Figure 6 shows numerical results for a simulation of the system (2.5) with rates given by (3.1) where η = γ and θ = e −γ .The presence of log terms in (4.14) means the solution converges to the similarity solution slowly.The fact that the numerically computed solutions converge to the predicted similarity solutions in both this case and in the above (Case IV(a)) suggests that the similarity solutions are attractors for the dynamics in the large-time dynamics limit of (2.5) with (3.1) and η = γ.4.4.3.Case IV(c): η = γ < 0, θ > e −γ .Since in this case the solution will approach the steady-state (4.7), we make the substitution ψ(r, t) = c r (t)/θc 1 e γ(1−r) .We then find which has a similarity reduction of the form ψ = g(y/t) with y = e −γr of the form with solution The large-time asymptotic solution is then c r (t) ∼ e γ(1−r) c 1 g(e −γr /t) for r = O(log t) as t → ∞. )) the integral is expressible in terms of elementary functions and has the same form as those displayed in Table 1.
In this case fragmentation dominates at large aggregation numbers, causing the system to evolve towards the equilibrium solution, and fragmentation also dominates the kinetics in the large-time limit.The system evolves according to the partial differential equation where ψ = c r (t)/c eq r .When η < 0, a similarity solution of this is available, as described below, but for η > 0 the convergence is uniform as in Case F(v) of [13].
For η < 0, (4.19) can be solved by a similarity reduction; as in Case IV(a) the variable y = e −ηr and χ = y/t with ψ = g(χ) yields with solution

.21)
Thus the full solution at large times is When η = −2/(n + 1) or η = −2(n + 1 2 ), the power of x in the integrand of (4.21) is a multiple of one half and the exact solutions of Table 1 are applicable.Such a solution is illustrated in Figure 7, the parameters correspond to γ < η with rc = 46, and hence to the equilibrium solution illustrated in Fig 1(c)-(d).In this case the system evolves to the steady-state solution (3.6), which for large aggregation numbers r has the asymptotic behaviour c sss r ∼ (J/ac 1 )e −γr and is thus divergent (since γ < 0).At large times the kinetics can be determined by considering ψ(r, t) = c r (t)/c sss r which satisfies ∂ψ ∂t = and the transformations χ = y/t and ψ = g(χ) yield and hence

Summary
In this section we have derived the large-time asymptotic solutions for the system of Becker-Döring equations (2.5) with aggregation and fragmentation rates given by (3.1) -that is where the rates depend on the cluster size in an exponential fashion.We have found six types of solution, depending on the size of the parameters γ and η, the regions of validity of each parameter combination is illustrated in Figure 2. In cases I, II and III, either the aggregation rate, or the fragmentation rate or both rates are size-independent and such cases have been considered before, so we simply quoted the results from [13].In the remaining cases novel similarity solutions were obtained.In Case V where γ < 0 and η > γ the system approaches equilibrium according to c r (t)/c eq r being a function of the similarity variable χ = e −ηr /t.In Case VI, where γ < 0 and η < γ the large-time kinetics of the system is again controlled by a similarity variable, this time χ = e −γr /t and the system approaches a steady-state solution according to c r (t)/c sss r = g(χ).On the borderline between these two cases, there is Case IV, where γ = η < 0 and here the size of θ = ac 1 /b plays a role in determining the large-time kinetics, there being three subcases.For θ < e −γ , the system converges to equilibrium with c r (t)/c eq r being a function of χ = e −γr /t, and forθ > e −γ the system converges to steady-state according to c r (t) = c sss r g(e −γr /t).The shape of the similarity functions g(χ) are given by integrals, which in a few cases can be explicitly evaluated in terms of exponentials and the error function.At the border between these two subcases, where η = γ < 0 and θ = e −γ there is a special case where c r (t) ∼ (1/ log t)E 1 (e −γr /ac 1 γ 2 t).

Discussion
We have shown that there are families of Becker-Döring models with exponentially sizedependent rate-coefficients that are interesting for a number of reasons.Firstly, the models describe equilibrium solutions which correspond to many physically relevant systems such as micelle and/or vesicle size distributions as well as those having the form of classical nucleation free energy profiles, both with and without a nucleation barrier (associated with the presence or absence of a critical cluster size).These are illustrated in Figure 1.The form of rate coefficients analysed here, namely those which vary exponentially with cluster size is currently finding experimental support even in very simple systems: for example, the work of Wonczak et al (2001) suggests that the fragmentation rates of argon clusters decay exponentially with increasing size.
Secondly, when the aggregation and fragmentation rates depend exponentially on the cluster size, then the concentrations in the constant monomer form of the Becker-Döring equations converge to equilibrium or steady-state via one of a large family of self-similar distributions.(The only exceptions to this rule are the borderline cases where either aggregation or fragmentation rates are size-independent).We have illustrated the range of the two-parameter family of similarity solutions open to the system in Figures 3, 4 and 7.The self-similar behaviour is robust and ubiquitous; it occurs for wide ranges of parameter values, not just at specially tuned systems.We have illustrated the selfsimilar behaviour by numerically solving the full system of differential equations and rescaling the data to show that the data curves accumulate onto a single distribution (Figures 5-6).This shows that the behaviour is stable, in that the similarity solutions acts as attractors in the large-time limit.
In the analysis of previous forms of rate coefficients, self-similar distributions have only been found for borderline cases, for example, in Case I when θ = 1, which separates two different generic behaviours, one in θ < 1 and the other for θ > 1.Here the opposite is the case; we find self-similar behaviour in all of Cases IV, V and VI, of Section 4 with only the borderline Cases I, II and III showing more complicated behaviour.
In conclusion we have shown that the Becker-Döring schemes with rate coefficients which depend on size in an exponential fashion exhibit self-similar kinetics in their approach to equilibrium and steady-state solutions.Such models can be used to describe a number of important systems in classical nucleation theory and the self-assembly systems of physicochemical and biological interest.

Figure 1 .
Figure1.Exponential rate coefficients.Graphs of log Q r and Q r against aggregation number r for the cases r c > < 1 and γ > < η.The parameters γ, η determine the sizedependence of the aggregation and fragmentation rates as described by equation (3.1); r denotes the aggregation number and Q r the partition function, which is related to the chemical potential by kT log Q r = µ r , such that Q 1 = 1.The quantity log Q r is plotted in the left-hand column of figures: at constant temperature, log Q r is directly proportional to the free energy (chemical potential) of clusters of size r.Equation (3.4) defines the quantity r c = rc, which denotes the other cluster size of the same chemical potential as the monomer, i.e.where Q rc = 1.

Figure 5 .
Figure 5. Graphs of the numerical results for Case IV(a), ψ = c r (t)/c eq r is plotted against the rescaled aggregation size χ = e −γr /t to illustrate convergence to the similarity solution g(χ) at larger times.The parameters used are γ = η = −0.1, a = 1, c 1 = 1, b = e γ (1 − 2γ)/(1 + 2γ), and curves are plotted at times t = 3 τ with 3 ≤ τ ≤ 8, for system size 1 ≤ r ≤ N = 1100.The dashed lines show the numerical results, and the solid line the expected self-similar solution.

Figure 6 .
Figure 6.Graphs of the numerical results, g = c r (t) log(t)/c eq r is plotted against the rescaled aggregation size χ = e −γr /t to illustrate convergence to the similarity solution g(χ) at larger times in Case IV(b).Here, a = 1, c 1 = 1, b = e γ , γ = η = −0.1, the system size is 1 ≤ r ≤ N = 1100 and results are shown for 1 < t < 10 40 in steps of 10 4 , the dashed lines show the numerically results, and the solid line the expected self-similar solution.