Crossover in coarsening rates for the monopole approximation of the Mullins–Sekerka model with kinetic drag

The Mullins–Sekerka sharp-interface model for phase transitions interpolates between attachment-limited and diffusion-limited kinetics if kinetic drag is included in the Gibbs–Thomson interface condition. Heuristics suggest that the typical length-scale of patterns may exhibit a crossover in coarsening rate from l(t) ˜ t1/2 at short times to l(t) ˜ t1/3 at long times. We establish rigorous, universal one-sided bounds on energy decay that partially justify this understanding in the monopole approximation and in the associated Lifshitz–Slyozov–Wagner mean-field model. Numerical simulations for the Lifshitz–Slyozov–Wagner model illustrate the crossover behaviour. The proofs are based on a method for estimating coarsening rates introduced by Kohn and Otto, and make use of a gradient-flow structure that the monopole approximation inherits from the Mullins–Sekerka model by restricting particle geometry to spheres.


Introduction
Phase separation by spinodal decomposition or heterogeneous nucleation produces complicated patterns whose evolution is driven by dissipation of surface energy while the volumes of the separate phases are preserved. The spatial structure is observed to coarsen in ways roughly characterized by a typical length-scale l(t) that grows as a power law t α for some exponent α > 0. It is a challenge to understand this kind of transient dynamic behaviour in spatially extended systems with complex morphology, and to explain how different physical mechanisms yield different values of the coarsening rate exponent α.
We study a variant of the Mullins-Sekerka model, a standard sharp-interface model that describes the late stages of phase-separation dynamics in a two-phase mixture. This particular variant alters interfacial motions by including a kinetic drag term that appears in the literature on rapid solidification from a melt. Heuristic arguments suggest that a phenomenon of temporal crossover occurs in this system. The length-scale should grow proportional to t 1/2 in the (early time) regime of 'attachment-limited' kinetics, and proportional to t 1/3 in the (late-time) 'diffusionlimited' regime. Our main result is to establish rigorous bounds on energy decay that partially justify this understanding in the so-called monopole approximation.
The monopole approximation involves constraining the geometry of the minority phase to be a finite system of spherical particles in R 3 . This is a reasonable approximation when the volume fraction of the minority phase is small and its particles are well separated (see [1,2] for a rigorous treatment). If the ith particle has radius R i = R i (t) > 0 and (fixed) centre x i , the (non-dimensionalized) interface evolution law reduces to a system of ordinary differential equations (ODEs) in which the time derivativesṘ i are determined from the equations Here, β > 0 is the kinetic drag coefficient and θ(t) is independent of i and is determined so that total particle volume is conserved in time: This system applies until such a time as some particle's size vanishes. Then the system is restarted with the remaining particles after relabelling. The dynamics of the monopole model remain complex due to the arbitrary arrangement of particles in space. Presumably, spatial correlations between particles develop in a time-dependent fashion and affect coarsening dynamics in ways that are not well understood. Our analysis employs the method introduced by Kohn and Otto [11], who treated two variants of the Cahn-Hilliard equation. The outcome is to establish rigorous power-law bounds for decay of a normalized energy E(t) for the system, in a time-averaged sense. These bounds are (i) universal, in that they apply to every solution; (ii) one-sided, in that slower coarsening is possible (unstable equilibria may exist, for example), but faster coarsening is impossible; (iii) independent of system size and spatial complexity; and (iv) independent of statistical assumptions about the system.
In the present problem, energy is proportional to surface area, and we let This quantity is 1 3 of the ratio of total surface area to total volume of particles, and it scales as inverse length. Our main results establish two different kinds of lower bounds on the decay of E that remain valid independent of initial conditions, as long as particles do not collide. The time of validity of these bounds involves two different dual quantities, L 1 (t) and L 2 (t), which scale as length. (These quantities are defined in § 4; see (4.3) and (4.13).) Theorem 1.1.
(i) Given 1 < p < 2, there exist positive constants c 1 andĉ 1 depending only on p such that, for any solution The meaning of these estimates will be discussed later in more detail, but we note here that they are consistent with the following heuristic expectations. Suppose that after non-dimensionalization R(t) is a typical length-scale (particle radius); we expect Then, as long as R(t) β, typically β + R i ≈ β in (1.1), and formal scaling analysis leads us to expect Eventually R(t) β, whence typically β + R i ≈ R i and formal analysis suggests R(t) ∼ t 1/3 . Crossover should occur when (t/β) 1/2 ≈ β ≈ t 1/3 , or t ≈ β 3 . Within a constant factor independent of β, the lower bound provided by (1.4) is better than that provided by (1.5) for T < β 3 , and vice versa for T > β 3 .
The Kohn-Otto method has been applied successfully to treat a considerable number of other models (see [4,[7][8][9][12][13][14]20,21] for examples and further developments). One reason for interest in the present results is that they relate to the case of coarsening by pure attachment-limited kinetics, which has resisted a complete analysis so far. In this case, the interface motion law states that the normal velocity v is proportional to the deviation between mean curvature κ and its averageκ. With the sign convention we will use, and in appropriate units, this means v =κ − κ.
(1.6) As Otto has pointed out (personal communication), for this law there is no lower bound on energy decay that is completely independent of solution morphology. A system of identical particles and 'holes' has average mean curvature zero (κ = 0), and no universal bound is possible for coarsening by pure mean curvature motion. Yet there are some positive results for (1.6) that establish bounds like (1.4), consistent with length-scale growth proportional to t 1/2 as heuristics suggest. Such bounds are valid when all particles are spheres, for example (with no 'holes'). This follows from the result in [7, theorem 4.3] for systems of either finitely or infinitely many spheres; a simple treatment is shown in § 5. The present paper extends this treatment to study coarsening behaviour in the system (1.1) when crossover to diffusion-limited kinetics comes into play. (In two spatial dimensions, Dai [6] has recently established t 1/2 bounds when all particles are convex regions with bounded eccentricity in a certain sense.) The remainder of the paper is structured as follows. In the next section we describe the Mullins-Sekerka model with kinetic drag and indicate its connection with the law of motion by deviation of mean curvature in the limit as β → ∞. In § 3 we show that the monopole model arises from the gradient-flow structure of the Mullins-Sekerka model by geometrically restricting particles to be spheres. The inherited gradient-flow structure plays a key role in establishing the rigorous coarsening bounds in § 4.
Section 5 deals with the simpler Lifshitz-Slyozov-Wagner (LSW) mean field model with kinetic drag, given bẏ This comes from (1.1) by ignoring the cross-interaction terms involving |x i − x j |. It is a model for very dilute systems in which the distance between particles becomes large compared with their radii. The ODE (1.7) remains weakly coupled by the mean field θ(t), which is determined by conservation of total volume i R 3 i as before.
The coarsening-rate estimates (1.4) and (1.5) are easy to establish for this LSW model (see § 5). Our main reason for considering (1.7) regards its tractability for numerical simulation. In § 5.2 we will present numerical simulations for this model that exhibit crossover of length-scale coarsening from (t/β) 1/2 to t 1/3 . It would be desirable to compute solutions of the monopole model (1.1) to see whether such crossover behaviour occurs as predicted heuristically. However, observing such behaviour by direct simulation of (1.1) appears to be quite difficult due to the large numbers of particles needed for good statistics, the global coupling of the equations and the great length of time anticipated for power-law coarsening to increase length-scales sufficiently (at least two orders of magnitude).

Mullins-Sekerka model with kinetic drag
The model that we start from is formulated as follows for a two-phase mixture in a domain Ω ⊂ R n . Let G t ⊂ Ω denote the region occupied by one of the phases at time t, and let Γ t = ∂G t denote the interface separating the two phases. The system in non-dimensional form is then Here u is a normalized chemical potential or concentration difference, n is the outward unit normal to the interface Γ t , v is the normal velocity of the interface Γ t , κ is the mean curvature of Γ t (taken to be positive when G t is convex) and β > 0 is a constant kinetic drag coefficient. The jump of the normal derivative of u across Γ t is [n · ∇u] = n · ∇u + − n · ∇u − , where u + and u − are the limits of u on Γ t coming along n from the outside and inside, respectively. One should also apply an appropriate boundary condition on ∂Ω, for example, the no-flux condition n · ∇u = 0 on ∂Ω.
For this model, the total volume of G t is conserved in time.
The standard Mullins-Sekerka system has β = 0. Including kinetic drag in the Gibbs-Thomson relation (2.2) provides a formal interpolation between a (faster) regime in which the main rate-limiting mechanism of domain growth is the attachment of monomers at the boundary, and a (slower) regime in which diffusion between different components (particles) becomes dominant. We can illustrate the effect of β by the following heuristic scaling argument. We rescale the chemical potential u, length x and time t in terms of typical scales U , X and T for these quantities by Under this scaling we find that (2.2) becomes (after dropping tildes) Presuming u, κ and v are now O(1) quantities, it is reasonable to take U = 1/X. Then we can distinguish two limiting cases: X β and X β.
Here we may take T = βX 2 and find X 2 /U T = X/β ≈ 0. Then [n · ∇u] ≈ 0, where the harmonic function u should be approximately a constant in space, equal to the average mean curvatureκ, since total volume of G t is conserved. That is, the dynamics are given approximately by mean-curvature fluctuation, as in (1.6). This limiting model is then invariant under rescaling with X = (T/β) 1/2 , which is the length/time-scale relation one then expects heuristically for 'statistically self-similar' coarsening.
Here it is reasonable to take T = X 3 . Then βX/U T ≈ 0 and the model reduces to the usual Mullins-Sekerka model, which is invariant under rescaling with X = T 1/3 .

Gradient-flow structure
In this section we aim to explain the gradient-flow structure of the monopole model in (1.1). Following [5], we show that this structure is obtained by restriction of a corresponding one for the full Mullins-Sekerka model with kinetic drag, to interfacial geometry consisting of a collection of spheres with fixed centres. This method of geometric restriction seems a rather powerful way of obtaining reduced or simplified models while preserving useful structure. For example, it was used in a study of diblock copolymers to model the evolution of radii and centres [10]. The gradientflow structure of the monopole model will be used in the next section to obtain the coarsening estimates that we seek.

Mullins-Sekerka model as gradient flow
The Mullins-Sekerka model was described formally as gradient flow for surface area on a 'manifold' of sets with smooth boundary in [17] for the case without kinetic drag (β = 0) and with periodic boundary conditions. Here we recall how this works and indicate an extension to the case β > 0 without boundary conditions, meaning one phase comprises finitely many particles in all of R 3 .
First, we recall generally that the ingredients of a gradient flow involve a Riemannian manifold M with metric g and a functional A : Here dA(y) is the cotangent vector induced by the differential of A at y. For Mullins-Sekerka flow, we consider M to be the 'manifold' of all bounded sets G ⊂ R 3 with smooth boundaries and having a prescribed volume. Elements of the tangent space T G M correspond formally to (smooth) normal-velocity fields v : ∂G → R with zero mean. To any such v ∈ T G M we associate a harmonic potential u = J G v defined as the unique solution for the Neumann jump boundaryvalue problem Technically, u can be found as the minimizer of which is complete by a critical Sobolev embedding theorem. Here we omit discussing details [5].
We describe a metric as follows. Given v andṽ ∈ T G M, we let u = J G v and u = J Gṽ and define Note that integration by parts using (3.3) forṽ yields As energy functional A : M → R we take half of the total surface area: Then the differential dA(G) (formally an element of the cotangent space T * G M) is well known to correspond to the mean curvature of ∂G: for any smooth curve s → G s with G 0 = G and arbitrarily specified normal velocity fieldṽ on ∂G 0 , Suppose that t → G t is a solution trajectory for gradient flow of A on M with the metric in (3.5), with normal velocity field v on ∂G t . By (3.1) and (3.6) this Sinceṽ is arbitrary with zero mean on ∂G t , we can infer that the quantity is constant on ∂G t , and θ = θ(t) is a function of time alone. Averaging (3.9) over ∂G t , since v has zero mean we may write It then follows that with u = J G v + θ instead, equations (2.1) and (2.2) of the Mullins-Sekerka model with kinetic drag are satisfied. Note that the spatial constant θ is then the limit of u(x) as |x| → ∞.

Monopole model as constrained gradient flow
Interfaces can develop singularities during Mullins-Sekerka dynamics: for example, pinch-off of dumb-bell shaped regions and coalescence of neighbouring domains. To avoid the difficulties caused by such singularities, we now consider a constrained simpler geometry: one phase consisting of finitely many non-overlapping spherical particles. Note that singularities can still occur, as small particles can vanish in finite time and growing spheres collide. When one or more particles vanish, the system simply evolves from its current state with fewer particles. We will not address collisions here, since they can be neglected if particles are well separated (see [19] for a study of the effect of collisions).
In this section we aim to describe the monopole model (piecewise in time between vanishing events) as gradient flow for the restriction of the energy function A to a finite-dimensional submanifold N of M. The submanifold consists of sets G that consist of the union of a finite number of non-overlapping balls B i = B(x i , R i ) with fixed centres x i and variable radii R i , with total volume prescribed as before.
The tangent space T G N corresponds to zero-mean normal velocity fields constant on the boundary of each ball.
We may identify elements of T G N with vectors v = (v i ) with components v i ∈ R, and write The metric is identified as follows. In this situation, we can represent the harmonic potential u = J G v as a superposition of monopoles, having the form To be precise, let ψ(x) = min(1, 1/|x|) and let and ψ(x) → 0 as |x| → ∞. By superposition, the harmonic potential u = J G v for G = i B i takes the form where, for all i, We write the metric using (3.6) and the fact that by the mean value theorem we have Then it follows that On N , the energy functional takes the form 3.16) and the differential acting on an arbitrarily given velocityṽ is The general equation of gradient flow (3.1) is now written for allṽ ∈ T G N , meaning allṽ for which We conclude that the quantity is independent of i and depends only on t.
Since v i =Ṙ i this yields exactly the monopole approximation (1.1) to the Mullins-Sekerka model with kinetic drag. Voorhees [22] gives a detailed discussion of the monopole approximation without kinetic drag (β = 0). With u = J G v + θ instead, we note that the equations governing this approximation may be written An interesting implication of the present analysis is that θ and all the velocities v i are indeed determined uniquely by (3.19) together with the constraint whenever the balls B i are non-overlapping. It may not be easy to see from direct examination that the governing linear system of equations is non-singular. However, the gradient structure makes it evident: the metric is a positive definite quadratic form on the finite-dimensional tangent space. Thus, the gradient-flow equation, requiring g G (v,ṽ) + dA(G),ṽ = 0 for allṽ ∈ T G N , (3.23) clearly determines v = (v i ) ∈ T G N uniquely from G = G t , the collection of balls with centres x i and radii R i .

Coarsening rates for the monopole model
In this section we establish the coarsening-rate estimates (1.4) and (1.5) for the monopole approximation of the Mullins-Sekerka model with kinetic drag, in the setting of finitely many disjoint spherical particles. We will also discuss the meaning of these estimates, with the aim to describe conditions under which coarsening-rate crossover might be expected to occur and might be observed either in computation or experiment. We will apply the strategy of Kohn and Otto [11], which involves three key ingredients: an interpolation inequality relating the 'inverse length' E (the normalized energy in (1.3)) to a dual 'length' L; a dissipation inequality that controls the growth of L in terms of energy dissipation; and ordinary differential inequalities that have already been established [7,11].
The dissipation inequality will relate to the basic energy dissipation identity for the gradient flow equation (3.23), namely Since, for the monopole model, we find (with v = (v i ) = (Ṙ i ) and u = Jv up to any constant) that

Attachment-limited bound on coarsening rate
The two estimates (1.4) and (1.5) will be proved using two different dual lengths. The first dual length L 1 is defined as This was used in [7] to prove upper bounds on coarsening rates for systems of spherical particles evolving under the pure attachment-limited kinetics in (1.6). The desired interpolation inequality between E and L 1 is an immediate consequence of the Cauchy-Schwarz inequality The time derivative of L 1 satisfieṡ Rescaling time as t = βt/8, we obtain the desired dissipation inequality: Based on this and the interpolation inequality (4.4), the ordinary differential inequalities of [7, lemma 4.2] apply directly to provide the following estimate: with constants γ 1 andγ 1 > 0 depending only on p. Unscaling time, we obtain (1.4). This proves theorem 1.1(i).

Diffusion-limited bound on coarsening rate
We will establish the second coarsening bound (1.5) using a dual length L 2 that is defined through a kind of specialization of the arguments of [4]. Let us introduce an auxiliary function φ, found as the solution of the boundary-value problem [n · ∇φ] = 0 on ∂G, (4.8) where χ G is the characteristic function of G = B i , the collection of balls Then we compute Define L 2 by Then (4.14) Hence, with C 1 = (8π/15) 1/2 we obtain the interpolation inequality Taking the time derivative of L 2 2 , we obtain Hence, with C 2 = 2π, we obtain the dissipation inequality Applying [11, lemma 3] after a trivial rescaling, we obtain the estimate (1.5) on the coarsening rate of the 3D monopole approximation. This proves theorem 1.1(ii).

Discussion
In this section we discuss the meaning of the coarsening estimates (1.4) and (1.5), particularly concerning the physical interpretation of the time restrictions on these estimates in terms of the crossover time.
As mentioned in § 1, up to a constant independent of β, the (attachment-limited) bound in (1.4) is better than the (diffusion-limited) bound in (1.5) for T < β 3 , and vice versa for T > β 3 . In order for the bound (1.4) to be meaningful, it is necessary that the bound be valid well before the crossover time. That is, we need i.e. we need the characteristic radius of particles to be initially much smaller than β. This is consistent with the heuristic arguments in § § 1 and 2, where we indicated that, as long as R(t) β, the system is expected to behave like volume-preserving mean curvature flow, with R(t) growing proportional to (t/β) 1/2 .
The bound in (1.5) is the sharper one in the later stages anyway, so the validation time T 2 :=ĉ 3 2 L 2 (0) 3 is less significant. But, for (1.5) to become valid roughly at or before the crossover time, we need T 2 β 3 , or (4.20) Interpreting L 2 (0) as a typical length is complicated, however, by the balance between the two terms appearing in the numerator of (4.13). Let The relative size of these terms can be related to the Debye screening effect, interpreting the numerator of (4.13) as roughly the electrostatic energy of a uniform distribution of charge on G.
Following [17], we roughly estimate this effect as follows. Suppose that R is a typical radius of particles, and suppose that h is a typical distance between the nearest particles. We consider the simplest situation, when all particles occupy a cubic region Ω with side length a, Ω being divided into small cubes of side length h and the particles being located on vertices of small squares. In this situation, the total particle number N = (a/h) 3 and (roughly, up to constants) It is argued in [17] that the screening length ξ scr is (4.24) From these considerations we can make a few points. First, if the ratio II/I is large, then the interpolation inequality (4.4) is very pessimistic. The bounds in (1.5) can be expected to be roughly optimal only when II/I is of order 1 or less, and this means heuristically that the screening length should be at least of the order of the system size: ξ scr a. (4.25) If this is already true at the initial time, then heuristically L 2 (0) ∼ R(0), and so (4.20) holds provided that R(0) β. If, however, (4.25) fails to hold and we have ξ scr a instead, then Note that, since NR 3 ∼ const. by volume conservation, (4.23) indicates that the quantity in (4.26) is roughly constant in time and equal to √ φ, where φ is the volume fraction of the particles in Ω. Now, to have T 2 β 3 , (4.20) requires a stricter condition than (4.19), namely, . (4.27) If this fails and β 3 T 2 instead, we anticipate a gap between the crossover time β 3 and the validity time T 2 . In the gap, heuristically we expect coarsening with R ∼ t 1/3 , but the bound (1 .5) is not yet valid. In this case, we expect that, about the time when ξ scr approaches a we will obtain R ∼ L 2 ∼ L 2 (0), and this is consistent with T 2 ∼ R 3 . What this means is that, heuristically, the (diffusionlimited) bound (1.5) can only be expected to become valid after the screening length becomes of the order of the system size. In this sense, improvements in the argument of § 4.2 would be desirable from a physical viewpoint.
We anticipate that it is likely to be very difficult to observe the expected crossover behaviour in numerical simulations of the monopole model. To obtain good statistics one should have to compute with very large numbers of particles for very long times. The reason can be indicated roughly as follows. With the normalization that the typical radius is initially R(0) ∼ 1, we expect it to be necessary to take β an order of magnitude larger and compute until R(t) is an order of magnitude larger than that, or R(t) ∼ β 2 . Since we expect R ∼ t 1/3 after the crossover time T ∼ β 3 , the time it takes for this is long: of the order t ∼ β 6 . Moreover, the initial number of particles must be taken as very large to end up with significant numbers of surviving particles, because the ratio of surviving particles to initial particles will be N/N (0) ∼ R(0) 3 /R 3 ∼ β −6 by volume conservation.

LSW model with kinetic drag
Because of the difficulty of simulating the monopole model directly, we consider here a simpler model, the Lifshitz-Slyozov-Wagner mean field model (1.7) In this model the particles interact with each other only through a spatial mean field θ = θ(t) which is determined by conservation of total volume as before. From The LSW model is justified in the extremely dilute limiting case when the capacity density of the particles approaches zero, which is a stronger condition than a disappearing volume fraction [15,16]. This LSW model has a gradient-flow structure, described in [18], that we can also obtain by simplification from the monopole model: letÑ be the collection of vectors R = (R i ) of particle radii with fixed total volume. The tangent space T RÑ corresponds to all possible vectors v = (v i ) of normal velocities such that Define a metric on T RÑ by 3) Then the system (1.7) with (5.2) is the gradient flow for the rescaled surface area with respect to the metric g.

Coarsening rate estimates
Now we establish bounds on the coarsening rates for the LSW model. Corresponding bounds for diffusion-limited and attachment-limited LSW models were separately established in [7]. But the argument here is simpler, as we make use of the gradient-flow structure. As before, we still work with the volume-averaged interfacial area Also, we need a dual length-scale. The situation here is slightly simpler than for the monopole model. For both the attachment-limited and diffusion-limited regimes, it suffices to consider By the Cauchy-Schwarz inequality we obtain the interpolation inequality EL 1 as before. Now consider the dissipation relation. We have combining (5.6) and (5.7) we obtain both kinds of dissipation inequalities that appeared in the previous section: Using the interpolation inequality and each of these dissipation inequalities in turn with the ordinary differential inequalities in [11, lemma 3]

Numerical simulations
To illustrate the coarsening behaviour, we performed numerical simulations for the LSW model with various values of β. We use the forward Euler method to discretize in time. The time steps ∆t n are chosen to be no bigger than dt = 0.5.
It is not convenient to directly discretize equation (1.7) because this equation indicates that the radii of disappearing particles shrink to zero at infinite speed, which is very hard to simulate accurately. On the other hand, the rescaled particle volumes x i := R 3 i shrink to zero at a finite speed even if β = 0. So we will discretize the following evolution equation for particle volumes x i : where Let x n i be the volume of the ith particle at time t n . The scheme is as follows: where θ n is determined by (5.13) with x i replaced by x n i . Note that particles can disappear in finite time and we need to estimate these disappearing times accurately. To do so, we compute the minimum time step that makes some particle size vanish in the system (5.14), and take ∆t n to be the minimum of this value and dt. Then we discard those particles with radius near zero from the system and continue the scheme with the surviving particles.

Initial data
To motivate our choice of initial values we mention that the LSW model is usually formulated to consider the evolution of a particle volume distribution f (t, x), which satisfies a transport equation Here θ(t) is defined as in (5.13) but with summations replaced by integrals with respect to the distribution f . We see that equation (5.12) is the corresponding characteristic equation for this partial differential equation.
Our arguments indicate that, in the early stage, when typical particle radii are much smaller than β, the LSW model behaves similarly to a volume-preserving mean curvature flow, governed by The corresponding transport equation for the particle volume distribution is The corresponding solution for two dimensions was described and used in [3] to compare with experiments concerning island growth on silicon surfaces. Since the goal of our simulation is to observe crossover behaviour, it is natural to start with a particle size distribution that comes from a self-similar solution for the volume-preserving mean curvature flow. We first discard the tail (1 − δ, 1] of the distribution near x = 1, since f 0 (x) exponentially decays to zero as x approaches 1, then uniformly divide the interval [0, 1 − δ] into subintervals of size h and use the node points as a set of particle volumes {x 0 0 = 0, x 0 1 = h, x 0 2 = 2h, . . . }. For each x 0 i , we assign a weight m i := f 0 (x 0 i ) which indicates the relative number of particles of volume x 0 i that we have. We choose δ = 0.15 and h = 10 −4 .

Results
The simulation approximates, for each i, the volume x i of the m i particles at later times. So the scheme is exactly (5.14), with the exception that θ n is determined by weighted sums:  A small but important point in analysing the data is that, since the systems are translation invariant in time, a temporal power law E(t) ∼ t −α should be more precisely represented as E(t) ∼ c(t + t 0 ) −α for some constant c and time translation t 0 . To best match the form of the solution in (5.18) and (5.19) at early times, we take t 0 = 1 2 β. Figures 1 and 2 show the behaviour of the log of total surface area E(t) plotted against log(t + t 0 ) for the two values β = 20 and 30. It is clear that E behaves similarly to (t + t 0 ) −1/2 at early times and similarly to (t + t 0 ) −1/3 ∼ t −1/3 later.