## Abstract

Often during phase growth, the rate of accretion, on the one hand, is determined by a competition between bulk diffusion and surface reaction rate. The morphology of the phase interface, on the other hand, is determined by an interplay between surface diffusivity and surface reaction rate. In this study, a framework to predict the growth and the morphology of an interface by modelling the interplay between bulk diffusion, surface reaction rate and surface diffusion is developed. The framework is demonstrated using the example of Cu–Sn intermetallic compound growth that is of significance to modern microelectronic assemblies. In particular, the dynamics and stability of the interface created when Cu and Sn react to form the compound Cu_{6}Sn_{5} is explored. Prior experimental observations of the Cu_{6}Sn_{5}–Sn interface have shown it to possess either a scalloped, flat or needle-shaped morphology. Diffuse interface simulations are carried out to elucidate the mechanism behind the interface formation. The computational model accounts for the bulk diffusion of Cu through the intermetallic compound, reaction at the interface to form Cu_{6}Sn_{5}, surface diffusion of Cu_{6}Sn_{5} along the interface and the influence of the electric current density in accelerating the bulk diffusion of Cu. A stability analysis is performed to identify the conditions under which the interface evolves into a flat, scalloped or needle-shaped structure.

## 1. Introduction

Many physical phenomena involving phase growth are driven by an interplay between bulk-diffusion of one or more reactant species to the phase interface, reaction rates at the interface as well as the surface diffusion of the product along the phase interface. Examples of such phenomena include dendritic growth in Li-ion batteries [1,2], growth of brittle intermetallic compounds in solder joints [3], tissue growth [4–6] and growth of compounds under chemical-vapour deposition [7]. The interplay between the aforementioned driving forces can not only affect the growth rate of the materials involved, but can also change the interface morphology from being flat, to scalloped to even dendritic. The growth of these materials can have significant physical as well as physiological effect. For example, the growth of dendrites on the negative electrode in the Li-ion battery can make it unsafe for commercial use [1] and moderating the growth of malignant tissues can aid in recovery of diseased individuals [6]. Recently, cracking of micrometre-scale solder joints has been attributed as the cause for an aircraft crash [8]. Another example of morphological instability is the formation of dendrites during synthesis of graphene on a copper substrate using chemical vapour deposition [7]. Tailoring the morphology of graphene is crucial to realize its special electronic properties. Yet another example emphasizing the importance of morphological instabilities is electroplating. Surface roughening of the plated metal depends on the current density, and higher current densities lead to dendritic or nodular deposits [9]. Thus, there is a significant need at the present time for predicting the growth and morphology of an interface under a given set of environmental conditions.

The competition between surface growth rate and surface diffusion has been analytically studied by Mullins & Sekerka [10], who derived an interface stability criterion to determine the morphology of the interface. This stability analysis has been developed further by other researchers for studying the solidification problem [11,12]. In general, to study the role of bulk diffusion, computational models are necessary as complex geometry and boundary conditions make it difficult to predict species concentrations analytically. At the present time, a computational framework to study the interplay between bulk diffusion, interface reaction and surface diffusion does not appear to have been demonstrated in the prior literature. The demonstration of such a modelling framework using the example of the growth of the Cu_{6}Sn_{5}–Sn interface is the goal of this study.

The dynamics of the formation of intermetallic compounds (IMC) and the stability of their interfaces have significant practical relevance. For instance, the mechanical behaviour of the ubiquitous solder joints in microelectronic assemblies are determined by the shape and the extent of Cu–Sn intermetallic compounds. Intermetallic compounds are generally stiff and brittle [13,14]. As a result, their mechanical behaviour is at odds with the role of solder joints in providing compliant interconnects that absorb the differential thermal expansion between the silicon chip and the polymeric substrate. In general, thermo-mechanical fatigue cracks grow very close to the interface with the IMC, but in the bulk solder material. At higher strain rates, cracks typically propagate through the IMC. Thus, both the thickness and the geometrical shape of the IMC have bearing on the reliability of the joint. The growth of the intermetallics is accelerated by the extent of current density [15,16], and smaller joints experience higher current densities making them even more susceptible to failure by excessive IMC growth.

The IMC in the Cu–Sn system consists mostly of Cu_{6}Sn_{5} and partly Cu_{3}Sn [17,18]. Cu_{6}Sn_{5} starts forming first and grows at a faster rate when compared with Cu_{3}Sn. The rate of bulk diffusion of Cu through the Cu_{6}Sn_{5} is believed to be significantly lower than the rate at which reactions at Cu_{6}Sn_{5}–Sn interface proceed [18,19]. Experimentally, the rate of Cu_{6}Sn_{5} growth has been found to be sub-linear with respect to time, strongly indicating that the process is diffusion limited [20]. The growth rate is proportional to *t*^{1/3} or *t*^{1/2} depending on whether grain-boundary diffusion or bulk diffusion is dominant [21]. Also, the Cu_{6}Sn_{5} interface has been experimentally observed to take a variety of shapes varying from flat to scalloped or even rod/needle-shaped [17–19]. These observed morphologies are shown in figure 1. When the solder joint is in a molten state, the Cu_{6}Sn_{5} layer is known to grow as scallops instantaneously. However, after solidification and ageing, the Cu_{6}Sn_{5} interface relaxes into a flat morphology. The rod- or plate-shaped crystals appear to form from all grains and do not seem to be influenced by grain anisotropy. The number and length of the rod or plate-shaped crystals is dependent on the rate of cooling during the assembly operation [17]. In general, slower cooling leads to larger number of rod- or plate-shaped structures. Further complicating matters, as the technological trends drive solder joints to smaller size, it is not uncommon for the intermetallic compounds to occupy 20–40% by volume of the joint.

In general, an interface growth simulation must account for bulk diffusion of Cu through the Cu_{6}Sn_{5} intermetallic compound, and the rate at which reactions occur at the interface, in addition to surface energy of the interface and the curvature-dependent diffusion along the interface. A simulation taking into account all of the above-mentioned mechanisms for the conditions leading to the various morphologies of the Cu_{6}Sn_{5} compound is largely missing in the literature. Zhou & Johnson [22] and Park *et al.* [23] appear to be the only ones attempting to numerically model the growth of the Cu_{6}Sn_{5} intermetallic compound. But, a study of the stability of the interface leading to the varied interface morphologies is largely absent in these studies.

In this study, a model for the growth of the Cu_{6}Sn_{5} compound based on an interface curvature-dependent reaction rate is used. A linear stability analysis is performed to determine the parameter that governs the shape of the interface. Next, the interplay between surface reaction, surface diffusion and bulk diffusion for the formation of Cu_{6}Sn_{5} is captured numerically using a diffuse interface model [24]. As, in general, physically realistic simulations are critically dependent on the availability of accurate material properties, and as material properties such as the bulk diffusivity of Cu through Cu_{6}Sn_{5}, surface reaction rate, surface diffusivity of the Cu_{6}Sn_{5} molecule are not readily available, the approach used in this study is to identify the combinations of the material properties that would lead to the formation of experimentally observed morphologies, while including in the simulations all the physically important mechanisms. Finally, using the developed numerical code, the effects of surface reaction, surface diffusion, bulk diffusion and current density on the size and morphology of the interface are studied. Conclusions are drawn on the shape of the interface and its evolution.

## 2. A model for Cu_{6}Sn_{5} intermetallic compound growth

The reaction between Cu and Sn that forms over 95% of the modern solder alloys results in the formation of Cu_{6}Sn_{5} and Cu_{3}Sn. Experiments and thermodynamic analysis of the Cu–Sn phase diagram indicates that precipitation of Cu_{6}Sn_{5} occurs earlier than Cu_{3}Sn [23]. The bulk of the intermetallic layer is known to be Cu_{6}Sn_{5}, and it is the morphology of this interface that is studied here. The growth of Cu_{3}Sn is not studied in order to simplify the model.

For the reaction leading to the formation of Cu_{6}Sn_{5} to occur, Cu must diffuse through the IMC. The diffusion of copper in the IMC is governed by its chemical potential, which is a function of the concentration of copper as well as the stress at a material point. Additionally, the diffusion of Cu can be accelerated by the presence of high current density (electromigration) or high thermal gradients (thermomigration). The thickness of Cu_{6}Sn_{5} can potentially decrease during the assembly process due to dissolution of the compound into the molten solder; however, this is not considered in the developed model.

It is proposed here that the observed morphology of the Cu_{6}Sn_{5}–Sn interface is due to competition between curvature-dependent reaction rate and surface diffusion of Cu_{6}Sn_{5} during the formation and growth of the compound [25,26]. Figure 2 summarizes the essential features of the model for the intermetallic compound growth. The proposed mechanism consists first of the inter-diffusion of Cu atoms through the Cu_{6}Sn_{5} intermetallic compound to the interface followed by reaction with Sn at the interface to form the intermetallic. This compound is then diffused along the interface driven by the surface chemical potential gradient. In a critical difference with other studies in the literature that have attempted to model the growth of Cu_{6}Sn_{5}, it is argued that the curvature of the interface affects the chemical potential of the interfacial copper and thereby alters the equilibrium constant and the kinetics of the reaction. Additionally, the specific values of diffusivity and reaction constants, whose extreme values depend on the temperature and on whether the reacting Sn is in liquid or solid state, are also proposed to determine the interface morphologies. The governing equations of the model are developed in the following subsections.

### (a) Electric potential and current density

As solder and Cu_{6}Sn_{5} are electrical conductors, there is no charge accumulation at any point in the bulk of the joint. Therefore, the electrical current or charge flux is divergence free:
**J**_{e} denotes the current density. The current is driven by gradients in the electrical potential. The current density is related to the gradient in the electrical potential, *ϕ*_{e} through the electrical conductivity *σ*:
*Γ*_{ϕe}, and Neumann boundary condition on total electrical flux, _{6}Sn_{5}–Sn interface, thus,
*Γ* is the Cu_{6}Sn_{5}–solder interface.

### (b) Bulk diffusion of Cu in Cu_{6}Sn_{5}

Potentially, three possible diffusion paths are present in the system: the inter-diffusion of Cu though the Cu_{6}Sn_{5}, the self-diffusion of Cu_{6}Sn_{5} and the self-diffusion of Sn. Of these, only the bulk diffusion of Cu through the Cu_{6}Sn_{5} is modelled here. Thus, the diffusion problem is solved only in the Cu_{6}Sn_{5} phase.

For diffusion, only the gradients in the chemical potential are important, and hence the chemical potential only needs to be determined with respect to an arbitrary datum. The datum for the chemical potential of Cu in Cu_{6}Sn_{5} is typically chosen to be that in bulk compound at equilibrium concentration, when not being subject to stress. With respect to this datum, the variation of the chemical potential with Cu concentration and stress is given by [27]
*μ*_{ΩCu} is the chemical potential of Cu in the bulk of Cu_{6}Sn_{5}, *ρ*_{Cu} is the concentration of Cu, *ρ*^{eq}_{Cu} is the equilibrium concentration of Cu, *p* is the hydrostatic stress and *ν* is a coefficient related to the difference in the atomic volume between the Cu atoms and Cu_{6}Sn_{5} molecules.

The chemical potential gradient is not the only cause for the diffusive flux; the current density can also cause a diffusive flux through electromigration. The driving force for electromigration is commonly expressed as [28]
*Z** is a parameter that relates the diffusive flux to the applied electrical current determined experimentally or through *ab initio* models [29,30].

Taking the gradient of the chemical potential in equation (2.6) and adding to the electromigration driving force of equation (2.7), the total driving force for diffusion is obtained as
*D*/*kT*)*ρ*_{Cu}, where *D* is the interdiffusivity of Cu in Cu_{6}Sn_{5}.

Finally, the species mass balance in the bulk for Cu yields
*ρ*_{Cu} is obtained using the above bulk diffusion equation with the appropriate boundary conditions, which are discussed next.

### (c) Boundary conditions at the interface

At the interface, the mass balance relation must be written for both Cu and Cu_{6}Sn_{5}. The interfacial mass balance equation can be derived by considering the mass balance on an interfacial pill-box control volume and then localizing the equations [24]. The coupled relations defining the normal velocity of the interface are

The Cu atoms from the bulk enter/leave the interface at a rate of **J**_{Cu}⋅**n**. At the interface, Cu reacts with the Sn to form Cu_{6}Sn_{5}, and the rate of production of Cu_{6}Sn_{5} is *r* moles m^{−2} s^{−1}. For every molecule of Cu_{6}Sn_{5} that is formed, six atoms of Cu are consumed. As the Cu_{6}Sn_{5} forms, the interface advances at a velocity whose normal component is *V*_{n}. **h**_{ΓIMC} is the surface flux of the Cu_{6}Sn_{5} molecules driven by the gradients in the surface chemical potential of Cu_{6}Sn_{5} given by [33]
*W* is the strain energy density, *γ* is the surface energy coefficient for the Cu_{6}Sn_{5} interface with solder and *Ω*_{IMC} is the volume of one molecule of Cu_{6}Sn_{5}. *κ* is the curvature d^{2}*z*/d*x*^{2}, defined with respect to the coordinate system shown in figure 2.

Thus, the gradients in the surface chemical potential come about due to gradients in the surface curvature *κ* as well as the gradients in the strain energy *W*. The surface flux, **h**_{ΓIMC}, is driven by gradients in the surface chemical potential [34]:
*M*_{Γ} is the interface mobility.

Reactions defined by equations (2.11) and (2.12) relate *ρ*_{IMC}, *ρ*_{Cu}, *J*_{Cu}, *V*_{n} and *r*. To reduce these unknowns further, additional assumptions are necessary, which are obtained by considering the kinetics of the reactions. Broadly, two important reactions occur at the interface:
*et al.* [25] is used here for the second reaction given in equation (2.16):
*Ω*_{Cu} is the volume of a Cu atom and *γ* is the interface energy. The above equation results when one assumes the reaction to be first order, making *r* proportional to *ρ*_{Cu}. Further, the rate of the reverse reaction is assumed to be small compared with the rate of the forward reaction. The exponential term accounts for the fact that the reaction rate depends on the chemical potential of the reactant Cu [27]. As the chemical potential of Cu at the interface varies with curvature, the equilibrium constant and the kinetics of the reaction become dependent on the curvature thereby altering the rate of the reaction. At a location where the interface curvature is negative, the chemical potential of Cu is higher as described by equation (2.13), increasing the rate of formation of Cu_{6}Sn_{5} and vice versa.

Thus, substituting equations (2.14) and (2.17) into equations (2.11) and (2.12), we get

Although the above two equations are still dependent on the earlier mentioned unknowns, they can be further resolved by considering the relative rate of the two reactions. In general, the slower reaction becomes the rate determining step.

At equilibrium, when *V*_{n}=0, one expects from equation (2.18) the condition that *ρ*_{Cu}≪*ρ*_{IMC} on the interface is made here to approximate at any time *ρ*_{Cu}*V*_{n}≈0, leading to the following condition:

An alternative but simpler form of the above boundary condition is possible when the system is reaction limited. The rate of reaction given in equation (2.15) may be modelled using a form that is frequently used in the literature to describe the rate at which bulk atoms diffuse into the interface [35,36]:
*A* is a constant proportional to the bulk diffusivity and *μ*_{ΓCu} is the surface chemical potential of a Cu atom. In general, when λ≪*A*, the system is reaction limited and when λ≫*A*, the system is diffusion limited. When the system is reaction limited, using the constitutive relation of equation (2.21) in equation (2.20) yields
*μ*_{ΩCu}≈*μ*_{ΓCu}=−*Ω*_{Cu}*γκ*, equation (2.6) suggests that the limiting concentration at the interface will be of the form

Finally, exploiting the fact that *ρ*_{Cu}≪*ρ*_{IMC}, the number density *ρ*_{IMC}≈1/*Ω*_{IMC}, and given the atomic volume *Ω*_{IMC}, equation (2.19) can be rewritten as an expression for the interface velocity, *V*_{n}:

## 3. Stability of the interface

In many physical problems, it is observed that the interface between two phases evolve either as dendrites, cellular structure or as a flat surface. Mullins & Sekerka [10] developed a method to study the stability of the interface. In this section, a similar approach is followed to determine the dimensionless parameter that governs the stability of the interface of Cu_{6}Sn_{5}.

The expression for the velocity of the interface was derived above in equation (2.24). From the expression, we can observe that if the interface was perfectly flat to begin with, the curvature *κ* will be zero, and, therefore, the interface will continue to remain planar. In order to study the stability of the interface, we assume that it is initially perturbed by an infinitesimal sinusoidal wave of the form
*δ* is assumed to be small, and the coordinate system shown in figure 2 is used. If the perturbation decays due to the evolution law given in equation (2.24), the interface is said to be stable. On the other hand, if the perturbations grow with time, the interface is said to be unstable.

For an interface defined by a curve *z*(*x*), the curvature of the perturbed interfaces is
*M*_{Γ} to be constant along the interface, the following expression for *V*_{n} is obtained from equation (2.24):
*x*, an expression for *M*_{Γ}=*D*_{Γ}*h*/*kT*, where *D*_{Γ} is the surface diffusivity and *h* is the interface thickness. Thus, the above equation can be rewritten as

If

From the above expression, one can identify the following dimensionless parameter that determines the stability of the interface:
*ζ*<1, the interface is stable and vice versa. It is seen that when *D*_{Γ} is large, the system tends to be stable, while when λ is large the system tends towards instability. This is intuitive as the interface is expected to be flat when surface diffusion dominates and vice versa.

It is of interest to note that the stability of the interface is not influenced by the surface energy in equation (3.6). This is strictly true only when the perturbation is small. When the perturbation is large, the first-order Taylor expansion that was made while deriving equation (3.4) is not valid. Further, the curvature-dependent reaction rate is also valid only within the linear reaction rate regime [25]. In addition, the effect of the frequency of perturbation also plays a role in determining the stability of the interface. Specifically, when the frequency of the perturbation is smaller than

Finally, in the above expression, only *ρ*_{Cu} remains to be computed to fully determine the stability of the system. In general, *ρ*_{Cu} will not be constant along the interface and will be curvature dependent. The exact nature of the relation will depend on whether the system is reaction limited or diffusion limited. Next, the stability parameter value, and its implications are discussed in more detail.

### (a) Scalloped versus flat interface

Of importance is the variation of the stability parameter *ζ* with the frequency *ω*. The variation of *ω* is seen in figure 3. There exists a value of *ω* at which a given material system is least stable. However, increasing *ω* beyond this peak value stabilizes the system.

The initial perturbation of the Cu_{6}Sn_{5} interface can be represented as a Fourier series with *ω* varying from 0 to *ω*’s will decay out, while lower frequency perturbations will grow with time. In particular, the perturbation corresponding to a frequency *ω**.

As mentioned earlier, the Cu_{6}Sn_{5} interface is observed to transform into a flat surface once the solder has solidified. This can be explained by the fact that at low temperatures the ratio of surface reaction to the surface diffusivity is small, thereby making the interface stable.

### (b) Scalloped versus needle-shaped interface

When *ζ*>1, depending on the degree of instability, the interface evolves either into a scalloped structure or a needle-shaped structure [20]. The range of perturbation frequencies that cause instability lies between *D*_{Γ} increases, the interface becomes increasingly unstable, with a needle-shape developing when the surface diffusivity becomes very small in comparison with the surface reaction rate, which is to be expected at high temperatures. It is also known that an interface evolves into a needle-shaped one when the surface diffusivity is zero [11,39]. Figure 4 illustrates the evolution of a sinusoidal perturbation under a law of the form
*et al.* [40].

### (c) Copper concentration at the interface

To fully specify *ζ* in equation (3.6), the concentration of Cu at the interface is necessary. The concentration of Cu at the interface depends on the applied boundary condition. As mentioned in §2c, the boundary condition depends on whether the system is reaction limited or diffusion limited. When the system is diffusion limited, the concentration of Cu is obtained by solving the Neumann boundary condition

On the other hand, when the system is interface reaction limited, the concentration of Cu at the interface is determined by the Gibbs–Thomson condition (equation (2.23)). This implies that the concentration of Cu at the negative curvature regions is greater than in the positive curvature regions tending to decrease the stability of the interface as indicated by the relation of concentration to stability parameter in equation (3.6). Thus, when the system is reaction limited the interface tends towards instability even with the same values of λ and *D*_{Γ}. Consequently, the interface frequently destabilizes into a dendritic structure.

The concentration of Cu at the interface is also strongly determined by the direction of electrical current. The bulk diffusion of Cu is affected by the direction and magnitude of the electric field. When the direction of current causes to increase the flux of Cu towards the Cu_{6}Sn_{5} interface, the Cu available for reaction increases and, therefore, the Cu_{6}Sn_{5} grows at a relatively faster rate and vice versa. This has been experimentally observed and numerically simulated in prior studies [16,41]. This phenomenon is also demonstrated later in this paper in §5a.

## 4. Numerical implementation

A stand-alone parallel, adaptive finite-element implementation built on the ‘libMesh’ library [42] was used to solve for the field quantities and to evolve the interface in this study.

### (a) Formulation

The primary challenges to developing physical models of moving interfaces are the representation of the interface and the application of boundary conditions at the moving interface. There are two thermodynamical treatments of the evolving interface that are commonly used in the literature. The first is the Korteweg idea of the interface where the interface is treated as a region of sharp gradients in the field parameters. Alternatively, the interface can be considered to be an infinitesimally thin region associated with the *interfacial* quantities. The equations of motion can then be developed in terms of this description of the interface as has been done in §2. In this work, the latter of the above two approaches in which phase field method is used as a purely numerical tool is used. The values of the phase field is used to indicate the Cu_{6}Sn_{5} phase and the Sn phase with a transition region with thickness of the order *O*(*ϵ*). The evolution of the interface is controlled by a degenerate Cahn–Hilliard equation:
*f*(*ϕ*) used here is the biquadratic function (figure 5*a*). The asymptotic analysis also shows that the solutions of equation (4.2) lead to a hyperbolic tangent variation of the phase-field function over the interface ranging from −1 to +1, with the extreme values representing the two phases, and values in between representing the interfacial region. This function can now be used to construct smoothed approximations of the Dirac-*δ* function and the Heaviside step function shown in figure 5*b*,*c*:

The above functions can be used to modify the material properties for the solution of the driving fields. The Dirac-*δ* function can be used to apply boundary conditions at the moving interface. The differential equations of the model modified with the Heaviside step function and the Dirac-*δ* function can be seen in the second column of table 1.

### (b) Solution procedure

As a general procedure the electrical, thermal and stress fields are solved for first, followed by the solution of the Cahn–Hilliard equation to update the interface. This is then followed by the solution of the diffusion equation. The solution for the electrical and stress problems are computed by the finite-element method. Specifically, the stress problem is solved with a standard displacement based finite-element method. The form of the Heaviside step function equation (4.4) can lead to a singular system matrix. To ameliorate this, a small number is added to the heaviside function during assembly. The model for the diffusion of Cu atoms in the bulk require the computation of the gradient of pressure. Quantities related to stress are computed in a post-processing phase after the displacements are determined. These values are found at the Gauss points, but, as the stress quantities depend on the gradients of the primary displacement variables, the level of continuity and differentiability is reduced. For example, with linear *C*^{0} elements, the displacements are continuous at the nodes, while the gradients of the displacement are discontinuous across element boundaries. This makes the computing of the gradients of stress quantities difficult. The stress quantities are, therefore, computed at the Gauss points and then extrapolated to the nodes by means of a global least-squares solution.

The modified Cahn–Hilliard equation used in this work is a non-convex nonlinear differential equation. A fully implicit solution of the nonlinear equation is computationally expensive as the mesh needs to be very densely refined in the interfacial region. In this work, a splitting method based on [43,44] is used. The phase field function is split into the sum of a convex and a concave part (figure 6). The convex part is treated implicitly, while the concave part is treated explicitly based on the previous time-step. This allows the problem to be time-stepped in a very efficient manner, as by choosing the convex part to be a quadratic function, each time-step is a linear solution step. The other challenging equation to solve is the diffusion equation. This is an advection–diffusion equation and frequently needs stabilized solutions when the advective term dominates. However, the diffusive term is significant in all the cases that are considered here. The equation is, therefore, solved with standard finite elements, and time-stepped using a backward Euler scheme.

### (c) Adaptivity

In order to reduce the computational cost of the phase field method, the problem is solved over an *h*-adaptive mesh. The adaptivity features in this code are provided by ‘libMesh’. In this study, the mesh is adapted solely on the basis of the phase field value. Elements are marked for adaptivity based on the average value of the phase field parameter over the element. This is computed as
*ϕ*^{av}_{el} is between −0.97 and 0.97, the element is marked for refinement, otherwise, the element is marked to be coarsened. This method has been preferred to the alternative of using an error estimator as numerical examples have shown the solution to be more reliable with time-stepping.

## 5. Simulations of interface growth and stability

The material properties used in the simulations are tabulated in table 2. In the pictures containing the simulation results, black denotes the Sn phase and lighter grey represents the Cu_{6}Sn_{5}. It should be noted that the Sn could exist in either the solid form or the liquid form depending on the temperature. Temperature is not explicitly included in the simulations, but is present implicitly as its effect is to modify the material properties such as the diffusivities and the reaction rates.

The boundary conditions that were applied in the simulations are shown in figure 7. In all of the simulations, an unlimited supply of copper is assumed at the bottom boundary as would be the case in microelectronic assemblies. A wetting angle of 90° is chosen for the interface–sidewall contact; this improved the stability of the code as sharp corners at the edges caused the concentration solution to diverge. Electrical boundary conditions were applied by imposing an electrical potential at the top and bottom surfaces as shown in figure 7.

### (a) Diffusion-limited versus reaction-limited growth

As mentioned earlier, experimental observations of IMC growth appear to indicate a diffusion-limited system inferred from a sub-linear rate of growth. In a reaction-limited system on other hand, the concentration of copper at the interface will remain a constant and the growth rate will be linear in time. In this section, each of the earlier discussed boundary conditions corresponding to a diffusion- or reaction-limited system are applied in independent simulations, and the expected parabolic or linear growth is demonstrated to result from the simulations.

In figure 8, a plot of the normalized thickness of Cu_{6}Sn_{5} layer as a function of time is shown. It is seen that the growth rate is parabolic as is expected for a system that is limited by bulk diffusion and is linear in the case where the Gibbs–Thomson condition is imposed. Also seen in figure 8 is the growth dependance on the direction of electric current, which is described next.

The rate of Cu_{6}Sn_{5} growth is strongly affected by the addition of electric current due to electromigration. A strong dependence of the growth rate on the direction of current has been reported in the literature [45]. As mentioned previously, the Cu_{6}Sn_{5} interface is known to evolve into a flat surface once the solder solidifies. This was explained earlier as being due to the domination of surface diffusion over surface reaction rate upon solidification of the solder. In the simulations, to study the rate of Cu_{6}Sn_{5} growth with respect to the direction and magnitude of current, the ratio of λ to *D*_{Γ} was set to 10. The value of *Z** in the Cu_{6}Sn_{5} was set to 1.0. The electrical potential of the top surface was set to 0.0 and the potential of the bottom surface was varied as indicated in figure 9. The electric current drives copper diffusion, leading to greater or lesser supply of copper (depending on current direction) at the interface when compared with the case without electric current. When the current direction causes to increase the supply of Cu, the Cu_{6}Sn_{5} growth is enhanced and vice versa.

### (b) Unstable interface and copper concentration variation on interface

The simulations were initialized with a sinusoidal perturbation to the Cu_{6}Sn_{5}–solder interface and then allowed to evolve. The reaction rate, λ was set to 6 and the surface diffusivity, *D*_{Γ} was set to 0.01, so that the interface did not stabilize to a flat surface. As seen in figure 10, the sinusoidal perturbation leads to formation of the scalloped pattern, similar to what is observed in experimental studies.

A plot of the concentration field in figure 11 clarifies that this curvature dependent growth is not due to the Gibbs–Thomson effect, but, rather, it is due to the reaction kinetics. The concentration of Cu at portions of the interface that are growing most rapidly (i.e. regions of negative curvature) is in fact lower than in regions of positive curvature as predicted in §3c. This is opposite to what one would have observed had the Gibbs–Thomson effect manifested.

### (c) The effect of curvature on stability

The effect of the curvature of the interface can be analysed by observing the interface velocity as a function of the frequency of the sinusoidal perturbation of the interface. It is seen in figure 12 that surface diffusion dominates when the interface curvature is high, while interface reaction dominates when curvature is low.

The initial profile of the physical interface can be modelled as a Fourier series with components of various frequencies. This result here implies that the size of the scallops that form in experiments will be related to the frequency of perturbation for which there is a maximal interface instability as discussed in §3a. It has been observed experimentally that the observed mean scallop width in experiments is ≈4.5 μm [46]. This width is known to be a function of processing temperature, which is to be expected as λ and *D*_{Γ} will vary exponentially with the inverse of temperature.

The dependence on the interface perturbation frequency is more clearly seen in figure 13, where the initial Cu_{6}Sn_{5} interface is modelled as a superposition of two sinusoidal perturbations of different frequencies. The perturbation with the higher frequency decays quickly, while the one with lower frequency grows into scallops. This provides an explanation for why the observed Cu_{6}Sn_{5} scallops are roughly the same size in all experimental studies. The size appears to be the result of late-stage growth dynamics rather than the initial nucleation pattern as has been proposed in the literature [23].

### (d) The effect of reaction rate and surface diffusivity on stability

The effect of the reaction rate λ and the surface diffusivity *D*_{Γ} are complementary. As was shown in equation (3.6), it is the ratio of these two parameters that determines the steady-state morphology of the interface. In figure 14, a simulation result with surface diffusivity of 0.1 is contrasted against a surface diffusivity value of 0.01 that was used in the simulation shown in figure 10. All simulations shown in figure 14 possessed the same initial frequency as that shown in figure 10. It is expected that once the solder solidifies, the ratio of the reactivity to the surface diffusivity causes a flat interface to become the stable solution.

If, instead of increasing the surface diffusivity value relative to the nominal simulation of figure 10, it were decreased to a value of 0.005, the effect is to elongate the scallops as shown in figure 14. In general, the code becomes unstable when the surface diffusivity is reduced below this value, challenging the simulation of needle formation.

### (e) The role of initial configuration on interface stability

In order to assess if the observed scalloping phenomena is dependent on the initial configuration of the interface, the following simulations were run. The interface was set to be completely flat initially and the contact angle at the sidewalls was varied. It is seen in figure 15 that in the case where the contact angle at the sidewalls was 90°, the interface remains flat even though the material parameters were such as to cause instability. This is because there are no destabilizing perturbations on the interface.

By contrast, when the contact angle at the sidewalls was changed to 80°, the interface spontaneously destabilized into the scalloped morphology. In practice, the interface is unlikely to be devoid of any perturbations. Deviations from perfect flatness can arise due to the presence of grains, nucleation sites, sidewall contact angles or surface roughness. When the material parameters are conducive for scallop formation, these perturbations will grow with time.

## 6. Summary

A framework for modelling diffusion–reaction systems, which is at the core of numerous phase evolution problems, was developed and described in this study. The example problem of intermetallic growth in a Cu–Sn system was analysed in detail. The growth of Cu_{6}Sn_{5} is driven by interdiffusion of Cu through the Cu_{6}Sn_{5}, the surface reaction at the interface to form Cu_{6}Sn_{5} and surface diffusion of the Cu_{6}Sn_{5}. Experimentally, the Cu_{6}Sn_{5} interface is known to take on a variety of morphologies. A mechanism based on an interface curvature-dependent reaction rate is demonstrated in this paper as an explanation for the observed interfacial morphologies, in opposition to a mechanism based on the Gibbs–Thomson effect. The latter effect is shown to lead to linear growth rates of the intermetallic compound, which is at variance with experimental observations. A stability analysis was performed to identify the critical parameter grouping that governs the equilibrium shape of the interface. Finally, simulations were carried out using a diffuse-interface computational tool demonstrating the effect of various material parameters. The growth rate of the Cu_{6}Sn_{5} was found to be parabolic as has been observed experimentally. The direction and magnitude of current was found to alter the growth of the Cu_{6}Sn_{5} layer significantly.

## Authors' contributions

A.U., S.S. and G.S. designed research. A.U. performed research. S.S. wrote the phase field program. A.U. and G.S. wrote the paper.

## Competing interests

The authors declare no conflict of interests.

## Funding

The authors acknowledge financial support from the Semiconductor Research Corporation under task 1292.090.

- Received February 23, 2016.
- Accepted May 23, 2016.

- © 2016 The Author(s)

Published by the Royal Society. All rights reserved.