On the mathematical modelling of cellular (discontinuous) precipitation

Cellular precipitation is a dynamic phase transition in solid solutions (such as alloys) where a metastable phase decomposes into two stable phases : an approximately planar (but corrugated) boundary advances into the metastable phase, leaving behind it interleaved plates (lamellas) of the two stable phases. The forces acting on each interface (thermodynamic, elastic and surface tension) are modelled here using a first-order ODE, and the diffusion of solute along the interface by a second-order ODE, with boundary conditions at the triple junctions where three interfaces meet. Careful attention is paid to the approximations and physical assumptions used in formulating the model. These equations, previously studied by approximate (mostly numerical) methods, have the peculiarity that \begin{document}$v,$\end{document} the velocity of advance of the interface, is not uniquely determined by the given physical data such as \begin{document}$c_0$\end{document} , the solute concentration in the metastable phase. It is hoped that our analytical treament will help to improve the understanding of this. We show how to solve the equations exactly in the limiting case where \begin{document}$v=0$\end{document} . For larger \begin{document}$v$\end{document} , a successive approximation scheme is formulated. One result of the analysis is that there is just one value for \begin{document}$c_0$\end{document} at which \begin{document}$v$\end{document} can be vanishingly small.

Abstract. Cellular precipitation is a dynamic phase transition in solid solutions (such as alloys) where a metastable phase decomposes into two stable phases : an approximately planar (but corrugated) boundary advances into the metastable phase, leaving behind it interleaved plates (lamellas) of the two stable phases.
The forces acting on each interface (thermodynamic, elastic and surface tension) are modelled here using a first-order ODE, and the diffusion of solute along the interface by a second-order ODE, with boundary conditions at the triple junctions where three interfaces meet. Careful attention is paid to the approximations and physical assumptions used in formulating the model.
These equations, previously studied by approximate (mostly numerical) methods, have the peculiarity that v, the velocity of advance of the interface, is not uniquely determined by the given physical data such as c 0 , the solute concentration in the metastable phase. It is hoped that our analytical treament will help to improve the understanding of this.
We show how to solve the equations exactly in the limiting case where v = 0. For larger v, a successive approximation scheme is formulated. One result of the analysis is that there is just one value for c 0 at which v can be vanishingly small.
1. Introduction. One of Paul Fife's lifelong interests was the mathematical modelling of dynamic phase transitions in solids. Here we consider such a phase transition, the separation of a spatially uniform metastable phase of a two-component alloy into two separate stable phases. For example, an alloy consisting of 30%Zn atoms and 70%Al atoms can separate into an Al-rich phase in which the concentration of solute Zn atoms is at most 5% (depending on temperature) and a Zn-rich phase in which at the concentration of solute Zn atoms is at least 60% . Thermodynamic data tell us what new phases are possible, but we would also like to understand the shapes, sizes and positions of the new phases.
Of the various possible configurations, one of the most straightforward, from the theoretical point of view, is cellular precipitation (also known as discontinuous precipitation) in which an approximately planar (but corrugated) boundary advances into the metastable phase, leaving behind it interleaved plates ("lamellas") of the two stable phases, as illustrated schematically in Fig. 1. The as yet undisturbed metastable phase is labelled 0; the two (relatively) stable phases into which it separates after the front has passed are labelled α and β.
This paper is concerned with the mathematical modelling of cellular precipitation, a problem which has so far received little attention from mathematicians. Previous treatments [5,10,3,4] encountered the difficulty that the solution of the equations used is not unique, and an arbitrary criterion was used to pick out the 'right' solution. In an effort to understand this situation better, we provide more careful derivations of some of the basic formulas of the model and try to solve them without some of the rather inaccurate approximations that were found necessary in the earlier papers.
2.1. The geometry. The geometry is assumed two-dimensional, i.e. everything is independent of a 'depth' coordinate z perpendicular to the plane of Fig. 1. The triple junctions (places where three phases meet) are assumed to be straight lines perpendicular to the plane of the diagram and to lie in a plane. The intersection of this latter plane with the plane of the diagram is a straight line. This line will be used as the x axis -see Fig. 2.
The cross-sections of the interlamellar boundaries, shown straight in Fig. 1, may in practice be curved, as shown in Fig. 2. Fig. 2 also shows the notation that will be used for describing the geometry of the specimen. We take the origin of coordinates at one of the triple junction points and position the x axis so that it contains all the triple junctions. The distance from the origin shown to the next triple junction on its right is denoted by X α and the distance to the next triple junction on its left is denoted by X β . If the specimen is taken to be infinite, the structure will be periodic, with period X α + X β . The sign convention for angles is that the angles marked as in Fig. 2, with origin at a triple junction having α phase just to the right it and β phase just to the left, are positive. All arc lengths are denoted by s and measured from the nearest triple junction.
A full list of symbols and their definitions is given in the Appendix.
y axis Figure 2. Notation for the coordinate axes and for angles and arc lengths 2.2. Thermodynamics. The specimen is assumed to consist of two types of atoms, A (the solvent) and B (the solute). The mole fraction of solute (percentage of B atoms) in a macroscopic neighbourhood near any given point will be denoted by c. The mole fraction in the initial metastable phase will be denoted by c 0 ; after the phase separation, the mole fraction of solute in the α (A-rich) phase will be denoted by c α , and in the β (B-rich) phase by β. Both c α and c β may depend on position, and c 0 (which does not depend on position) lies between them : The thermodynamic properties of such an alloy system can be summarized in a free energy diagram like the one in Fig.3. The curves shown are two pieces of the graph of the function f (·), where f (c) is the Helmholtz free energy per unit volume at solute concentration c. We denote by c eq α and c eq β the concentrations of the two phases when they are in equilibrium with one another across a plane boundary. By a well-known construction, c eq α and c eq β are the c-coordinates of the places on the graph where a common tangent touches the two curves, as shown in Fig. 3. The slope of the common tangent will be denoted by µ eq , so that 3. The equations of our model.

3.1.
The forces acting on an interface. The interfaces (grain and phase boundaries) are assumed to move in the positive y direction (as shown in Fig. 1) at a uniform velocity v. The force that moves them is of thermodynamic origin, arising from pressure differences (denoted here by p) between the phases on opposite sides of the interface. In addition, we must take into account elastic forces, arising from  . The thermodynamic construction for phase equilibrium. The common tangent has slope µ TJ and touches the curve at the points with abscissas c eq α and c eq β .
the difference in the volumes occupied by the two types of atom. The importance of elastic forces for moving grain boundaries was pointed out by Hillert [9] in 1983; a 1997 paper co-authored by Paul Fife [6] showed that for the related phenomenon of diffusion-induced grain boundary motion the elastic force was crucial. The elastic force is an essential part of the theoretical treatment of cellular precipitation given by Brener and Temkin [4]; indeed all the physical assumptions used in the present paper are essentially the same as those used by Brener and Temkin. A formula for p, taking into account the elastic forces, will be derived in sub-section 3.4. An additional force on the interface is due to surface tension. Surface tension exerts a force per unit area, in direction normal to the interface, with magnitude equal to γ times the curvature, where γ denotes the surface energy per unit area. We shall assume that the velocity component of the grain boundary in the normal direction is equal to the net force force per unit area tending to move the grain boundary, multiplied by a a mobility factor M, which for simplicity we take to be the same for all three types of interface considered in the model. Using the sign conventions shown in Fig. 2, this modelling assumption gives the equations v cos θ = (p + γdθ/ds)M for the 'corrugated' 0α and 0β interfaces (3) v sin ϕ = (p + γdϕ/ds)M for an αβ interface At the triple junction the forces exerted by the three interfaces meeting there must be in mechanical equilibrium; therefore, by the sine rule, the surface tensions in these interfaces are proportional to the sines of the dihedral angles opposite to them: where δ 0 is the dihedral angle in the 0 phase, γ αβ is the surface energy per unit area in the αβ interface, and so on. The slopes of the tangents to the interfaces at a triple junction with α phase to its right are related to the dihedral angles (see Fig.4) by the formulas where, for each arc, the argument of a function such as ϕ(·) is the arc length measured from the nearest triple junction. 3.2. The transport of solute atoms. In order to form the new phases, solute atoms must move from place to place as the corrugated interface passes. This transport takes place by diffusion. It will be assumed (following, for example, Cahn [5]) that the diffusivity is much greater in the interfaces than in the bulk material, so that the interfaces are the highway along which solute atoms travel. For a quantitative decription, we shall use the description of diffusion in an elastic solid due to Larché and Cahn [11]. According to their treatment, diffusion is driven by the diffusion potential, µ which, if the material is isotropic and the composition and strain depend on only one of the three position coordinates, is equal to where f (c) is the Helmholtz free energy per unit volume at solute concentration c, f (c) means df (c)/dc, which is called the stress-free diffusion potential and is analogous to the chemical potential, c is the local concentration, c ∞ is the concentration in the same phase a long distance away, Y is an elastic modulus which for the geometry used here is equal to E/(1 − ν); in this last formula E is Young's modulus, ν is Poisson's ratio, and η is the coefficient in Vegard's law, according to which the strain due to the misfit when B atoms are dissolved in the material at zero stress is proportional to c In this theory the flux of solute atoms is the vector where B is a mobility coefficient (not to be confused with the mechanical mobility M in eqns (3,4)). This mobility is related to the bulk diffusivity D by so that (8) reduces to Fick's law J = −D∇c when the stress is zero (in which case, by (7), µ = f (c)). But if the stress is not zero, eqn (7) implies dµ/dc = f (c)+2Y η 2 and Fick's law takes the more general form The diffusion equation implied by (8) and (10) is where In a region small enough for the interface to be treated as planar, take Cartesian coordinates s, z, n where s is distance in the tangential direction, z is distance from the plane of the cross-section in Fig. 1, and n is distance away from the interface, positive in the upwards direction. The interface will be treated as a channel of width λ, advancing in the n direction at rate v cos θ and in the s direction at rate s sin θ. We assume that the diffusivity in the channel is large, but that the diffusivity elsewhere in the material is very small. To quantify these assumptions more precisely, we assume that the diffusivity depends on n in the folowing way where D 0 and D 1 are constants. At a suitable moment we shall take the limit λ → 0. Assuming a steady state, c and µ will be functions of n t := n − vt cos θ and s t := s − vt sin θ only (t being the time), so that ∂/∂t = −v cos θ ∂/∂n − v sin θ ∂/∂s. On taking t = 0, so that n t = n and s t = s, the diffusion equation (11) takes the form Introducing the re-scaled variable ν := n/λ, eqn (14) becomes, on multiplying by λ and then using (13) −v cos θ ∂c ∂ν − λv sin θ ∂c ∂s = λ ∂ ∂s D * ∂c ∂s On taking the limit λ → 0 in (16), and assuming that all derivatives are of order 1 or less (rather than, say, of order 1/λ) we obtain The boundary conditions to be satisfied at large ν are where c s (s) and c g (s) are, respectively, the concentrations in the shrinking and growing grains at a distance from the grain boundary that is large compared to D 0 /λv but small compared to the radius of curvature of the grain boundary. The solutions satisfying these conditions are where the function g(·) is to be determined from the boundary conditions near the interface. Thus, the solute diffuses ahead of the advancing interface a distance of order D 0 /λv, as first shown by Balluffi and Cahn [2] for diffusion-induced grain boundary motion. Behind the interface there is no elastic strain and no diffusion in the normal direction. Now we consider the interior of the channel, where |ν| < 1 2 . The analogue of (16) is now where the definition of D * 1 is analogous to that of D * 0 in (17). In the limit of small λ the coefficient of 1/λ 2 must go to zero, so that the equation tells us that D * 1 ∂c/∂ν is independent of ν. In other words, according to the definition (17), (1 + 2Y η 2 /f (c))∂c/∂ν is independent of ν or, going back still further, to (7). ∂µ/∂ν is independent of ν.
The approximate solutions (21, 22) imply that both the partial derivatives on the right are of order 1 (indeed one of them is actually zero); therefore the ones on the left must be of order λ 2 , and are zero in the limit of small λ. Thus, the constant value of ∂µ/∂ν inside the channel must be zero. It follows that the diffusion potential inside the channel is independent of ν and is the same on the two sides of the channel.
To evaluate the integral in (26), we use the facts that µ is independent of ν in the channel, and that c is independent of ν behind the interface so that, by eqn (7) µ(s) = f (c g (s)) in the growing grain. Hence, by the definitions of D * 1 and µ (eqn (7)), we have Thus, if we make the approximation f (c) ≈ f (c g ), eqn (26) simplifies, in the limit λ → 0, to This equation is due to Cahn [5].

3.3.
A sum rule. If we integrate eqn (29) with respect to s along the intersection of the 0α interface with the plane of the diagram the result is where c (s) := dc/ds and 2S α is the arc length of the interface. Using the assumed symmetry shown in Fig. 2 and the fact that cos θ = dx/ds, eqn (30) simplifies to where 2X α denotes the thickness of a lamella of the α phase, measured betwen two triple junctions on opposite faces of the lamella. Adding together (31) and the corresponding formula for the 0β interface, and assuming for simplicity that D 1 is the same for all three types of interface, we obtain a result which can be simplified further using the formula which expresses the fact (noted by Brener and Temkin [3,4]) that the net flow into the triple junction along the three interfaces that meet there must be zero. In this way we obtain the sum rule This sum rule expresses the overall rate of accumulation of solute atoms in the two new phases as the sum of the rate at which the old phase is being eaten up and the rate of 'recycling' arrival along the αβ interface.
3.4. The pressure. To use eqns (3) and (4) we need a formula for the local pressure force on an interface. To calculate this we shall treat the interface as locally planar, as in the preceding section. The force on the interface is the space rate of decrease of the total Helmholtz free energy (including elastic energy) when the interface moves, the solute atoms redistributing themselves so as to bring about the required concentration changes. Solute atoms must be supplied or removed as the interface passes by. These atoms enter or leave the relevant part of the specimen by travelling along the interface, whose chemical behaviour will be represented, for the pressure calculation, as a (stress-free) reservoir whose stress-free diffusion potential µ equals the diffusion potential in the interface. As shown in the paragraph containing eqns (24,25) the diffusion potential in the channel depends only on s and is continuous across the interface, so we can treat µ at the interface as a function of s only. The free energy of the reservoir is equal to a constant plus µ times the number of solute atoms in it and therefore, since mass is conserved, equal to a constant minus µ times the number of solute atoms in the specimen. Thus the total Helmholtz free energy of specimen and reservoir together, which will be called here the augmented free energy, at a place where the solute concentration is c, is a constant plus per unit volume of specimen, where w denotes the elastic energy per unit volume. According to Larché and Cahn [11], the formula for w when (as here) the composition and strain depend on only one of the three Cartesian position coordinates is where the notation is defined after eqn (7). The derivative of f (c) + w with respect to c is the diffusion potential, given in eqn (7). The augmented free energy in the growing grain at the point just behind the grain boundary with arc coordinate s is, by eqn (34), per unit volume. In the shrinking grain, at the corresponding point just in front of the grain boundary, the augmented free energy per unit volume, calculated using eqn (34) and (35), is The pressure is the difference between these two augmented free energies To use this formula we need to know µ and c i Since the stress in the growing grain is zero, the diffusion potential in the interface is equal to the diffusion potential in the growing grain so that, by (7): At the corresponding point just across the grain boundary the diffusion potential is the same but the solute concentration, which se shall call c i , satisfies

OLIVER PENROSE AND JOHN W. CAHN
by virtue of (7). Thus the equation for c i is This equation can be solved using our usual approximation of treating the relevant part of the free energy curve of the shrinking grain (see Fig. 3) as a parabola. The equation then becomes so that Using this formula, along with eqn (39), in the right hand side of the pressure equation (38), we obtain With a view to approximating this expression by means of a Taylor expansion in powers of c g − c eq g and c s − c eq s , let us evaluate the relevant partial derivatives of p. The derivatives we need are those of The results of the Taylor expansion are different according to whether the phases on the two sides of the interface are the same or different. If the two phases are the same, so that c eq g = c eq s =: c eq , then the zero-order approximation to p is zero and so are the equilibrium values of both the first derivatives in (45); the second-order approximation for p 1 turns out to be To this must must be added the lowest-order approximation to the quadratic term in (44), whose calculation makes use of the assumption that the free energy curve of the shrinking grain is a parabola, so that f (c g ) − f (c s ) ≈ (c g − c s )f (c eq s ). The resulting formula for p is If Y η 2 f (c eq g ) this reduces to Hillert's [9] formula for the elastic driving force in diffusion-induced grain boundary motion; a less laborious derivation of a formula essentially the same as (47) can be found in ref. [15]. Now consider the alternative case, where the phases on the two sides of the interface are different. Since c s , the solute concentration in the shrinking grain, may vary strongly with position near the interface, an expansion of c s in powers of c s −c eq s may not be appropriate. We can still expand the right side of (44) in powers of c g − c eq g , however. This expansion gives, using the formula for ∂p/∂c g from (45) 3.5. Summary of the equations.
• The force equation (3) v cos θ = (p + γdθ/ds)M with boundary conditions (6) for the 0α and 0β interfaces where S α or S β denotes the half-length of the relevant arc in Figs 1 and 2.
with the boundary condition where µ TJ denotes the diffusion potential at the triple junction -it is the same for all three interfaces meeting there. For the 0α and 0β interfaces there is the further condition • The triple junction net flow condition (32) or, alternatively, the sum rule (33) • The pressure formulas (47) and (48), if the phases on opposite sides of the interface are the same, and if they are different.

4.
Solving the equations. The system of equations summarized in subsection 3.5 has generally been tackled by numerical methods, involving some questionable approximations such as replacing cos θ by 1 in Cahn's eqn (29). The approach was to fix c 0 and then try to calculate the other variables. The method does not yield a unique solution, but rather a one-parameter family of solutions, which could be parameterized by the value of v. Experimentally, however, the value of v for given c 0 and material parameters does appear to be unique. The physical mechanism by which v is selected is not well understood. Our aim here is to understand the mathematics of the model better; for example there might be some instablitiy that rules out the unwanted values of v.
Our approach to the mathematical problem is different from those used previously. We approach it from the opposite direction: rather than treating c 0 as given we fix the values of two other variables, ϕ(0) and v, and try to express the remaining variables, including c 0 , in terms of them.
4.1. The αβ interface. We begin by looking at the equations for the αβ interface, which trails behind the advancing front shown at the top of Figs 1 and 2. Using the notation ϕ := θ − 1 2 π (see Fig 2) and the approximation sin ϕ = ϕ(1 + O(ϕ 2 ), eqns (49) and (51) become If ϕ > 0 as shown in Fig 2 the shrinking grain belongs to the α phase and the growing grain to the β phase; if, on the other hand, ϕ < 0 then it is the growing grain that belongs to the α phase. In either case the pressure is given by eqn (56), a formula which can be approximated using a first-order Taylor expansion when c g and c s are both close to their thermal equilibrium values. The resulting formula can be written In the same spirit, we can approximate the right-hand side of (58) to lowest order in c g − c eq g and c s − c eq s , which for that expression simply means replacing c g and c s by their equilibrium values, so that (58) is approximated as Combining this with (59) we find that Eliminating p between (61) and (57) we obtain the following constant-coefficient linear differential equation for ϕ, A similar equation was obtained by Mullins [13] for a problem involving a solidvapour interface. Mullins' equation is also used by Brener and Temkin [3]. Equation (62) has solutions of the form ϕ = e −qs where q is a solution of the characteristic equation The cubic equation (63) has one positive solution and two others which are either both negative or complex with negative real parts. Only the positive solution of (63) leads to a solution of (62) that remains bounded as s → ∞. For small v the second term on the left is negligible and so the solution that is bounded as s → ∞ is ϕ = ϕ(0)e −qs (64) where q = f (c eq g )(c eq g − c eq s ) 2 Dλγ To proceed to the next stage of the calculation, we need formulas for the diffusion potential and the flow at the triple junction, when v and ϕ(0) are known. The appropriate solution of eqn (60) with ϕ given by (64) is from which, using (39), the diffusion potential at the triple junction can be calculated : the last formula being exact if the free energy curve is a parabola. For the triple junction flow condition (32) we need the derivative of c g which, by (66), is where Γ := γ f (c eq g )(c eq g − c eq s ) f (c eq g )(c eq g − c eq s ) 2 Dλγ It may be worth noting that, since ϕ(0) and c g − c s have the same sign, µ TJ − µ eq is positive (or zero) and c αβ (0) is negative (or zero), regardless of the sign of ϕ.

4.2.
The 0α and 0β interfaces. The next step is to solve the equations for the 0α and 0β interfaces, using the boundary conditions determined by the given value of ϕ(0) and the formulas (66) and (68). We consider first the easiest case, in which v = 0, and then show how to extend the results to larger values of v 4.2.1. The v = 0 case. Eqn (51) shows that, when v = 0, each of c α and c β must vary linearly with s, and by the symmetry assumed in Fig. 2, the linear function must be a constant. The value of this constant is determined by the condition that the diffusion potential at the triple junction is the same for all the three interfaces meeting there 1 Eqn (67) shows that when v = 0 the diffusion potential has its equilibrium value µ eq , and hence by (2), the concentrations in the growing grains at the 0α and 0β also have their equilibrium values. Thus, the constant value of c on the 0α interface is c eq α , and on the 0β interface it is c eq β . Denoting the zero-velocity values of the variable c by c (0) , we can write this concisely as Here we are using the notational convention that in an equation such as (70), containing variables such as c, p, θ etc with no subscript, either of the subscripts α and β can be inserted; thus eqn (70) stands for the pair of equations c α (s) = c eq α , c β (s) = c eq β , meaning c (0) (s) = c eq α on the 0α interface c (0) (s) = c eq β on the 0β interface (71) The zero-velocity pressure differences across the 0α and 0β interfaces are also constants, being given, according to (55) and (56). by Setting v = 0 in the force equation (49) we find that so that each arc is circular, with radius 1/κ, i.e. γ/p is the total length of the arc when v = 0. The condition for determining S, assuming the symmetry shown in Fig. 2, is This result enables us to calculate the thicknesses of the lamellas, to the current order of approximation, using a formula already implicit in eqn (33): as in (6). At this stage we would like to be able to use the sum rule (33) to determine a value for c (0) , but this is impossible for the moment because, since we are taking v = 0, eqn (33) just says 0 = 0. To get informatiohn about c 0 we need to go beyond the v = 0 case.

4.2.2.
Improving the approximation. We can improve on the v = 0 solution (or any other approximate solution)by the following successive approximations prodedure which, it is to be hoped, will converge if v and ϕ(0) are sufficiently small. The idea is that if the nth order approximations to variables such as c(s), p(s), θ(s) and S (denoted by c (n) (s), etc.) are accurate to n th order in the two 'independent' variables v and ϕ(0), meaning that the error is o(v n +ϕ(0) n ) as v → 0 and ϕ(0) → 0, then the next approximation will be accurate to (n + 1)th order (meaning that the error is o(v n+1 + ϕ(0) n+1 ). The approximation scheme consists of the following steps: • Given the nth order approximations to c(s), p(s), θ(s) and S for both the 0α and 0β interfaces, we obtain the next approximation to c(s) (treating c 0 as given for the time being) by solving the following approximate version of (51) with the boundary conditions f (c (n+1) (0)) = f (c (n+1) (2S (n) )) = µ TJ with µ TJ as given in (67), i.e.
µ TJ = µ eq + γq c eq g − c eq s ϕ(0) The explicit solution of this boundary-value problem is c (n+1) (s) = c eq + v Dλ ). If c (n) , θ (n) and S (n) are all accurate to nth order (as defined above) then the new approximation c (n+1) will be accurate to (n + 1)th order because of the factor v in front of the integral. • The next approximation to p can now be calculated from the following formula, which is based on (55) and (56): • The next approximation to θ(s) is obtained from the following formula, based on eqn (49) with the initial condition from (6), i.e.

OLIVER PENROSE AND JOHN W. CAHN
• Using this solution, we use the end condition with θ(0) as in (88), to find S (n+1) . • To get a better approximation, increase n by 1 and repeat the preceding sequence of steps • At the end of the calculation, we obtain the best available approximation to c 0 by using the sum rule (33) with c αβ (0) given by (68) 4.2.3. Illustration of the procedure. As an illustration, let us calculate a first-order approximation, using the v = 0 solution in eqns (70) to (76) as our zero-order approximation.