Numerical Simulations of Electromigration and Stressmigration Driven Void Evolution in Solder Interconnects PUBLIC ACCESS

[+] Author and Article Information
Subramanya Sadasiva

 Purdue University, West Lafayette, IN 47907subramanya@purdue.edu

Ganesh Subbarayan1

 Purdue University, West Lafayette, IN 47907ganeshs@purdue.edu

Lei Jiang

Intel Corporation,Portland, OR 97124lei.jiang@intel.com

Daniel Pantuso

Intel Corporation,Portland, OR 97124daniel.pantuso@intel.com


Corresponding author.

J. Electron. Packag 134(2), 020907 (Jun 11, 2012) (9 pages) doi:10.1115/1.4006707 History: Received September 20, 2011; Revised February 13, 2012; Published June 11, 2012; Online June 11, 2012

Understanding the effect of high current density on void formation and growth and relating the size of the void to the resulting electrical/mechanical failure is a critical need at the present time to ensure reliable functioning of flip-chip packages. In general, toward this end, the modeling and simulation of geometrical evolution of current induced voids have been relatively few. Simulations considering the coupled effects of mass transport through mechanisms of surface and bulk diffusion under the influence of electrical, thermal, and stress fields in solder joints leading to eventual electromigration failure do not appear to be common. In this study, we develop a phase field model for the evolution of voids under electrical, thermal, and stress fields in a flip-chip solder interconnect. We derive the equations of motion for the void accounting for energetic contributions from the active factors of surface energy, stress, and electric potential, considering both surface diffusion and transfer of the material through the bulk of the material. We describe the implementation of this model using a finite element code written in the PYTHON language, coupled with a commercial finite element solver from which we obtain the electrical, thermal, and stress fields driving the void motion. We demonstrate the implemented methodology through simulations of void evolution in flip-chip solder joints under the effects of mechanical/electrical fields and surface/bulk diffusion.

The decrease in the size of the flip-chip solder joints has resulted in a corresponding increase in the current densities prevalent in the joints. This has led to an increase in the importance of failure due to electromigration in flip-chip solder joints [1]. Electromigration is caused by the momentum exchange between electrons and ions that forces material diffusion in the direction of electron flow. This leads to the formation of voids near the cathode that might lead to an open failure. There is also the likelihood of the formation of hillocks at the anode that might cause a shorting failure. The formation and motion of the voids are also driven by the existence of stress and temperature gradients. A unique problem that occurs in solder joints is the failure due to current crowding. The current density in the signal trace leading up to the joint is significantly higher than the current density in the solder joint. At the interface between the two, there is a region of the solder joint, where the current density is an order of magnitude higher than the rest of the joint (≈105 A cm−2 at the entrance versus ≈104 A cm−2 near the middle [1]). This causes the void formation to be very rapid near the current crowding region (see Fig. 1).

The computational simulation of the formation and propagation of voids in solder is a challenging problem. There are a multitude of interacting physical effects such as the electrical, thermal, and stress fields that determine the void shape and its rate of evolution. A significant attempt at modeling electromigration in solder joints has been through the use of damage mechanics to simulate the loss in load-bearing capacity of the solder joints [2-3]. These models enable one to find locations of void formation and in modeling the motion of very small, diffuse, voids, whose boundaries are not explicitly captured. These methods do not explicitly model the geometry of the voids, but infer the effects of the formed voids through the damage they cause to the load-bearing capacity. However, geometrically speaking, the removal of material from a region results in the formation of new surfaces. Further, energywise, every new surface formation has an associated surface energy cost [4] that needs to be accounted for in any complete simulation of the motion of voids.

An alternative approach to modeling electromigration has been to capture the divergence of atomic flux directly within a finite element code [5-6]. In these approaches, the modeling of motion and growth of larger voids appear to rely mostly on element deletion schemes, which approximate the boundary of the voids accurately in the limit of mesh refinement [6]. Furthermore, these methods also ignore the surface energy associated with void boundaries or surface diffusion along the void boundary, and rely on the underlying physics of electromigration being captured entirely in the bulk through the divergence of flux.

As the voids evolve, the underlying geometrical domain on which the boundary value problems are solved changes. One way to handle this is by explicitly tracking the motion of the void surface. There are very few attempts in the current literature aimed at explicitly modeling the evolution of electromigration-induced voids [7-8]. The main challenge with explicit geometry tracking methods is handling the constantly changing problem domain, which necessitates remeshing with a mesh-based numerical scheme such as the finite element method. Also, in an explicit scheme, the processes that cause topological changes such as the splitting of the voids require special treatment and are in general very difficult to model numerically.

An alternative to this method is to use a phase field model for the void boundary. In phase field models, the existence of a boundary is indicated by a change in an order parameter. Quantities that are defined along the surface of the boundary are smeared over a narrow region surrounding the boundary. The main advantage of using a phase field model is the ability to retain a fixed finite element mesh for the simulations even as the geometrical shapes of the voids evolve. Also, changes in the void connectivity such as the merging and splitting of voids are handled automatically and efficiently. This method has been used earlier for the modeling of void evolution in Al–Cu interconnects through surface diffusion [9].

In this study, we develop a phase field model for simulating the motion of pre-existing voids in flip-chip solder joints due to diffusive mass transport along the surface of the voids (and therefore preserve the volume of the void) as well as diffusive mass transport between the voids (and therefore do not preserve the volume of the individual voids). We implement a finite element solution methodology, including the coupling to a commercial finite element code that is used to estimate the electrical, thermal, and stress fields driving the void evolution.

The structure of the rest of the paper is as follows. We first describe the basis for our governing equation derived from a statement of the minimization of a free energy rate. We follow that with a description of the phase field model that is used for the solution of the void propagation problem. Finally, we illustrate the utility of this method by studying the motion of the void in a flip-chip solder joint, under general loading conditions.

Minimization of Free Energy Rate

Considering the system shown in Fig. 2, the evolution of a near-equilibrium system can be described by the minimization of the dissipation or of a free energy rate. The free energy rate Ψ· in the void evolution problem has three sources, namely, surface effects, electrical field, and due to mechanical load Display Formula

In the above rate, the free energy of the void surface, Γ can be written as Display Formula
Thus, the rate of change of the surface free energy is given by [10] Display Formula
The surface free energy measures the energetic cost of adding or removing material from the surface of the void. The role played by surface energy in solids is slightly different from that in liquids. In liquids, the surface energy affects both the stress equilibrium of the material and the cost of adding and removing the material. However, in solids the surface energy does not play a significant role in the determination of the deformation equilibrium. The inclusion of a surface term in the expression for the energy makes it possible to model scale effects.

Assuming that a charge density ρc exists everywhere in the domain, the net electrical potential energy is given by ρc φe , where φe is the electric potential. The total electrical energy of the system at steady state is given by Display Formula

The rate of change of this quantity, assuming equilibrium is reached for electrical potential much faster than the rest of the system, is obtained by taking the material time derivative [11] of Eq. 4Display Formula
where [[]] is the jump in the quantities across the interface defined as [[a]] = a+  + a− , referring to the quantity on either side of the interface. Assuming that the material velocity is significantly slower than the time taken to reach steady state, we can assume that the system is always at steady electrical state. Hence, the spontaneous change term, (ρcφe)t0. Also, if the external boundaries of the system are fixed, Γext(ρcφe)vndΩ=0. Further, for a material–void system, one side of the jump term can be ignored. Hence, the free energy rate due to electrical work can be approximated in electromigration-induced void evolution problems as Display Formula
Similarly, the contribution to free energy from mechanical loads is Display Formula
where φm is the strain energy density. To evaluate the mechanical energy contributions, from the balance of mechanical energy on the system, we get Display Formula
where D is the rate of deformation tensor. Assuming linear elasticity (σ = E: ɛ), and therefore the strain energy density φm=(12)ɛ:E:ɛ, we have, after similar considerations for equilibrium of the stresses Display Formula
Arguing, as before, that the system reaches equilibrium instantly relative to the times involved in mass transport, we have φmt0. Assuming fixed external boundaries, and no loads in the void, we get Display Formula
Combining Eqs. 3,6,10, we get Display Formula
Following the arguments from Ref. [12], and introducing surface and bulk mobilities Ms and Mb , we get the following equations for the velocity of the front.

For surface diffusion, assuming a local conservation of mass Display Formula

Display Formula

For the bulk diffusion, we assume that the mass is transferred only from one void surface to another. This leads to an inherently nonlocal equation for the bulk diffusion equations. The equation for the velocity of the void surface is given as Display Formula

where λ is a Lagrange multiplier that enforces conservation of mass throughout the system. This constraint can be written as Display Formula
Finally, the motion under coupled surface and bulk diffusion effects is thus described by Display Formula
This is similar to the Eq. (1.14) derived in Ref. [13]. In Ref. [13], the authors derive the equations using the principles of rational thermodynamics. While the first term in the two equations is the same, the interpretation of the second term differs. While the authors in Ref. [13] use the second and third terms as modeling evaporation and condensation effects from and to the surface, here, on the other hand, we use the terms to mean any rapid transport of material through the bulk of the material.

Diffuse Interface Equations

The equations derived above give us the evolution equations, based on a sharp interface description of the void interface. These have the problem of being difficult to handle when the shape of voids evolve or when there are changes in the topology of the voids since these would require remeshing of the domain. A related challenge is on utilizing the calculated void motion defined on the sharp boundary in a commercial finite element analysis software used to analyze the electrical/mechanical behavior of the system. A diffuse interface approximation for the void interface functions by smoothing the surface over a notional thickness of order ε and indicating the surface through the level sets of a phase field variable u that varies from − 1 inside the void to + 1 outside the void (see Fig. 3). The surface energy (Eq. 2) is written as Display Formula

where cF is a constant scaling parameter defined as in Ref. [14].

The gradient energy was first used in Ref. [15] to study spinodal decomposition in alloys. It was later extended to study other problems with moving boundaries in a diffuse sense, such as solidification [16]. The gradient energy allows the treatment of surface quantities in a diffuse volumetric sense. Locations with nonzero gradients indicate the presence of the interface. The normal to the interface can be inferred by (u/u|u||u|). The term f(u) in Eq. 17 is a bulk free energy term, with two minima that correspond to the two phases, in this case, the material and the void. The presence of (f(u)/f(u)ɛ) in the energy functional penalizes solutions that are not ±1.

A good choice for the potential function is the biquadratic function f(u)=(1-u2)22 shown in Fig. 4, cF=43, and the notional thickness of the interface is 2ɛ. Using this description, the equations for the motion of the void can be written in terms of partial differential equations in u. We have, for the surface diffusion problem Display Formula

Display Formula
where Ms (u) is restricted to the interfacial region where −1 ≤ u ≤ 1. The equation for diffusion through the bulk is given by Display Formula
Display Formula
The combined motion of the void surface due to both the surface and the bulk effects is obtained by adding Eqs. 18,20. Equation 18 is a specialized version of the Cahn–Hilliard equation [15], while Eq. 20 is a version of the Allen–Cahn equation with an additional constraint on the phase field variable [17].

The equations for the void motion are solved in a two-step fashion solving first for the surface diffusion and then solving for the diffusion through the bulk. The partial differential equation for the surface diffusion is 4th order in space. Usually, these equations require C1 continuity across element boundaries when using finite element analysis. However, by splitting the governing equations into two 2nd order equations, we can solve these equations using C0 elements following the approach outlined in Ref. [18]. As there are no external fluxes, we can assume Neumann boundary conditions everywhere on the boundary.

In order to make the solution procedure robust, we scaled the length by the dimension of the solder joint, r, and the time by an arbitrary scale, r4Msγ. Further, we scaled the electric potential by the applied potential, and the strain energy with the maximum strain energy observed in the structure. On scaling, Eqs. 18,19,20,21 become Display Formula

Display Formula
Display Formula
Display Formula
where ke=ρcφ0eγ, km=φ0mγ, and kd=Mbr2Ms. kd controls the relative effect of surface diffusion with respect to bulk diffusion. The other parameters control the relative effect of the quantities with respect to the motion due to the reduction of surface energy. Specifically, ke controls the relative effect between the electric potential and the surface energy, km controls the relative effect between the strain energy density and the surface energy, and kd controls the relative effect of bulk and surface diffusion velocities.

Discretized Forms of Governing Equations

Multiplying the equations for surface diffusion (Eq. 18) by a test function w and integrating by parts, we get Display Formula

Display Formula

Discretizing the above equations using C0 continuous linear triangular elements (for convenience, since it assures constant derivatives), and assuming that dudt=ui+1-uiΔt, we get Display Formula

Here, A is the lumped mass matrix for the system, BΩNiNjdΩ and CΩNiM(ui)NjdΩ. Similarly, for the bulk diffusion equations, by multiplying Eq. 20 by the same test function w, we get Display Formula
Display Formula
Again, discretizing using C0 linear triangular elements and using an Euler backward scheme, we get the following matrix equations: Display Formula
Here, a are the diagonal terms of the lumped mass matrix A.

The coupled effect of bulk and surface diffusions are accounted for by solving the surface diffusion equation first and then solving the bulk diffusion equation. This is equivalent to solving Eq. 28 first and then solving Eq. 31 using the value of ui from the solution to Eq. 28.

The solutions for the strain energy, electric potential, and temperature values are obtained by solving the corresponding boundary value problems in a commercial finite element code, using C0 finite elements on the same mesh used for the solution of the phase field problem. The void is indicated by modifying the material properties for Young’s modulus E and electrical conductivity η as Display Formula

E(u)=E0(1+u)2,   η(u)=η0(1+u)2

Challenges to Numerical Solution

As the Cahn–Hilliard equation is a fourth order nonlinear equation, time integration is difficult. Simple explicit solutions such as the forward Euler method are precluded as the time step required for numerical stability is ≈ch4 , where h is the size of the element and c is an arbitrary constant. Hence, a semi-implicit method is used for time integration. The mobility and the forces are treated explicitly and are computed and held constant from the end of the time step, while the phase field variable itself is treated implicitly and intreated for. A further restriction on the time-step size is imposed due to the nonconvex nature of f(u). This causes the solution to converge uniformly to either +1 or −1 for large values of time steps. In our studies, we found that the numerical time steps necessary for the convergence of the non-linear solver were on the order of 10−4 to 10−5 . As this problem is solely due to the numerical difficulties presented by the Cahn–Hilliard equation, the stress and electrical field solution was not recomputed at every time step, but only computed once every few steps of the solution of the diffusion equations. This was also beneficial in reducing the number of times that the commercial solver needed to be called for computing the electrical and stress fields.


The boundary value problem to obtain the electric field values was solved using a commercial finite element code, using thermoelectric elements. Similarly, the stress problem was solved using plane stress elements. The finite element analyses relied on linear triangular elements. This was linked to the solver for the phase field equations that was written in PYTHON language using SCIPY and SCIPY [19] for the numerical data structures and solvers, respectively. The linear systems were solved using the UMFPACK sparse direct solvers. The code developed, pyPhase, consists of approximately 1200 lines of PYTHON code. The PYTHON code includes capabilities to read in information from commercial finite element tools and creating and running simulations. It has the capability of handling arbitrary geometries of solder joints and also assemblies that include multiple joints. The code is independent of the choice of commercial finite element software. The structure of the code is described in Fig. 5.

To validate the model and its implementation, it was tested on rectangular line geometries. In a line that is under the influence of surface diffusion, circular voids would collapse to slits under axial loading, while under electrical loading, the void would move toward the cathode. This is shown in Figs.  67.

We also tested the solution of the coupled solution including both bulk and surface diffusion. In the absence of any external loading such as electric potential gradients or stress, the system is expected to evolve purely by reduction of the surface free energy. While under surface diffusion, multiple noncircular voids are expected to relax to separate circular voids, with bulk diffusion allowed, smaller voids are expected to lose volume and eventually disappear. This can be tested by varying the relative rate of surface and bulk diffusion.

In Fig. 8, the relative rate of surface and bulk diffusion kd is varied from 0.1 to 10 from left to right, as can be seen; at high values of kd , the bulk diffusion effect dominates and the smaller void disappears. On the other hand, when the value of kd is lowered, the area of individual voids is conserved and there is only a relaxation due to the surface diffusion equation.

The 2D model (Fig. 9) was constructed in a commercial finite element tool, according to dimensions from Refs. [(20-21) ]. The material properties used were those of tin (Table 1). The chip side and board side pads were assumed to be made of copper. We assumed that the chip side was at ground, while the board side was maintained at 5 V. The solder joint itself was assumed to be under isothermal conditions. A spatially varying shear displacement was applied on the chip side to simulate the effect of the thermal expansion mismatch. This displacement varied from a value of 0% to 10% of the pad diameter across the diameter of the joint. The interface between the solder joint and the pads was treated as diffusion barrier. This prevents the void from propagating into the pads.

Given below are the results of the simulations on the flip-chip solder joint. The problem was scaled and nondimensionalized for numerical stability as described in Eqs. 22,23,24,25. Specifically, we chose r, the characteristic length, to be the midsection diameter of the solder joint. As described earlier, we also used a normalized timescale τ=r4Msγ.

Surface Diffusion

We first demonstrate our model on examples that are limited solely to surface diffusion. In this case, Eq. 20 is deactivated and not solved. The first example illustrates the effect of strain energy alone on the void motion (Fig. 1). In the next example, we simulate the motion of a void surface due to the electrical effects alone (Fig. 1). We see that in these cases the strain energy density and electrical effects cause the void to evolve in opposite directions. The evolution under an electrical field leads quite quickly to a pancake shaped void that has the potential to cause electrical failure. This is due to the fact that the electric current pushes the particles away from the current crowding region. Hence, it is expected that the voids in the solder migrate to this location. For the case with strain energy driven diffusion, the reason for the direction of the void surface motion is that even though the mechanical loading is a horizontal shear force. The direction of the strain energy gradient is tangential to the surface of the void. This causes the void to tend to collapse in that direction (see Fig. 1).

Combined Bulk and Surface Diffusion

In real applications, both surface and bulk diffusions are significant. This leads to a void growth analogous to Ostwald ripening effect observed in heterogeneous microstructures where one void grows at the expense of another. Shown in Fig. 1 is a case where there is a significant diffusion through the bulk. It is observed that the system evolves by transferring the mass to the void closest to the current crowding corner. The farther void reduces in size, until it is completely dissolved. The reason for this location of the final void configuration is the same as the explanation for the position of the void for the case with only surface diffusion.

Simulating Electromigration in Assemblies of Solder Joints

The electromigration in real solder joints is subject to effects that are visible only in assemblies of solder joints. One of these is the effect of the direction of the current in the current crowding region on the propagation of voids. In this example, an assembly of two solder joints is considered. The dimensions and loading of this assembly are described in Fig. 1. The entire joint is loaded by applying a potential difference of 10 V across the two solder joints. Additionally, a shear force is simulated on the two solder joints by applying a deformation of 10% of the midsolder section diameter on the joints (13.2 μm) in each of the joints. By loading the assembly in this fashion, one is able to consider the effect of the direction of the current in the current crowding region on the motion of the voids since, for instance, the top left corner of the first joint and the top right corner of the second joint experience opposite current directions, while maintaining the mechanical load direction to be the same.

In the first case (Fig. 1), the effects of the strain energy driven diffusion are turned off by setting km  = 0. The direction of the current in the left solder joint causes the void to move into the current crowding region and form a pancake like void in that region, while in the right joint, the direction of the current is reversed and the void is pushed into the solder joint. Figure 1 illustrates the effect of mechanical load alone on the evolution of the void, and provides the baseline for the study of the effect of strain energy distribution on the electromigration of the void. As before, to understand the reason for the evolution of the void into the bulk of the solder joint, it is illustrative to study the strain energy distribution in the solder joint (Fig. 1). In Fig. 1, the combined effects of both electromigration and stress-induced void evolution are illustrated. As can be seen, the strain energy gradient causes the void to evolve in a direction counter to the direction of the electromigration. In this case, it slows down the evolution of the void toward the current crowding region. This suggests that the severity of failure due to the electromigration of voids might depend significantly on the position of the solder joint in the package, depending on the type and severity of the mechanical loading seen by each joint.

We demonstrated, in this work, a phase field method to simulate the motion of pre-existing voids in the flip-chip solder joints. We detailed the sharp interface equations for the motion of a void in a material evolving under the influence of surface energy, strain energy, and electric potential. This is related to diffuse interface equations that reduce in the limit of ɛ to the sharp interface equations, and tracks the motion of the interface. The equations presented here account for the surface diffusion and also the bulk diffusion between multiple voids. The numerical solution to these equations is implemented into a numerical framework that uses a commercial finite element code to solve for the stress and electrical potential fields. We use this system to study the dynamics of void motion in a two-dimensional model of a flip-chip solder joint assembly. We showed the effects of introducing a void in the current crowding region of the solder interconnect. It was found that the electromigration void evolution depends on the directionality of current, with a pancake void forming in the cathodic region and the void migrating into the interior of the joint away from the interface near the anodic region. Further, the stress-induced void migration direction is opposite to that of current-induced void migration, and as a result, failure is likely to depend on the location of the solder joint in the array, the presence of significant shear loading, and the direction of current. Electromigration in solders is a highly complicated process and while the model shown here gives us an indication to the likely direction of the motion of the void, the role of plasticity, and the initiation of voids needs to be accounted for in the models in a thermodynamically consistent way. This is currently under investigation.

The authors are grateful for support from Intel Corporation. The authors would also like to thank Dr. Yu Xiao for discussions.

vn =

normal velocity at a point on the boundary of the void

u, f(u) =

phase field order parameter, phase field potential function

cF =

velocity scaling parameter

Ψ =

free energy

γ =

surface energy density

E =

Young’s modulus

η =

electrical conductivity

φe , φ0e =

electric potential in domain, applied electric potential

φm , φ0m =

strain energy density in domain, maximum strain energy density

κ =

mean curvature

ɛ =

notional interface thickness

ρc =

charge density

Ms =

surface diffusion mobility

Mb =

bulk diffusion mobility

kd =

scaling parameter accounting for the relative magnitudes of bulk and surface mobility

λ =

Lagrange multiplier for global conservation of mass

τ =

scaled time

Copyright © 2012 by American Society of Mechanical Engineers
View article in PDF format.



Grahic Jump Location
Figure 1

Formation of pancake voids in the current crowding region [7]

Grahic Jump Location
Figure 2

Schematic showing the geometry and boundary conditions needed to solve the electromigration problem

Grahic Jump Location
Figure 3

Diffuse description of the interface

Grahic Jump Location
Figure 4

The biquadratic potential used in the simulations

Grahic Jump Location
Figure 5

Schematic structure of the code

Grahic Jump Location
Figure 6

Migration of circular void under electric potential gradient

Grahic Jump Location
Figure 7

Collapse of circular void into a slit under axial tension

Grahic Jump Location
Figure 8

As the value of kd is increased, the effect of bulk diffusion is increased and the bulk diffusion starts to dominate

Grahic Jump Location
Figure 9

Schematic diagram of the solder joint over which the diffusion problem is solved

Grahic Jump Location
Figure 10

Example showing surface diffusion due to strain energy alone. Here, kd  = 0, ke  = 0, and km  = 40.

Grahic Jump Location
Figure 11

Example showing surface diffusion due to electric potential gradients alone. Here, kd  = 0, ke  = 40, and km  = 0, as time progresses the void rearranges into a pancake formation, similar to that shown in Fig. 1.

Grahic Jump Location
Figure 12

Distribution of strain energy near the surface of the void

Grahic Jump Location
Figure 13

Example showing significant diffusion through the bulk of the material. Here, kd  = 5, ke  = 20, and km  = 20.

Grahic Jump Location
Figure 14

Schematic showing the configuration and loading of an assembly with two solder joints. The shear load is exactly symmetric.

Grahic Jump Location
Figure 15

Evolution of voids in the left and right solder joints showing effect of the direction of the current on the motion of the void. Here, kd  = 0, ke  = 100, and km  = 0.

Grahic Jump Location
Figure 16

Evolution of voids in the left and right solder joints showing effect of the strain energy gradient on the development of the void. Here, kd  = 0, ke  = 0, and km  = 100.

Grahic Jump Location
Figure 17

Strain energy distribution around the void region. The extremely high and low values have been removed in order to show the distribution.

Grahic Jump Location
Figure 18

Evolution of voids in the left and right solder joints showing effect of the strain energy gradient and electrical potential on the development of the void. Here, kd  = 0, ke  = 100, and km  = 100.


Table Grahic Jump Location
Table 1
Material properties used in the simulations


Some tools below are only available to our subscribers or users with an online account.

Related Content

Customize your page view by dragging and repositioning the boxes below.

Related Journal Articles
Related eBook Content
Topic Collections

Sorry! You do not have access to this content. For assistance or to subscribe, please contact us:

  • TELEPHONE: 1-800-843-2763 (Toll-free in the USA)
  • EMAIL: asmedigitalcollection@asme.org
Sign In