## Abstract

The moving interface problem in a one-dimensional binary α/β diffusion couple is studied using sharp and diffuse interface (Cahn–Hilliard) approaches. With both methods, we calculate the solute field and the Kirkendall marker velocity and displacement fields. In the sharp interface treatment, the velocity field is generally discontinuous at the interphase boundary, but can be integrated to obtain a displacement field that is continuous everywhere. The diffuse interface approach avoids this discontinuity, simplifies the integration and yet gives the same qualitative behaviour. Special features observed experimentally and reported in the literature are also studied with the two methods: (i) multiple Kirkendall planes, where markers placed on the initial compositional discontinuity of the diffusion couple bifurcate into two locations, and (ii) a Kirkendall plane that coincides with the interphase interface. These situations occur with special values of the interdiffusion coefficients and starting couple compositions. The details of the deformation in these special situations are given using both methods and are discussed in terms of the stress-free strain rate associated with the Kirkendall effect.

## 1. Introduction

The Kirkendall effect involves both diffusion and deformation. In general, the deformation of a material can be defined by the change in position with time of a set of inert markers (real or imagined) embedded in the material. The Kirkendall effect refers to changes in the marker positions that are observed as markers move in response to diffusion in the material. For diffusion that occurs primarily in a single direction, markers that are initially aligned in a plane perpendicular to the diffusion direction tend to move in concert, defining marker planes as the material deforms unidirectionally.

Fundamental to understanding the Kirkendall effect is the definition of independent diffusion fluxes for the components with respect to this set of markers (intrinsic fluxes). As shown in detail by Stephenson (1988), the deformation that accompanies diffusion is described by a stress-free strain rate. Although the Kirkendall effect has been known for 60 years (Smigelskas & Kirkendall 1947; Darken 1948), a number of interesting theoretical and experimental effects continue to be revealed which warrant further consideration, including: (i) multiple Kirkendall planes, (ii) stability of individual Kirkendall planes, and (iii) discontinuity in the Kirkendall velocity at moving interphase boundaries.

Diffusion couples are widely used for experimental and theoretical studies. For a binary alloy composed of components A and B maintained at a constant temperature, the initial state consists of two end members joined with an initial planar discontinuity at a position we take to be *z*=0. The end members are uniform and of different compositions that may also be of different phases. When the couple is held at a fixed temperature, diffusion occurs over time in the *z*-direction. Experimentally, markers are typically placed only on the initial discontinuity of the diffusion couple; however, they can be placed throughout the sample to measure the entire field of deformation field in the *z*-direction (e.g. Frerichs & Ruth 1997; van Dal *et al.* 2000). The time rate of change of the position of the entire set of markers in the sample determines the ‘Kirkendall velocity’ as a field defined throughout the sample.

In models, one has a choice of describing this velocity field as a function of position *z* and time *t* (Eulerian, laboratory or spatial description) or as a function of position at time zero *Z* and time *t* (Lagrangian, reference or material description). The former is commonly used for fluids, the latter for solids. In this paper, we use a spatial description for the concentration and velocity fields, but the computation of displacements will require reconsideration of the material description.

We note that in the analysis of diffusion couples, there are several quantities of interest, each of which is associated in some way with the initial planar discontinuity. One or more Kirkendall planes issue from *z*=0 at *t*=0 and represent the motion of inert markers tracking the deformation, i.e. a Kirkendall plane is a material surface. In this paper, we treat end members that are of different phases so that a phase boundary also issues from the point of initial discontinuity, but it usually does not coincide with a Kirkendall plane once diffusion occurs, since the interphase boundary is not a material surface. We will see that, depending on the relative size of the diffusivities of A and B, markers may overtake or recede from the phase boundary, or in some circumstances even accumulate on the interface. In addition, the initial discontinuity is also the site of the Matano plane, which plays a role in the detailed analysis of diffusion (Philibert 1991). The Matano plane is defined by appropriate equal-area constructions based on the solute profile (see appendix A), and remains at *z*=0 for all times. Thus, the Matano plane is also not a material surface. The phase boundary, the Kirkendall plane(s) and the Matano interface are generally associated with different locations in the sample once diffusion takes place.

This paper is stimulated by a variety of theoretical and experimental efforts. Kirkendall markers (i.e. those markers placed initially in the plane *z*=0) are occasionally found on two different planes after heat treatment. The possibility of such multiple Kirkendall planes was first described theoretically for single-phase diffusion couples by Cornet & Calais (1972), using similarity solutions for the composition and Kirkendall velocity fields in a binary diffusion couple. A necessary ingredient for multiple Kirkendall planes is a strong variation of intrinsic diffusion coefficients with composition; for example, multiple Kirkendall planes can occur if the difference in diffusivities changes sign with composition across the diffusion couple. The possibility of multiple Kirkendall planes was further discussed theoretically by van Loo *et al.* (1990), using a partitioning (physico-chemical) analysis of multiphase diffusion couples that is useful for phases of narrow composition range. In multiphase diffusion couples, the Kirkendall velocity is usually discontinuous at moving interphase interfaces. This requires that the moving interface is an extremely effective source/sink of vacancies. In addition, the likelihood of finding multiple Kirkendall planes is higher than in single-phase couples because the fast diffusing species may be different in different phases of the same binary system.

Experimentally, van Dal *et al.* (2000) observed two Kirkendall planes in different phases of a multiphase Ni–Ti couple that met the theoretical conditions prescribed by Cornet & Calais (1972). In addition, van Dal *et al.* (2000) showed experimental results for single-phase couples of Ni–Pd and Fe–Pd with solitary Kirkendall planes, a result that also agreed with the criterion established by Cornet & Calais (1972). The Ni–Pd and Fe–Pd experimental results were reanalysed by Höglund & Ågren (2001) using the Dictra1 code (Ågren 1982; Andersson *et al.* 2002), confirming the conclusions of van Dal *et al.* (2000). Two Kirkendall planes were observed for a Au–Zn multiphase couple by van Dal *et al.* (2001) and for a Co–CoSi_{2} multiphase couple by Paul *et al.* (2004*b*). Double Kirkendall planes were found within a single phase, NiAl, by Paul *et al.* (2004*a*). In this case, the identity of the fast diffusing species switches on either side of perfect stoichiometry due to the change in defect structure of the ordered B2–NiAl compound.

van Dal *et al.* (2001) observed another interesting effect dealing with the so-called stability of Kirkendall planes. The term ‘stability’ in this context refers to whether markers placed on the plane of the initial discontinuity of a diffusion couple remain sharply concentrated on a plane or whether they become dispersed in the diffusion direction. They related stability or instability to the sign, negative or positive, respectively, of the spatial derivative of the Kirkendall velocity at the position of the Kirkendall plane.

The purpose of this paper is to describe the peculiarities of the deformation in a one-dimensional binary diffusion couple with a moving α/β interface under various conditions. We compute the trajectories (position versus time) of a full set of markers, focusing on the cases with multiple Kirkendall planes and Kirkendall planes that coincide with the moving interface. Computation of the trajectories has also been performed by Höglund & Ågren (2001) for a single-phase diffusion couple and by Larsson *et al.* (2006) for a couple with a moving interface. Sharp and diffuse interface computational methods will be used in the present paper. In particular, the sharp interface motion will be described using two methods: similarity (Boltzmann–Matano) solution for an infinitely long diffusion couple and a discrete method for a finite-length diffusion couple. The Dictra code will be used for the discrete calculations. A discrete method will also be used for the diffuse interface calculations, again for a finite-length diffusion couple. The FiPy code (Guyer *et al.* 2007) will be used for the diffuse interface calculations. The use of two models, sharp and diffuse, and three computational methods will shed light on the details of the deformation near the moving α/β interface and in situations where the deformation is multi-valued, i.e. where multiple Kirkendall planes exist. Possibilities of stress generation due to diffusion in constrained geometries will also be discussed.

## 2. Models

### (a) Diffusion and marker velocity

To permit application to both sharp and diffuse interface approaches, we seek expressions for the standard Darken (1948) approach in terms of gradients of the chemical potentials, *μ*_{A} and *μ*_{B}, of the components A and B. Conservation laws for the concentrations *c*_{A} and *c*_{B} (moles per unit volume) have the form(2.1)Here and are the fluxes of components A and B, respectively, measured in the laboratory frame. It is customary in treatments of the Kirkendall effect to define another pair of fluxes and with respect to a velocity field and posit constitutive relations for them. The velocity field is measured in the laboratory frame and characterizes the motion of a hypothetical set of markers embedded in the material. If the material is crystalline, these markers move with the lattice, which in the Darken treatment is defective due to vacancy sources and sinks. The fluxes and are called intrinsic (or lattice) fluxes and are related to those measured in the laboratory frame by(2.2)

The conservation laws (2.1) then have the form(2.3)In a binary system, the concentrations, *c*_{A} and *c*_{B}, and mole fractions, *X*_{A} and *X*_{B}, satisfy the relations (Lupis 1983)where and are the partial molar volumes of species A and B, respectively; and *V*_{m} is the molar volume. These relations hold only when the partial molar volumes are constant. Variable partial molar volumes are treated by Philibert (1991) and Sekerka (2004) and complicate the analysis considerably.

Multiplying the conservation law for each species in equation (2.3) by its corresponding partial molar volume, and adding, gives(2.4)If the fluxes and velocity are assumed to be only in the *z*-direction, this equation can be integrated to give(2.5)where the velocity and fluxes are now scalars and we have assumed that both decay in the far field.

If this expression for *v* is used in the one-dimensional version of equation (2.3) for species B, we obtain(2.6)The intrinsic fluxes are assumed to be proportional to the gradients of their chemical potentials, which in one dimension are(2.7)where the intrinsic (or lattice) mobilities *M*_{A} and *M*_{B} can depend on the concentration.

Making use of the Gibbs–Duhem equation, the marker velocity field is then given by(2.8)and the diffusion equation is given by(2.9)where the interdiffusion mobility is defined by(2.10)It is now clear that(2.11)

As is well known, marker motion occurs due to different mobilities and different partial molar volumes, as seen in equation (2.8). For the rest of this paper, we will assume that and thus marker motion will be due solely to the difference between *M*_{A} and *M*_{B}. Thus, with , we arrive at the equations used in this paper; *viz*. the marker velocity field is given by(2.12)and the diffusion equation becomes(2.13)where, now, the interdiffusion fluxes are given by(2.14)

More generally, we note that it is only the divergence of the velocity field which is determined by the solute conservation laws. If the velocity field is not one-dimensional, additional equations are required to determine the flow (Boettinger *et al.* 2005). We also note that equation (2.12) is identical to that used by Wu *et al.* (2004) when specialized to the case of a binary alloy. This equation is correct only when the velocity is constrained to be in the diffusion direction. The more general case is described in appendix B.

This velocity field characterizes the deformation of the material during diffusion. If one envisions diffusion in a crystalline material, the markers are (coarse-grained) positions in an imperfect lattice that is capable of emitting or absorbing vacancies at sufficient rates to maintain an equilibrium chemical potential of vacancies everywhere. The same equations also apply to diffusion in fluids or amorphous solids. Here the velocity describes the net convective motion of the substance as indicated by the motion of small inert particles in a fluid even though constituent components A and B have different diffusion mobilities. Indeed, Meyer's formula, describing the relation of an effective diffusion coefficient to the independent diffusivities of components, in the manner of Darken (1948), and describing the net convective motion, was developed in 1899 in the kinetic theory of gases (e.g. Kennard 1938, pp. 184–190).

### (b) Chemical potentials

The sharp and diffuse interface approaches to the moving interface problem require two different forms for the chemical potential difference. For the diffuse interface treatment, we employ a non-classical chemical potential that includes a gradient energy term with coefficient *K*, as in the approach of Cahn & Hilliard (1958),(2.15)where are the chemical potentials for a material with no gradients. We note that Langer & Sekerka (1975) employed a diffuse interface Cahn–Hilliard model to compute the composition profile in a diffusion couple with a moving interphase interface to study deviations from local equilibrium at the moving interface. For the sharp interface model, we let *K*=0. For the solution thermodynamics, we use a regular solution, such that(2.16)The regular solution parameter *Ω* is taken positive, *R* is the gas constant and *T* is the absolute temperature. Substitution of equation (2.16) into equations (2.12) and (2.13) gives the Kirkendall velocity and diffusion equation for the diffuse interface model.

### (c) Diffusion coefficients

For *K*=0, we recover the usual diffusion and Kirkendall velocity equations. The intrinsic diffusion coefficients *D*_{A} and *D*_{B} are defined by(2.17)The interdiffusion coefficient is defined by the relation(2.18)and, using equations (2.10) and (2.14), is given by(2.19)The term in brackets is the thermodynamic factor for the regular solution. The intrinsic diffusion coefficients are then(2.20)The diffusion equation (2.13) for *K*=0 becomes(2.21)and the marker velocity, equation (2.12), becomes(2.22)

## 3. Parameters

The materials parameters used for the calculations are shown in table 1. The partial molar volumes are taken as equal. The positive regular solution parameter *Ω* leads to a miscibility gap. At 1000 K, the composition limits of the two phases are and (figure 1). The use of a single free-energy function limits the diffuse interface analysis to two phases with different compositions, but the same crystal structure. The value of *K* was chosen such that the interface width has a reasonable value.

The intrinsic mobilities have been selected as linear functions of mole fraction. Two sets of intrinsic mobilities are used and denoted as *expanding* and *squeezing*, referring to whether the deformation of the diffusion couple is towards the ends of the diffusion couple or towards the centre. Figure 2 shows the mobilities for the ‘expanding case’ and figure 3 shows the diffusion coefficients derived from these mobilities and the regular solution model. The functions *M*_{A}(*X*_{B}) and *M*_{B}(*X*_{B}) have been constructed to make symmetric around *X*_{B}=0.5 and (*D*_{B}−*D*_{A}) antisymmetric around *X*_{B}=0.5 to elucidate the effects described in §1. For the expanding case, species A diffuses faster in the A-rich phase and species B diffuses faster in the B-rich phase. The reverse is true for the squeezing case.

## 4. Computational methods

### (a) Similarity solution

#### (i) Solute field

For an infinitely long diffusion couple, one can solve equation (2.21) using the Boltzmann–Matano transformation (Matano 1933) to compute a similarity solution for the solute field. This is possible for initial conditions of the form(4.1)where and are constants, with the boundary conditions(4.2)

The values of *X*_{B} are fixed on either side of the moving α/β interface *z*_{I}(*t*), according to(4.3)The flux continuity condition is given by(4.4)where *v*_{I} is the interface velocity.

Let denote a representative value of , to be used for non-dimensionalization. The diffusion equation admits a similarity solution of the form(4.5)The resulting ordinary differential equation (ODE) for is then(4.6)where the subscript ξ means differentiation with respect to ξ and where , with the far-field boundary conditions(4.7)

The similarity solution then satisfies(4.8)(4.9)where the phase boundary position has the form .

#### (ii) Velocity field

Equation (2.22) can be rewritten in terms of the similarity variable as(4.10)where(4.11)

#### (iii) Displacement field

Since the velocity field *v*(*z*, *t*) is described in terms of laboratory coordinates, the displacement field *u*(*z*, *t*) is also taken to be a function of the laboratory coordinates and time. The relation between the reference and laboratory coordinates is defined by the mapping *Z*=*Z*(*z*, *t*), giving the reference (initial) position of the marker that is currently at position *z*; the inverse mapping is *z*=*z*(*Z*, *t*), giving the current position *z* of the marker that was at position *Z* initially. The displacement field *u*(*z*, *t*) then satisfies(4.12)That is, *u*(*z*, *t*) is the difference between the current and initial positions of the marker that is presently at position *z* in the laboratory frame. We note that in most treatments of solids, the displacement field is instead considered to be a function of the reference position *Z* and time, which corresponds to the function *U*(*Z*, *t*)=*u*(*z*(*Z*, *t*), *t*) in our notation; the two approaches are equivalent as long as the mapping remains invertible.

The displacement and velocity fields satisfy the equation(4.13)We set(4.14)giving(4.15)with decay conditions on as *ξ*→±∞. This equation has a singularity at points where(4.16)and at those points, we also have , so that there. Since equation (4.12) gives(4.17)points where correspond to the Kirkendall plane *Z*=0 in the reference frame.

In fact, note that the function for all values of *ξ* satisfies equation (4.15) identically, and for this displacement each point *z* in the laboratory frame is mapped to the single point *Z*=0 in the reference frame. Also note that for ,(4.18)so that the strain *∂u*/*∂z* is not small in this case. The mapping of multiple points to a single point does not occur if the strain is small. We will make use of the solution in the case of multiple Kirkendall planes, as described in appendix C (electronic supplementary material).

With a different notation, equation (4.16) and the location of the singularity were obtained by Cornet & Calais (1972). Equation (4.16) converted to dimensional form gives the condition at the singularity as(4.19)an expression widely used in the experimental papers cited above. At a fixed time, the intersection(s) of a plot of the velocity profile *v*(*z*, *t*) and the line *v*=*z*/2*t* determine the value(s) of *z*, where Kirkendall markers placed on the initial interface are located after diffusion. Note that this construction assumes that the solute and velocity profiles have the self-similar forms given by equations (4.5) and (4.14).

With a self-similar form for the solute, velocity and displacement fields, the solution can be computed by solving ODEs rather than partial differential equations (PDEs). We employ a nonlinear shooting method based on an adaptive ODE solver in order to obtain solutions with a high degree of accuracy; details are given in appendix C (electronic supplementary material).

### (b) Dictra solution

Dictra (Ågren 1982; Andersson *et al.* 2002) is a code to perform one-dimensional diffusion simulations using the finite-difference method. It requires the existence of a Calphad (Kaufman & Bernstein 1970; Saunders & Miodownik 1998) type thermodynamic database for the phase(s) of interest to compute the thermodynamic factors of diffusion and a mobility database. Such a database can be constructed using the optimizer within Dictra to fit the various experimentally measured diffusion coefficients.

For the present paper, the thermodynamic description given in table 1 is treated in Dictra with two identical phase descriptions, entered as α and β, where α and β are the A-rich and B-rich phases, respectively. The positive regular solution parameter *Ω* is entered as a zero-order binary interaction term in the free-energy description. However, the linear composition-dependent mobility descriptions given in table 1 cannot be directly used in Dictra because the mobility is expressed as an exponential of a polynomial function of the mole fraction. We approximate the linear dependence of the mobility of each component by the function available within Dictra as(4.20)where are pure component terms and are zero-ordered binary interactions. The pure component terms are fixed by the end-member values of the linear functions in table 1. The values for are determined to best fit the linear behaviour at 1000 K. Values used for the expanding and squeezing cases are as follows:

The Dictra simulations are performed over the sample length 2*L* with no flux end conditions and the initial conditions described in §5. A 300-point geometric grid is employed, placing a higher density of points at the α/β interface. Starting from the left end of the sample with an initial grid spacing of *a*_{0}, each successive grid spacing is set to *a*_{n}=*a*_{n−1}*r*, where *r*=0.98 to the left of the initial interphase interface and *r*=1.02 to the right of it. During the simulations, an adaptive grid mesh is used, which results in the automatic addition and deletion of grid points as the interface moves. The maximum time step allowed is limited to two orders of magnitude less than the total time, to maintain sufficient accuracy for the mesh points in the immediate vicinity of the phase boundary. The diffusion equation is integrated implicitly. Outputs are the spatial profiles of *X*_{B} and the vacancy flux for selected times. The Kirkendall velocity is the vacancy flux times the molar volume.

Displacements were not computed from the Dictra-computed velocities. Such calculations must be performed outside of the Dictra computational environment, as described by Höglund & Ågren (2001).

### (c) Diffuse interface solution

With a diffuse interface approach, the interface width introduces an additional length scale, which breaks the scale invariance of the problem. Since the solution does not have a self-similar form, we must compute a numerical solution of the full PDE. Equation (2.13), with the substitution of equations (2.15) and (2.16), becomes(4.21)This equation is solved for *X*_{B}(*z*, *t*) with the FiPy finite volume package (Guyer *et al*. 2007); an example script is included in the electronic supplementary material. The computational domain consists of 2000 uniformly distributed mesh points over −*L*≤*z*≤*L*, with initial conditions(4.22)

The boundary conditions ∂*X*_{B}/∂*z*=0 and ∂^{3}*X*_{B}/∂*z*^{3}=0 are imposed at *z*=±*L*. Equation (4.21) is discretized to second order in space and solved with fully implicit first-order time steps using matrix decomposition. An adaptive time stepper is used to keep the relative truncation error of the linear solver below 10^{−5}.

After the determination of *X*_{B}(*z*, *t*), equation (2.12) with the substitution of equations (2.15) and (2.16),(4.23)can be used to compute the velocity field. We determine the reference coordinates (corresponding to arbitrary marker positions) in terms of the laboratory coordinates by solving the PDE (obtained from equations (4.12) and (4.13))(4.24)for *Z*(*z*, *t*), subject to and . The laboratory position of a particular marker that was at *z*=*z*^{*} at *t*=0 can be determined by interpolating the value of *z*(*t*) that gives *Z*(*z*, *t*)=*Z*(*z*^{*}, 0)=*z*^{*}. Unlike in the sharp interface models, the mappings *z*=*z*(*Z*, *t*) and *Z*=*Z*(*z*, *t*) are always invertible.

## 5. Results

Two sets of initial conditions were used for numerical evaluation,(5.1)The finite-length calculations were performed on a domain −*L*≤*z*≤*L*, with *L*=10^{−4} m.

### (a) Composition profiles and interface motion

Figure 4 shows composition profiles for , for the expanding case and for *t*=8×10^{2}, 8×10^{3}, 8×10^{4} and 8×10^{5} s computed by diffuse interface and Dictra methods and at *t*=8×10^{2} s for the similarity solution. Profiles for the other three situations were obtained but are not shown. The motion of the α/β interface to the right is clearly seen. Owing to the finite-size sample used for the diffuse interface and Dictra methods, the interface slows with time reaching a final equilibrium position, representing a lever rule mixture of two phases each of uniform composition. The profiles are identical using the three methods, with the exception that the α/β interface is spread for the diffuse interface calculation. It can be seen that the full interface width is approximately 2×10^{−6} m or approximately . Simulations were also performed with interfaces of twice this thickness and half this thickness without significantly changing the results. Even the shortest time in figure 4 is long compared with the diffusion time for the diffuse interface to form . In practice, we find that the interface width reaches 80% of its steady-state value in 30 s, as determined by the distance between the values of *z*, where and .

### (b) Velocity profiles

Figure 5 shows the velocity profiles for at *t*=8×10^{2}, 8×10^{3}, 8×10^{4} and 8×10^{5} s computed only by Dictra. At this scale, the difference in results from the other methods cannot be seen. Figure 5*a* shows the expanding case. The velocity is negative (to the left) in the α-phase and positive (to the right) in the β-phase, and thus the reason the term ‘expanding’ was used. The material moves outwards from the central region of the diffusion couple towards the ends. In the context of a vacancy mechanism for diffusion, vacancies are leaving the central region and flow outwards in both phases. Figure 5*b* shows similar information for the squeezing case. In contrast to the expanding case, the velocity is positive (to the right) in the α-phase and negative (to the left) in the β-phase, and thus the term squeezing has been used. The material is moving inwards towards the central region of the diffusion couple.

Note that, in both cases, the velocities are zero at the ends of the sample and that the velocities decay with time. Had we modelled the situation where the partial molar volumes are not constant, one end would have moved with respect to the other. We also could have selected mobility functions where the markers move to the left in both phases, but at different rates, or to the right, but the most interesting cases are the expanding and squeezing cases.

Figure 6*a* shows the central region of figure 5*a* for *t*=800 s. Note that results from the sharp interface models have a discontinuous velocity profile at the moving α/β interface, whereas that for the diffuse interface model has a continuous but sharply varying velocity profile. Other than this diffuseness, the results from the three methods are in agreement. Following the analysis of §2*a*, a plot of *v*=*z*/2*t* is also shown. The intersection of this line with the velocity profile for the similarity solution indicates the *z* position of the Kirkendall plane after 800 s of diffusion.

Figure 6*b* shows similar velocity plots, also for the expanding case, but with . The velocity profiles are similar, but the plot of *v*=*z*/2*t* intersects the velocity profile for the similarity solution and the Dictra solution at three locations. This is the graphical form of the criterion established by Cornet & Calais (1972) for multiple Kirkendall planes. Indeed, the end composition, , was chosen to give multiple intersections. In figure 6*a*,*b*, the plot of *v*=*z*/2*t* intersects the velocity profile at a *z* position(s), where ∂*v*/∂*z*>0.

Figure 6*c*,*d* shows similar information for the squeezing case for and , respectively. We note again the smooth nature of the velocity profile for the diffuse interface calculation and the discontinuous profile for the sharp interface calculations. In figure 6*c* the plot of *v*=*z*/2*t* intersects the velocity profile at a *z* position, where ∂*v*/∂*z*<0. The significance of the sign of ∂*v*/∂*z*<0 will be addressed in §6.

There is no significance to the fact that the plot of *v*=*z*/2*t* intersects the diffuse interface profile only once in figure 6*b*. Other choices of *K*/*Ω* would lead to sharper interfaces and multiple intersections, but it is still possible to obtain a continuous distribution of trajectories (see §5*c*) for markers originally placed close enough to the initial discontinuity of the diffusion couple.

### (c) Displacements

The most transparent way to display the deformation that occurs in these four cases is given in figure 7. The plots, obtained with the similarity method, show trajectories (position with time) of markers that are evenly spaced at time zero. Markers placed on the original discontinuity (Kirkendall markers) are shown in dashed green and trajectories of markers placed at other initial positions are shown in red. The α/β interface position is shown in thick blue. The trajectories show the mapping between the initial (or reference) position *Z* and the actual (or laboratory) position *z* at various times. In figure 7*a*,*b*, the trajectories move to the left in the α-phase and to the right in the β-phase. In figure 7*c*,*d*, the reverse is true. As the α/β interface passes through the position of a marker, the trajectory sharply changes direction. The trajectories in figure 7*a*,*c* are straightforward but not so for figure 7*b*,*d*. Note the presence of two green trajectories in figure 7b, corresponding to two Kirkendall planes. Between the two green trajectories, there are no trajectories no matter how close together their original positions are. In figure 7*d*, trajectories end on the α/β interface.

The examination of the results using the diffuse interface model (figure 8) shows some differences from the above description. First, the sharp change in direction of markers as the α/β interfacial region passes through them is smoothed. Second, in figure 8*b*, in contrast to figure 7*b*, a distribution of markers would exist throughout the central zone if they had been placed close enough to the initial discontinuity of the diffusion couple. Third, in figure 8*d*, in contrast to figure 7*d*, the markers are not swallowed in the interfacial region but get jammed very closely together near the interface. Fourth, due to the finite-size computational domain, the long-time trajectories differ from those of figure 7. The ‘interface position’ for the diffuse interface calculation was taken as the location where the composition was for graphical presentation; it is however diffuse.

More insight into the nature of the deformation and the various singularities of the deformation can be obtained from figure 9. The plots show the marker position versus their original position after 800 s of diffusion for the four cases studied. Far from the centre of the diffusion couple (not shown), the curves approach the dashed line at 45° corresponding to no deformation. The curves are continuous and smooth for the diffuse interface model. The curves for the sharp interface model exhibit various types of singularities. The Kirkendall plane is located at the *z* position corresponding to *Z*=0. The Matano plane is always located at *z*=0 because the diffusion equation was solved in the *z*-coordinate system with the initial discontinuity at *z*=0. The α/β interface position is also indicated. The curves in figure 9*a*,*c* are continuous with discontinuities in slope at the Kirkendall interface position and at the α/β interface position. Figure 9*a*,*c* is an example of common behaviour for deformation due to the Kirkendall effect. Over the diffusion time, a region of relatively large strain, |1−∂*z*/∂*Z*|, is left behind as the α/β interface moves from its initial position at the initial discontinuity of the diffusion couple at *t*=0 (*Z*=0) to its present position.

Figure 9*b*,*d* shows unusual but interesting cases of what can happen with the Kirkendall effect. In figure 9*b*, Kirkendall markers, *Z*=0, are spread over a range of *z* positions, while in figure 9*d* there is a range of initial positions that have the same final location.

## 6. Discussion

### (a) Markers

The actual or conceptual introduction of markers in a diffusion couple serves to define a strain field that is otherwise difficult to define for a binary diffusion process, where the microscopic motion of the two species can differ significantly at each point. For binary diffusion in a fluid system, a barycentric velocity field based on the local centre of mass motion of the species is often employed. For diffusion in a crystalline solid, the motion of the lattice defines a convenient velocity field that corresponds to the marker velocity considered here. The individual motion of the components (diffusion) is referred to that lattice. For a crystalline sample, the diffusion process is often accompanied by the motion of vacancies that are created and destroyed at imperfections in the lattice, leading to the creation and elimination of lattice points. The creation or elimination of lattice points causes changes in the distance between markers and thus strain; a stress-free strain! A marker then, strictly speaking, is not attached to any one particular lattice point, and thus serves as a coarse-grained description of the ensuing motion of this defective lattice. In amorphous materials, where no lattice is present, the same definition of marker position and velocity can be applied and used to define the intrinsic diffusion fluxes of the two species (Stephenson 1988). Here the deformation effect of diffusion can only be viewed through the idea of stress-free strain.

### (b) Computational methods/discontinuity of Kirkendall velocity at a moving interface

In the models that we consider, the velocities of Kirkendall markers are determined by the concentration field, which satisfies a nonlinear diffusion equation. We have considered a one-dimensional binary diffusion couple at a constant temperature, where each member of the initial couple represents a different phase of the two-component system. For simplicity, we have assumed a constant molar volume to eliminate the effects of volume change on marker motion (Dantzig *et al.* 2006). Trajectories of the markers can be obtained from the computed velocity field, and, in particular, the motion of the Kirkendall markers initially located at the couple interface can be examined. The velocity field varies from point to point, and the velocity field in each phase can differ in direction, leading to the complex behaviour of the Kirkendall markers. In our calculations, we have employed three different numerical procedures to determine the concentration field, each of which enjoys various strengths and weaknesses.

The phase field approach always yields a solution to the composition and velocity fields that are smooth throughout space, with no discontinuities. Its use, however, does require numerical solutions of the nonlinear Cahn–Hilliard PDE in space and time. The diffuse nature of the interface generates a smooth solute profile, and thus the marker motions are given by a smooth, one-to-one mapping from the reference frame to the laboratory frame.

In contrast, the velocity field is generally discontinuous at the phase boundary for sharp interface models. In most cases, this discontinuity of velocity presents no particular problem, and the velocity can be integrated to determine an invertible mapping between the reference coordinates and the laboratory coordinates with no discontinuity in the displacement field anywhere. As illustrated in figure 9, in some cases, the associated deformation field can be discontinuous, with a non-invertible mapping from the reference frame to the laboratory frame. This can lead to the cases of multiple Kirkendall planes or Kirkendall planes coincident with the moving α/β interface.

The sharp interface model allows a similarity solution in an infinitely long diffusion couple. The similarity solution can be computed with high accuracy by using adaptive solvers for ODEs; however, the integration needs to be started and stopped at points representing the phase boundary and Kirkendall planes, where appropriate jump conditions need to be applied explicitly. Jump conditions are automatically implemented for the finite-difference method employed by Dictra. The sharp interface approach reveals detailed features near the moving α/β interface.

### (c) Vacancies

In this paper, we have not included vacancies in our treatment. They are not essential, as the Kirkendall effect can be measured in amorphous solids. However, their consideration provides a more transparent way to understand the results of this paper.

To have different fluxes of components with mass on a perfect lattice, there must also exist a flux of vacancies such that the sum of the three fluxes on the lattice (intrinsic fluxes) is zero. A flux of vacancies on a perfect lattice would lead to spatial and temporal variations of vacancy content as diffusion occurs. If one considers an imperfect lattice with sources and sinks for vacancies, the vacancy content will attempt to maintain an equilibrium level. Adda & Philibert (1966) have formally treated the situation where this attempt is inadequate, i.e. a kinetic parameter is introduced where the source/sink rate is proportional to the deviation of the vacancy concentration from the equilibrium level. In standard treatments of the Kirkendall effect, as also herein, the vacancies are assumed to maintain an equilibrium concentration. Thus, the local rate of vacancy creation/annihilation is slaved to the intrinsic fluxes of the massive components. If we let *σ*_{v} be the local vacancy source, it is determined by the divergence of the vacancy flux as(6.1)Thus, ∂*v*/∂*z* and *σ*_{v} have the same sign and vacancies are created when ∂*v*/∂*z* is positive and vacancies are annihilated when ∂*v*/∂*z* is negative. When vacancies are created, the lattice swells by the addition of lattice points and, when vacancies are annihilated, the lattice shrinks by the removal of lattice points. Most of these vacant lattice sites become populated by the atoms diffusing into the locale. We note that for situations where one chooses to ignore vacancy considerations, the quantity ∂*v*/∂*z* is the one-dimensional form of the stress-free rate of deformation tensor associated with the Kirkendall effect.

The sign of ∂*v*/∂*z* and *σ*_{v} explains one of the observations described in §1 regarding marker stability (van Dal *et al.* 2001). If one were to consider that the exact positions of markers in a row perpendicular to the diffusion direction were slightly perturbed in the diffusion direction, they may be found at later times to have spread or become concentrated in the diffusion direction. Say we are given the marker velocity field *v*=*v*(*z*, *t*), expressed as a function of the laboratory coordinate *z* and time *t*. Then, the marker that is originally at the point *Z* in the reference frame has the position *z*=*z*(*Z*, *t*) in the laboratory frame, which satisfies(6.2)Consider the motion of two nearby markers, one initially at the point *Z* and the other at the point *Z*+Δ*Z*. Their current positions in the laboratory frame are at *z*(*Z*, *t*) and *z*(*Z*+Δ*Z*, *t*), and their current separation in the laboratory frame is(6.3)To leading order, the rate of change of their separation Δ*z* is then(6.4)where we have used the leading order approximation of equation (6.3) in the last step. Thus, we have(6.5)which says that the relative rate of change of the separation of nearby markers is given by the velocity gradient; if ∂*v*/∂*z*>0, the markers separate in time.

In a sharp interface model, we have seen that a discontinuity of Kirkendall velocity usually occurs at the moving α/β interface. In the context of the vacancy mechanism, this requires an interface source term related to the jump in velocity across the interface; *viz*. the interface vacancy source is , where is the marker velocity in phase ϕ, evaluated at the interface.

In the expanding case, the jump in Kirkendall velocity across the α/β interface is positive and the interface must serve as a source of vacancies and lattice points. This is the reason why the sample does not physically tear apart with velocities moving away from the central region on both sides. Conversely, in the squeezing case, the jump in Kirkendall velocity across the α/β interface is negative and the interface must serve as a sink of vacancies and lattice points.

### (d) Multiple Kirkendall planes/non-invertible mappings

Two of the four cases that we consider here feature instances in which the mapping between reference coordinates and spatial coordinates breaks down. In the expanding case of marker motion directed outwards from the phase boundary, a denuded zone between oppositely directed planes can appear, where a region in the laboratory frame is not covered by points from the reference frame. Conversely, in the squeezing case of marker motion directed towards the phase boundary, markers may accumulate at the interface, so that a finite region in the reference frame is mapped to a single point in the laboratory frame. These examples illustrate the complications that can occur in applying the formalism of continuum mechanics to problems of multicomponent, multiphase diffusion.

In the context of the vacancy model, one can see that when large numbers of lattice sites are inserted or removed, there is a breakdown in the usual invertible mappings of continuum mechanics from initial position to actual position of markers. When lattice sites are inserted after the beginning of the diffusion process, there will be material at a range of positions that appear to have had the same initial position (figure 7*b*). On the contrary, when lattice sites are removed, there will be a range of initial positions that correspond to a single spatial position (figure 7*d*).

### (e) Stress

In these models, the stress field that accompanies the marker-defined strain field is decoupled from the determination of the concentration field, so that the computed solute and displacement profiles are independent of any underlying stress–strain model. Given such a model, however, a stress field that is consistent with the solute and strain distributions can be computed. Since we assume a one-dimensional displacement field, the associated stress field generally requires tractions normal to the diffusion direction, which constrain the displacement field to the diffusion direction. The component of traction along the diffusion direction can be assumed to vanish. The resulting stress field is rather artificial, being designed to accommodate a one-dimensional displacement field in the couple. In practice, for samples with traction-free boundaries, the Kirkendall effect often leads to lateral bulging or contraction of the sample surface in the diffusion zone (Boettinger *et al.* 2005). In the interior of a thick sample (lateral dimension ), the displacement occurs primarily in the diffusion direction and has the form considered here and the transverse stress is as computed in appendix B.

## Acknowledgments

We thank S. R. Coriell, R. F. Sekerka and J. A. Warren for their helpful discussions.

## Appendix A Matano interface

Diffusion coefficients are often measured by performing a Boltzmann–Matano analysis (Matano 1933) of the solute profile in a diffusion couple. Given an experimental diffusion profile, the analysis often starts with the determination of the so-called Matano plane *z*_{M}, which is defined by the equation(A1)which represents an equal-area construction for the solute profile on either side of *z*_{M}. Here we are assuming an effectively infinite couple with to simplify the discussion. For the two-phase couple, the discontinuity of the integrand at the α/β interface must be taken into account in performing the integration. With the initial discontinuity located at *z*=0, the Matano plane is initially given by *z*_{M}=0 as well, and it is a property of the diffusion equation that the Matano interface actually remains at *z*_{M}=0 for subsequent times as long as the gradient vanishes at the ends of the interval. This can be verified (e.g. Sekerka 2004) by differentiating both sides of equation (A 1) with respect to time, replacing ∂*x*_{B}/∂*t* by , and performing the integration. Here the moving phase boundary generates additional terms upon differentiation that are seen to cancel by invoking the jump condition at the α/β interface. The vertical dotted line in figure 8 confirms that the empirically determined Matano interface remains at *z*_{M}=0 in the finite couple as well.

## Appendix B Stress due to diffusion in constrained geometry

An idealized geometry for one-dimensional displacements of the type we consider consists of an infinite cylinder with an axis along the *z*-direction. In a cylindrical coordinate system (*r*, *θ*, *z*), appropriate boundary conditions on the side of the cylinder (*r*=const.) are(B1)where is the normal vector to the side of the cylinder; is the surface displacement; and and are the tangential surface tractions on the sides of the cylinder; here *σ* is the Cauchy stress tensor and and are unit vectors in the *θ*- and *z*-directions, respectively. The vanishing tangential tractions allow displacement along the axial direction, while a non-zero normal component of the surface traction is required to constrain the displacement to the *z*-direction, so that . With these boundary conditions, *σ*_{rr}(*z*, *t*), *σ*_{zz}(*z*, *t*) and *σ*_{θθ}(*z*, *t*) are the principal components of stress, with *σ*_{rr}=*σ*_{θθ}. The stress equilibrium equations ∂*σ*_{jk}/∂*x*_{k}=0 then have a single non-trivial component, requiring ∂*σ*_{zz}/∂*z*=0, and the stress-free conditions that prevail as *z*→±∞ then require *σ*_{zz}=0 as well.

The expression for *v*(*z*, *t*) given in equation (2.22) is a special form of a more general three-dimensional equation (Boettinger *et al.* 2005),(B2)where is the traceless deviatoric stress and *E* and *μ* are the elastic and viscous moduli, respectively. The l.h.s. term is the rate of deformation tensor (Malvern 1969). Computations based on a similar model have been performed by Szabó *et al.* (2001). If is unidirectional, then these equations take the formOn summing, it gives(B3)which can be integrated to obtain equation (2.22). Thus, this equation is also consistent with the stress–strain relation (B 2).

For our situation, the stress components have *σ*_{33}=0, *σ*_{11}=*σ*_{22}=*σ*_{rr}, so that *σ*_{11}+*σ*_{22}+*σ*_{33}=2*σ*_{rr}. For simplicity, we set 1/*E*=0 and neglect elastic strains. Then, we have(B4)In terms of the dimensionless similarity solution, this gives(B5)An example for the squeezing case with is shown in figure 10. Both tension and compression are required on the sides of the cylinder *r*=const. in the diffusion zone in order to obtain a one-dimensional displacement along the axis. Note that the discontinuous slope of the velocity profile leads to a discontinuity of *σ*_{rr} at the phase boundary. The stress balance at the interface is still satisfied, however, since the traction on the surface *z*=const. only involves the stress components *σ*_{jz}, which vanish identically in this geometry.

## Footnotes

Electronic supplementary material is available at http://dx.doi.org/10.1098/rspa.2007.1904 or via http://www.journals.royalsoc.ac.uk.

↵Certain commercial equipment, instruments or materials are identified in this paper in order to specify the experimental procedure adequately. Such identification is not intended to imply recommendation or endorsement by the National Institute of Standards and Technology, nor is it intended to imply that the materials or equipment identified are necessarily the best available for the purpose.

- Received March 29, 2007.
- Accepted September 11, 2007.

- © 2007 The Royal Society