Exact solutions for cluster-growth kinetics with evolving size and shape profiles

In this paper we construct a model for the simultaneous compaction by which clusters are restructured, and growth of clusters by pairwise coagulation. The model has the form of a multi-component aggregation problem in which the components are cluster mass and cluster diameter. Following suitable approximations, exact explicit solutions are derived which may be useful for the verification of simulations of such systems. Numerical simulations are presented to illustrate typical behaviour and to show the accuracy of approximations made in deriving the model. The solutions are then simplified using asymptotic techniques to show the relevant timescales of the kinetic processes and elucidate the shape of the cluster distribution functions at large times. PACS numbers: 64.60.Ht, 82.60.Nh, 05.70.Fh


Introduction
It would be useful to have models of nucleation which describe the differences in both size and shape of growing clusters and yet are simple enough to be solvable analytically.Current exactly solvable models of coagulation only describe cluster masses.As well as elucidating the kinetics of aggregation and compaction, models involving size and shape would be useful in the testing of numerical simulations of systems such as those used by Xiong et al [24,25].An alternative approach which takes explicit account of the separately evolving size and shape of a typical cluster is given by Schild et al [15].Although this is useful, it only follows one cluster, so the method cannot output a distribution of sizes and shapes which can be tested against experimental observations.Typically one expects collision events, which allow aggregation, also to cause compaction of clusters; the long-term effect of this is to transform fractal aggregates (similar those observed in diffusion limited aggregation (DLA)) into more compact clusters.Modelling this by a singlecomponent coagulation process is complicated and requires many assumptions to be made [18].This approach is pursued by Vemury and Pratsinis [20], where such a one-component model is analysed in an attempt to find the self-preserving shape of cluster-size distribution.
Ideally the results described below should be calibrated against experimental data, or data from computer simulations, for example the work on sintering carried out by Akhter et al [2], where computer simulations are compared to experimental data.However, we would not expect exceptionally good agreement from the models solved in this paper, since these have no size or shape dependence in the aggregation kernels.More realistic kernels could be used, and then it would be interesting to compare numerical solutions of such models against experimental data and other simulation techniques.Kostoglou et al [13,14] and di Stasio et al [17] have also worked on modelling simultaneous coagulation and restructuring of clustershape.Another area where two-component aggregation problems naturally arise is that of charged clusters [3,4,19], where clusters are characterized independently by size and charge.
Other work on multi-component coagulation problems includes that of Elvingson and Wall [11,21], who developed a two-component version of the Becker-Döring equations to model the formation of mixed micelles; these are clusters formed from two-species of surfactant molecule.A similar model has been analysed by Wu [23].Multi-component Becker-Döring systems have been used in several models of the kinetics of vesicle formation [6,7,9].However, the situations under consideration in this paper require Smoluchowski [16] aggregation rather than the restricted stepwise growth of Becker-Döring models.An unusual multi-component coagulation which includes Smoluchowski-type aggregation arises in the modelling of river-flow [8], where to make progress on the analysis the system is again approximated by a single-component problem.In a few special cases of multi-component Smoluchowski aggregation, exact solutions are available [22], and the solution constructed in this paper relies on the ideas and methodology presented there.
In section 2 we derive a multi-component model of simultaneous coagulation and compaction which is of Smoluchowski-type.This model is solved in section 3 by means of generating function techniques.A numerical solution is also performed to allow us to analyse some of the errors made in the modelling assumptions.The large-size and large-time asymptotics of the exact solution is carried out in section 4-this allows some simplification of expressions.Finally a discussion of the results is presented in section 5.

Formulation of model
In our model we associate two parameters with each cluster: as in Smoluchowski's model of aggregation [16], we partition the distribution of clusters according to cluster mass; the novel feature of this work is that we also partition clusters according to their maximum diameter.Thus we denote a cluster of mass j and maximum diameter k by C j,k .
We then allow two processes to act on the distribution of cluster sizes: a restructuring of the cluster which transforms a fractal aggregate to make it more compact.This occurs through some geometric rearrangement of the cluster's constituent parts so as to reduce its maximum diameter, as illustrated in figure 1.To each compacting event we assign a transition rate γ j,k .Once a cluster is maximally compact (i.e. its maximum diameter has reached some minimum) this process will be assumed to have no further influence on a cluster.If we follow the spherical liquid drop model of a cluster, then the minimum diameter for a cluster composed of j monomers is k c (j ) such that 4  3 π 1 2 k c (j ) 3 = j/σ , where σ is the density of a monomer (since j is a measure of mass, j/σ is a volume).Thus k c (j ) = (6j/πσ ) 1/3 ; if we work in units in which the monomer has unit diameter, we find σ = 6/π ; and hence k c (j ) = j 1/3 .In other applications, where clusters may preferentially form rod-like or disc-like aggregates, some other functional form of k c (j ) may be more appropriate.Whilst we are interested in the full range of cluster sizes 1 j < ∞, the range of maximum diameters k is restricted to k c (j ) k j .This region of (j, k) space is illustrated in figure 2.
The second process which occurs is coagulation, by which two clusters combine.Clearly, the masses must simply sum, but clusters may combine in any orientation, so the maximum diameter of the aggregate may be less than the sum of the maximum diameters of the initial clusters.Formally, we have with q potentially taking any value from max{k, s} to k + s.In a mean field model, we should form some weighted average over all possible configurations.However, since we have a mechanism to reduce the cluster's maximum diameter, we take the 'worst' case scenario of the greatest possible value of q, and allow the restructuring mechanism (figure 1) to spread the resulting distribution over a range of diameters smaller than k + s.This mechanism is illustrated in figure 3. Magnetic or electrically charged particles tend to form extremely elongated structures during growth by coagulation (see, for example, Kammler et al [12]).
Taking the two mechanisms illustrated in figures 1 and 3 and applying the law of mass action to derive equations for the concentrations c j,k (t), we obtain dc j,j dt = F j,j − L j,j − γ j,j c j,j ,  where F j,k and L j,k are the rates of formation and loss of clusters C j,k through the usual aggregation and fragmentation processes, that is, (2.4)

Integrable model of coagulation and compaction
Whilst model (2.2) is interesting and can be solved numerically, our aim here is to construct a model which is explicitly and exactly solvable.Hence, we simplify equations (2.2).Firstly we specify the rate coefficients a r,s,j,k and γ j,k ; the simplest aggregation rates to solve are typically size-independent, thus we assume a r,s,j,k = a.In place of the lower limit k = k c (j ), we first approximate k c (j ) by 1 + ε(j − 1) with ε small; however, such models are in general still insoluble; to obtain an integrable system we take the limiting case and put ε = 0.This is equivalent to defining k c (j ) = 1.In section 3.3 we solve the systems numerically and analyse the differences between systems with k c (j ) = 1 and k c (j ) = j 1/3 .For the compaction rate γ j,k we assume γ j,k = γ (k − 1) for some constant γ , since this automatically becomes zero on the line k = 1, simplifying the mathematical formulation of the problem.Physically, it also has the advantage of assigning a high rate of diameter reduction to very 'wispy' aggregates (whose maximum diameter is close to their mass) and low rates of diameter reduction to clusters which are almost maximally compact.Thus as well as being mathematically convenient, we believe this to have good physical justification.For simplicity we do not include any size-dependence (j ) in γ j,k : the rate at which the maximum diameter reduces depends only on the maximum diameter.Combining these assumptions for γ j,k and a j,k,r,s , we obtain note that on the line k = 1, which represents the maximally compact clusters, the last term automatically vanishes, since no further compaction of these clusters can occur.The approximation of k c (j ) by unity increases the region of (j, k) space accessible to the model, as illustrated in figure 4. Also, note that provided that c j,k (t) = 0, for k > j, is satisfied at t = 0 then this condition is satisfied by the distribution at all later times.Thus we need to make no explicit specification of the condition k j in (2.5), or write out a special equation valid on k = j , since c j,j +1 = 0 automatically causes the penultimate term of (2.5) to vanish on k = j + 1.At t = 0, we assume the system is completely in monomeric form, that is, the initial data are c j,k = 0, for all j, k with the exception of c 1,1 = .Illustration of region of size (mass j ) and shape (maximum diameter k) parameter space included in our explicitly solvable model system of section 2.2.This corresponds to the case ε = 0 of the linear approximation to the lower limit illustrated in figure 2.

Solution of model
We will aim to solve the system using the generating function approach of Davies et al [10]; hence we introduce transform variables x, y and define and we also make use of an alternative generating function which is related to C(x, y, t) by G(x, y, t) = C(x, y, t) e x+y .Associated with these functions, we define functions which represent the total number of clusters in the system The initial conditions for C are C(x, y, 0) = e −x−y , whilst for G they take on the simpler form of G(x, y, 0) = ; both of these imply C 0 (0) = G 0 (0) = .
The equation for the generating function C(x, y, t) is The associated equation for C 0 (t) is Ċ0 = − 1 2 aC 2 0 which, when the initial condition C 0 (0) = is imposed, has the solution Solving by characteristics, with initial data on s = 0 parameterized by ξ, η and Expanding G(x, y, t) as a power series in both e −x and e −y , we find the full explicit solution for each individual concentration where This is the exact explicit solution of the problem originally posed in (2.5), with initial data of c j,k = 0 for all j, k except for c 1,1 = .Although our aim was to construct such a solution, and it will be useful for verifying numerical solutions of such problems, it is not clear exactly what behaviour is described by this function.Hence in the next section we will form approximations of it to show the kinetic phenomena it describes.

Analysis of moments
Using the generating function (3.8), we now find properties of the distribution, such as the behaviour of the first few moments.We define the joint moments by The number of clusters is given by M 0,0 = G(0, 0, t); however, higher moments are given by more complex formulae.Since C = e −x−y G = ∞ j =1 ∞ k=1 c j,k e −jx−ky , we have for example, M 0,1 = G(0, 0, t) − G x (0, 0, t) and M 2,0 = G xx (0, 0, t) − 2G x (0, 0, t) + G(0, 0, t).In particular, we have ) where we have defined the time-dependent quantity E(t) by (3.11).
From the moments (3.14)-(3.16), it is possible to derive quantities of macroscopic interest.For example, the average cluster size, J , can be derived in several ways: .17) In systems which do not undergo any sort of gelation behaviour, all these definitions should give rise to broadly the same kinetic behaviour.Note that J 1 and J 2 are independent of γ , as one would expect, since compaction does not alter cluster mass and neither does it affect the subsequent rate of cluster coagulation.However, J 3 depends on γ -this definition of average cluster size showing some influence of the restructuring history of clusters.For a wide range of parameters γ, a and most times t, the value of J 3 lies between J 1 and J 2 : at small times J 3 is close to J 2 , at larger times J 3 approaches J 1 .For small γ /a the crossover of J 3 from J 2 to J 1 occurs at large times, and for large γ /a the crossover occurs at small times.

Fractal dimension
In our definitions, the volume of a cluster C j,k scales with its aggregation number, thus V ∼ j , and the diameter scales with L ∼ k.For fractal clusters, the dimension D is defined by Using definitions (3.17)-(3.21),nine different fractal dimensions can be constructed, At small times, we expect the growing clusters to be linear in geometry, thus to have dimension close to unity.However, of the nine definitions two give rise to dimensions of two (D 2,1 , D 3,1 ) and two more to dimensions of one half (D 1,2 , D 1,3 ).This leaves five definitions of dimensions, which are plotted in figure 5. From this we see that two give quite low estimates, and one of these gives dimensions below unity, which we discount as unphysical.In the left-hand graph, where γ = 0.01, we expect compaction to occur on the timescale t = O(1/γ ); thus the second lowest curve also gives a compaction timescale which is unexpectedly long.The upper three curves all give qualitatively similar results.A similar outcome is observed when γ = 1, suggesting that the definitions D 1,1 , D 2,2 and D 3,2 should be preferred over the others.The upper curves in figure 5 exaggerate the compaction effect, as fractal dimensions greater than three are only possible because of the approximation we make when relaxing the constraint k (6j/πσ ) 1/3 to k 1.
A more accurate calculation of the fractal dimension may be achieved by noting that the dimension of a particular cluster c j,k is log k/ log j .The average fractal dimension of the whole population is thus given by an average of the form for some constants p, q, r, s.Unfortunately formulae such as (3.23) cannot be explicitly be evaluated given the form of our solution (3.10).

Numerical solution
In figure 5 we observe several of the curves rising rapidly to dimensions above three.The reason for this is that in section 2.2, when deriving a set of equations which are explicitly integrable, we replaced k c = j 1/3 by k c = 1.This significantly alters the calculation of the dimension of the more compact clusters.To assess the implications of this approximation, we have used Matlab to solve the system of ordinary differential equations (2.2) with γ j,k = γ (k − k c (j )) and a r,s,j,k = a in both the cases k c (j ) = 1 (the integrable case) and k c (j ) = j 1/3 .The outputs were used to calculate the average cluster sizes J p and diameters K q and the fractal dimensions D p,q (1 p, q 3).A sample of the results are shown in figure 6, where J 2 , K 2 and D 2,2 are plotted against time.
Whereas excellent agreement is seen for the average cluster size J 2 , differences are clearly visible in the maximum diameter, K 2 .For the system with k c = 1, K 2 has a more rapid decay than the system with k c = j 1/3 .The differences are more pronounced when the fractal dimension D 2,2 is calculated (right-hand graph in figure 6).Some of the difference can be corrected for a posteriori, as we shall now show.
The range of diameters 1 k j used in our integrable model can be mapped onto the range j 1/3 k j in the more realistic model through the scaling which is affine linear in k.The definitions of the moments can then be modified to with corresponding new definitions for the average sizes Ĵ q , Kq and Dp,q .In figure 6 it can be seen that these modified quantities lie closer to the corresponding quantities for the system with k c = j 1/3 than the system with k c = 1.

Asymptotics
In this section we return to the special case k c (j ) = 1 for which the explicit solution is available and we aim to describe in simpler terms the kinetics it describes.Using 5.1.7 of Abramowitz and Stegun [1], we rewrite the solution (3.10)-(3.11)as Although it is useful to have the exact explicit solution (3.10), the expression is too complex to give an intuitive feel of the dynamics it describes.In this section we make use of asymptotic approximations to give simpler functional forms of the solution at large times and for larger clusters (j, k, t 1).This procedure will also highlight any potential similarity solutions which may be approached.
It is well known that many aggregation phenomena exhibit self-similar scaling behaviour at large times and large cluster sizes.To show the connection with already solved models we write Using solution (3.10), we recover the classical solution S j (t) = 4 (2+a t) 2 a t 2+a t j −1 for the additive kernel.This has the large-time asymptotic form This result implies that the typical cluster size scales linearly with time; hence we introduce the scaled size variable η = j/t.From this we note that the aggregation number (mass) of the cluster does not depend on the cluster's diameter (k) or on the compaction rate (γ ).In more realistic aggregation kernels the aggregation rate (a j,k,r,s ) would depend on both the mass and diameter of the cluster and this extra effect could lead to some correlation between cluster mass and compaction rate (γ ).
There are two obvious special cases of (4.1)-(4.2) which may lead to particularly simple forms, and we examine these first: they are the cases where aggregation and compaction occur on vastly different timescales.The crucial asymptotic formulae we make use of are the power series obtained by expanding Ei(x) about x = 0, namely

Rapid compaction and slow coagulation (γ a )
This is the less interesting of the two special cases: clearly whenever any cluster c j,k is formed, it will be compacted down to c j,1 over a very short timescale.Over a longer timescale the distribution of cluster sizes will evolve, being dominated by c j,1 .

Rapid compaction, faster timescale
This expected behaviour is confirmed by solution (4.2).We put a ∼ O( 1) and γ 1, then formally define the initial rapid timescale by τ = γ t; however, since the initial conditions are fully compact there is no dynamics over this timescale.For completeness with later calculations, we note that the asymptotic form of E(τ ) = E(t) over this timescale is E(τ ) ∼ 1 − e −τ − a τ/2γ ; thus E(τ ) grows from zero towards a maximum of unity where it saturates.The small correction term suggests that over the longer timescale E(t) will start to decline.

Rapid compaction, slower timescale
Over the longer timescale (t = O(1)) each term in the quantity E(t) can be expanded giving confirming our earlier indications that E(t) declines over this slower timescale.The decay of E(t) at large times (t 1) is algebraic, with E(t) ∼ 2/a t as t → ∞.In this limit (4.2) can be approximated by (4.8) Thus we see that each extra power of k reduces the concentration by a factor of 1 2 γ a t 2 , which is extremely large since we are considering both γ 1 and the large t limit.The interesting asymptotic scaling of size with time is j ∼ t, for which we define η = j/t and hence obtain and Thus we see that the concentrations c j,1 (t) do indeed dominate the system, as expected, and exhibit self-similar growth in cluster size.

Slow compaction and rapid aggregation (γ a )
For the more interesting case we retain a ∼ O (1), and assume γ 1; thus the two timescales are that of aggregation t = O (1), and the slower one being T = γ t = O(1), equivalent to t = O(1/γ ) 1.

Slow compaction, faster timescale
On the faster timescale t = O(1), we find that the quantity E(t) can be approximated by and so is uniformly small (provided that log t 1/γ , which is certainly satisfied, since the new timescale introduced in section 4.6 is t ∼ 1/γ ).
Expanding the exact solution (3.10), for t ∼ O (1) with γ 1, we find (4.12) Since γ 1, the concentrations c j,k (t) are only of significant size (O( 1)) along the line j = k, where, at large times, we observe Some spreading of mass into the region k < j starts to occur at larger values of j, k and at larger times; however, spreading of the most numerous clusters from the line k = j to the compact state k = 1 does not occur until the longer timescale t = O(1/γ ) is entered.

Slow compaction, slower timescale
To analyse the slower (t 1 for which T = γ t = O(1)) timescale, we firstly consider the form of E(t) = E(T ) (4.1).For small T we have linear growth with E(T ) ∼ T and for large T we find E(T ) is small and decaying algebraically with E(T ) ∼ 2γ /a T .Numerical evaluations of the function show that there is a single maximum between these two limits, and using asymptotic analysis based on γ a we find the location (T c ) and height (E c ) of the maximum are given by 1, The quantity E(T ) has the same form in the limit γ a as E(τ ) has in the limit γ a , namely is zero at T = 0, rises to a maximum and then decays.When γ a the value of E(τ ) at the maximum is close to unity, whereas for γ 1, E(T ) is small even at the maximum T = T c .
To investigate the asymptotic behaviour of the solution (4.2) at large times and large cluster sizes, we introduce the scalings j = J /γ and k = θj = θJ /γ, θ = k/j (with θ ∈ (0, 1)) being the relative compactness of the cluster C j,k .This leads to where and θ c is the position of the maximum of H (θ, T ).The dominant contribution to the shape of cluster distribution function in (j, k) space is due to the term H (θ, T ).To simplify this term we form the Taylor series of H (θ, T ) around its maximum in the manner of Laplace's method [5].Solving H θ = 0, we find the relative compactness (θ ) of the most frequently occurring cluster type, θ = θ c (T ), Note that this is the same θ -value for all cluster sizes J .To quadratic terms, H (θ, T ) is approximated by Combining this with the prefactor given in (4.15), we find Over this timescale we see the transition from most of the mass being focused around the line k = j to the fully compact state where the distribution has its maximum around the line k = 1.
The position of the maximum being given by k = θ c (T )j = j e −γ t .The solution at large cluster sizes and at large times is illustrated in figure 7, where γ /a = 10 −2 , at times t = 4, 30, 70, 140, 400.In the top graphs, the mass can be seen to lie predominantly along the line k = j , which in successive graphs moves to k = 0.74j at t = 30, k = 0.50j at t = 70, k = 0.25j at t = 140, the ratios k/j agreeing well with e −t/100 .At t = 400, we find almost all the system's mass along the line k = 1; consistent with the prediction θ = k/j = e −4 ≈ 0.02.Simultaneous with this change in shape of the clusters, we observe a steady increase in size as the distribution evolves from a large and sharply peaked maximum at j = 1 to much lower concentrations over a broad range of sizes.
At even longer timescales, when t 1/γ , the system is dominated by fully compact clusters, that is, clusters of the form C j,1 .For large j we have the similarity solution c j,1 (t) ∼ 4 e −2j/a t /a 2 t 2 , ( and c j,2 (t) ∼ j e −γ t c j,1 (t), thus c j,2 (t) c j,1 (t).

Compaction and aggregation on similar timescales
Figure 8 shows the case where the rates of aggregation and compaction are similar, that is γ ∼ a .In this case there is no way to simplify the form of E(t).Simply deriving the large-time asymptotic solution of (4.1)-(4.2) will lead to a solution in which the maximally compact clusters dominate all others (c j,k (t) c j,1 (t) for all k 2) and the size-distribution has the self-similar form c j,1 (t) = 4 e −2j/a t /a 2 t 2 .
For large j, k with k = θj and time taken to be O(1), we have

.22)
As with (4.15), the dominant term in the expression for c j,k (t) is e jH (θ,t) and H (θ, t) has a maximum θ c (t) which is independent of j , that is, the transformation from extended (θ = 1) to compact (θ = 0) clusters occurs at the same time for all cluster sizes.This is given by solving H θ (θ, t) = 0 and leads to However, there is no simplification of the expression for H (θ c (t), t) and so no straightforward approximation for c j,k (t) is available.

Discussion
We have formulated a model of cluster growth in which both the size (mass) and shape (maximum diameter) of clusters are explicitly and independently taken into account.The form of the resulting model is that of a multi-component aggregation problem with an additional restructuring process which we have referred to as 'compaction' by which a cluster's maximum diameter is reduced while its mass is left unchanged.
The model is approximated by simplifying the range of maximum diameters allowed from (6j/πσ ) 1/3 k j to 1 k j , where k is the maximum diameter of a cluster of mass j , density σ and hence volume j/σ .Following certain assumptions on the form of the rate coefficients, we have obtained a model which is solvable explicitly using analytical techniques.The resulting solution can be used to check numerical solvers of multi-component systems.
As proposed, the model only allows compaction but generalizations which model processes by which spherical clusters are stretched could easily be incorporated, as illustrated in [14].Due to the fact that only compaction is included in the model, all mass will eventually end up on the curve k = (6j/πσ ) 1/3 and on k = 1 in the approximated model.We have discussed various ways that the average fractal dimension can be calculated; unfortunately the most accurate formula (3.23) does not lead to expressions which can be explicitly evaluated using our asymptotic solution.Instead we have shown that a cruder approximation based on the average cluster size and average diameter can be used to give an indication of the rate of compaction of clusters.
We have used Matlab to analyse the difference between systems where the maximally compact cluster has a maximum diameter of k c = O(j 1/3 ) and the explicitly solvable system where k c = 1.The differences are not large, in the former model the range of maximum diameters allowed at any cluster mass j is j 1/3 k j , which we have approximated by 1 k j .For large masses, j , the relative difference is O(j −2/3 ), which is small.However, by incorporating the range from k = 1 to k = O(j 1/3 ) we are losing some of the geometric information about the allowable structure of clusters.In calculations of the average maximum diameter and fractal dimension, these differences become noticeable, as can be seen in figure 5, where fractal dimensions of 3 and above are rapidly realized.We have illustrated how an a posteriori rescaling of the results by (3.24)-(3.25)can eliminate the majority of the discrepancy in the calculation of the average diameter and fractal dimension (see figure 6 for details).
We have chosen monodisperse initial data and simplified the resulting solution using asymptotics; this enables us to illustrate some of the kinetic features of simultaneous aggregation and compaction.For the combination of aggregation kernel and compaction rates adopted here, the large-cluster size asymptotics are particularly simple: the timescale over which compaction occurs is the same for all cluster sizes.We expect that when more general rate coefficients are employed, large and small clusters may restructure over different timescales.For rapid compaction and slow coagulation, the results are straightforward as clusters are always in their most compact form, and grow in size according to the usual selfsimilar solution.When aggregation is much faster than compaction the kinetics of the solution are more interesting.There is a faster timescale where self-similar growth is seen along the line k = j , that is, linear aggregates form in a self-similar fashion.Over a slower timescale the whole distribution of cluster sizes restructures to the more compact form, whilst continuing to grow in size following the same self-similar rule.For more general rates of growth and shape restructuring one would expect more complex rules, where small and large clusters changed their shape and fractal dimension over differing timescales.

Figure 3 .
Figure 3. Illustration of aggregation of clusters of arbitrary sizes (r and j ) and arbitrary shapes, this being described by each cluster's maximum diameter (s and k).

Figure 4 .
Figure 4. Illustration of region of size (mass j ) and shape (maximum diameter k) parameter space included in our explicitly solvable model system of section 2.2.This corresponds to the case ε = 0 of the linear approximation to the lower limit illustrated in figure2.

Figure 7 .
Figure 7.Left-hand column: plots of c j,k (t)/c 1,1 (t) from (4.2) against j and k.Right-hand column: plots of log c j,k (t) in black and in grey, the approximation (4.19).Descending in sequence, the plots illustrate the shape of the distribution at times t = 4, 30, 70, 140, 400 for the parameter values γ = 0.01, a = 1, = 1.