The Finite Element Analysis of Weak Spots in Interconnects and Packages

1. Introduction

Concerning the mechanical stress and the electrical-mechanical behavior ULSI multilevel metallization systems are more and more sensitive against influences of geometrical and material changes. The mechanical and electrical reliability of these metallization systems is influenced by this. The reliability of such metallization systems is investigated by thermal and thermal-electrical accelerated stress tests under high temperature load. This leads to degradation due to electro-, thermo- and stress migration. This is one major concern in reliability investigations. Generally measurements are time consuming, expensive and the time-to-market cycle is in the focus of interest too. The prediction of local weak spots in interconnects, vias and solder bumps by finite element simulations are a helpful procedure. Beside this the modern 3-d integration leads to more complex material compositions in the systems concerning the coefficient of thermal expansion (CTE) and other material properties. Higher applied currents on the interconnects and bumps result in Joule heating, high temperature gradients and mechanical stress gradients in the bump and metallization systems. The temperature gradients are much higher compared to systems with wide interconnect lines and bumps with large diameters like in conventional packages. Due to this in ball grid array (BGA) bumps as well as µ-bumps and small through silicon via (TSV) connections current induced migration effects as electromigration (EM) and as result of the high temperature gradients thermomigration (TM) can occur. In interconnects consisting of copper caused by stress gradients due to the different material properties under high temperature load also stress migration (SM) occurs.

In this chapter the degradation phenomena in dual damascene copper metallization structures as well as degradation in bumps, µ-bumps and TSV are presented. The degradation is current, temperature and mechanical stress induced under a high applied current and temperature load. The finite element analyses and the mass flux divergence calculation of these phenomena will show the suitability of the method in comparison with experimental results. For a prediction of the weakest spot the suitability of the finite element mesh as well as the modeling concerning the structure shape has to be investigated. Especially edges in the model can influence the quality of the results. The geometrical data of the different metallization or package structures can be taken from the layout as well as optical, scanning microscope or other analytical techniques. Especially for the prediction of the electromigration induced weakest spot in the system, the location of the maximum current density is a major indicator for the fault location. Due to this the maximum current density at localizations of structure inhomogeneity must be checked up. Out of this the maximum current density is calculated by conformal mapping to predict an optimized modeling concerning the shape of the edges and the use of a radius instead of edges. Geometrical variations like the thickness of the first and second metallization and a comparison of the different migration mechanisms will be presented. Concerning the mechanical stress of the DD-Cu metallization the process induced stress will be considered with different processing temperatures of the copper metallization. The influence of different dielectrics on the mechanical stress is also determined. Compared to DD-metallizations and the traces, the bumps are only exposed by electro- and thermomigration. The thermal-electrical-mechanical behavior of µ-bumps and TSV will be shown for a Wafer-on-Wafer (WoW) structure.

2. Calculation of the migration mechanism

To figure out the possibilities of determining the effects of migration in solder bumps, interconnects, via and conductive pathways shall be shown here. Migration, particularly in copper or aluminum metallization or migration of solder bumps of a flip-chip package, are presented in [Banas, Hou, Liu 2007, Liu2008, Tan, Wang]. Simulation algorithms are based on analogies and allow only partly- or no material or sizing variations.

Basically, migration processes are described in the metallization on the consideration of diffusion processes. These diffusion processes can lead to a change in dimension or change of geometries of the metallization [Wever]. On one side material loss or the hole formation at these locations leads to tensile stress, while it comes in places of material accumulation to a compressive stress. This leads to a reflux, also called the back stress in the material. There is a critical length existing, which compensates the electromigration related flux and the reflux. The migration of atoms leads to a change of the expansion in the metal, which then in turn changes the chemical potential. Since the chemical potential is very often influenced by a stress term, this leads to stress gradients and a reflux. In microelectronic applications the current flux pathways are subject to large thermal strain before any current load applies. This thermal strain, along with the strain caused by the electromigration can lead to a nonlinear reflux. Also, in the interconnects local self-heating take place, providing a contribution to the migration of the thermal gradients due to the local temperature increases. The components of the electrical, thermal and stress migration flux are added by superposition.

The diffusion coefficient with the activation enthalpy of ∆H is defined as follows:

$$D=D_0{e}^{-\Delta H/k_B{T}_{}}$$E1

The change of place takes place at cubic (fcc) metals such as Cu, Au and Ag across the holes. The activation enthalpy is additively composed of the enthalpy of formation and the activation energy for the hole formation. Both are nearly the same size. With equation (1) and D₀ as temperature independent diffusion constant as well as the melting temperature T_m [Heumann, Philibert] applies:

$$5\cdot {10}^{-6}<D_0<5\cdot {10}^{-4}[m²/s],dH\cong 1.5\cdot {10}^{-3}{T}_m[eV]$$E2

Due to this a diffusion coefficient of D(T_m)≈10^-8 cm²/s [Heumann] is resulting.

Parameter	Al	Cu	Ag
Tm	660	1083	960
∆V/Ω	0.71-1.3	0.9	0.7-0.9
E [eV]	1.234	2.15	1.921

Table 1.

Melting Temperature, Active Volume and Activation Energy [Philibert].

In the case of multiple components, for example solder bumps, all main materials must be involved in the material transport. In the example of a material of two components A and B, and the chemical or interdiffusion coefficients, there are the partial diffusion coefficients D_A and D_B [Wilkenson]. Using the Darken equation, the diffusion results as follows:

$$\upsilon =\left(D_A-D_B\right)\frac{\delta N_A}{\delta x}with\tilde{D}=D_B{x}_A+D_A{x}_B$$E3

From the migration velocity υ, the chemical diffusion coefficient $\stackrel{\text{̡}}{D}$ and the mole fraction, x_A and x_B of the considered components the partial diffusion coefficients D_A and D_B can be determined.

The direction of electrical transport in alloys such as leaded solder materials based on Sn, goes, depending on the proportion of Pb to the anode or the cathode. In terms of thermotransport moves the Pb in Sn to the cold part of the sample.

There are different diffusion paths in the material. These are the grain boundaries, surfaces to adjacent materials and the volume or bulk. In the aluminum as a metallization material primarily grain boundary diffusion occurs. In contrary in copper as metallization material predominantly interface or surface migration is occurring.

In the calculation of grain boundaries, the surface diffusion as well as the intermediate phase diffusion must be considered in thin layers. Atomistic transport along a grain boundary or phase boundary has a low activation energy and is therefore by orders of magnitude faster than in the crystal itself. The analytical models for the grain boundary diffusion δ D_GB are only valid for the self diffusion in pure metals [Kaur]. The part of electromigration at grain boundaries depends also on the effective width of the grain δ (d) concerning the mass transportation in relation to the average grain size.

A composition of the interconnect from almost poly crystalline parts and parts where bamboo structured grains occur, leads to the predominant grain diffusion on one hand and on the other to volume diffusion and thus to high gradients in the mass flux. If the mass flux is blocked by a large grain, then the material accumulates relating to direction of the electron or material flux, in front of the blocking grain, while it comes behind the blocking grain to a material loss. For grain boundaries in the area of 280nm the diffusion coefficient of the grain boundary diffusion is in the range of the bulk material D_GB $\not\simeq$ D_Bulk [Kaur].

Only above a defined threshold current density $\stackrel{\text{⃘}}{j}$th, resulting from the Blech effect and denoted as 'short-length effect', it comes to effective place change and thus to the mass transportation or material flux [Blech]. The mass flux is dependent on the atomic particle density N, the Boltzmann constant k_B, the local temperature T, the current density$\stackrel{\text{⃘}}{j}$, the specific resistance ρ, the diffusion process relating activation energy E_A, the diffusion coefficient D₀, and the effective charge eZ.

$$\overrightarrow{J_{EM}}=\frac{N}{k_{B}T}eZ*\left(\overrightarrow{j}-{\overrightarrow{j}}_{th}\right)\rho D_0\exp \left(-\frac{E_A}{k_{B}T}\right)$$E4

The thermotransport denoted as Soret-effect is a phenomenon of overlay of diffusion and heat conduction. They closely resemble the electrotransport with the difference that the considered system is not isothermal [Wever]. There is a flux of soluted atoms and the heat flux. These fluxes are described about the chemical potential and the thermal gradients. Without the occurrence of concentration gradients the equation for thermal flux is:

$$\overrightarrow{J_{TM}}=-\frac{NQ*}{k_B{T}^2}{D}_0\exp \left(-\frac{E_A}{k_{B}T}\right)gradT$$E5

Q* represents here the transported energy at a constant temperature and is commonly referred to as reduced heat of transport or transfer. The ratio Q* / k_BT² is referred as Soret coefficient. The heat of transport is the heat flux per unit of material without temperature gradients. Is the value of Q* > 0 a heat flux is generated to keep the soluted atoms isothermal, which takes place towards the dissolved flux. Is Q* < 0 the flux of dissolved particles and the heat flux are counter set. It follows that in an isothermal system, a density gradient produces a thermal flux and vice versa a temperature gradient leads to a material flux [Shewman]. The heat of transport is approximately equal to the activation energy for the material flux [Jaffe].

The differences of the coefficient of thermal expansion (CTE) between the metallization material and surrounding materials produce, depending on the ambient temperature a mechanical stress. This in turn leads to a material flux within the metallization. Under strain the enthalpy for example of tensile-stress formation in a grain boundary is reduced, which increases the concentration of holes [Heumann].

$$\overrightarrow{J_{SM}}=-\frac{N\Omega }{k_{B}T}{D}_0\exp \left(-\frac{E_A}{k_{B}T}\right)grad{\sigma }_H$$E6

The Ω represents the atomic volume and σ_H the hydrostatic stress. Out of the components of the electrical, thermal, and stress migration, the total mass flux in the metallization structure emerges through superposition from the equations (4, 5 and 6).

The formation of material accumulation or void formation requires a divergence in the flux of the total mass flux. Only divergences lead to a change in the density of the material. The calculation of the mass flux and mass flux divergence in metallization, traces and bumps in the past were done using analogies between the electrical and thermal and the thermal mechanical behavior or a variation of the materials or dimension were not possible [Banas, Hou, Liu 2007, Liu2008, Tan, Wang]. This lead to an overestimation of the temperature gradients and current density like proposed in [Ogurtani]. Without any applied current and temperature gradients the calculation of the stress migration is not possible.

With the new program code the divgrad of T and σ are calculated directly based on the simulation results. The calculation of the mass flux is done for each element under consideration of the neighbor elements. Out of this the stress gradients are calculated. The divgrad terms are calculated on base of the super elements under consideration of an Ansatzfunction. The calculated values for the different migration mechanisms are reloaded into ANSYS for graphic display. Out of this the stress migration can be calculated under SM stress test conditions without any applied current or temperature gradients. The simplified equations neglecting concentration gradients for the different migration mechanism are:

$$\mathit{div}\mathrm{}\stackrel{\text{⃘}}{J_{\mathit{EM}}}=\left(\frac{E_A}{k_B{T}^2}+\frac{{Ϗ}_{ϟ}}{1+{Ϗ}_{ϟ}(T-T_0)}-\frac{1}{T}\right)\circ \stackrel{\text{⃘}}{J_{\mathit{EM}}}\circ \mathit{grad}\mathrm{}T$$E7

$$\mathit{div}\mathrm{}\stackrel{\text{⃘}}{J_{\mathit{TM}}}=\left(\frac{E_A}{k_B{T}^2}-\frac{2}{T}\right)\circ \stackrel{\text{⃘}}{J_{\mathit{TM}}}\circ \mathit{grad}\mathrm{}T-\frac{\mathit{QND}}{k_B{T}^2}\circ \varnothing T$$E8

$$\mathit{div}\mathrm{}\stackrel{\text{⃘}}{J_{\mathit{SM}}}=\left(\frac{E_A}{k_B{T}^2}-\frac{1}{T}\right)\circ \stackrel{\text{⃘}}{J_{\mathit{SM}}}\circ \mathit{grad}\mathrm{}T-\frac{\mathrm{\Omega }\mathit{ND}}{k_{B}T}\circ \varnothing {ϡ}_H$$E9

In equation (7-9) N is the atomic concentration, the Boltzmann constant k_B, the local temperature T, the resistivity ρ, the atomic volume Ω, Q* the heat of transport, the activation energy E_A (taken from grain boundary and interface migration as strongest influence) and σ_H is the hydrostatic stress. The simulation and calculation sequence is shown in figure 1.

3. Modeling and simulation improvements

The numbers of elements, which are determining the mesh and due to this the density of the nodes have a strong influence on the accuracy of the simulation. Out of this the mesh of the investigated metallization or metallic material as main interesting point in the simulation plays a major rule. In the metallization itself the potential as well as the temperature is calculated. The metallization is surrounded by dielectric material and the traces by FR-4 or PCB. The bumps can be surrounded by underfillers made of different plastic materials.

Figure 1.

Simulation and calculation sequence for one cycle (static) and more cycle (dynamic).

Concerning the mechanical stress calculation the interfaces play also an important rule. Due to the different materials and their material properties the stress in the structures can be high and can therefore have a strong influence on the reliability of the structures. The calculation of the stress gradients and stress migration mass flux divergence is strongly affected by the accuracy of the mesh.

3.1. Advantages of FEM compared with other applicable methods

Simulation methods like finite difference, finite element and Monte Carlo method are common in the electrical engineering. Monte Carlo method is described for the calculation of time to failure (TTF) distributions induced by electromigration [Huang]. Also in the scope of semiconductor simulation it is widely used [Jacoboni]. The complexity of the metallization and bump structures in combination with the question of multiphysic investigations makes this approach inconvenient in use. The application of the finite difference method is described in investigations of the intermetallic phase growth for instance [Chao]. A model for electromigration calculated by concentration gradients is described in [Joo] and the void evolution and motion was investigated by [Averbuch]. The mass transport along interface caused by electromigration using the level-set method is described in [Li]. All these approaches consider only a part of the thermal-electrical and mechanical behavior.

Calculations in the scope of the coupled thermal-electrical-mechanical behavior of metallization structures and bumps can be sufficiently carried out by commercial programs like ANSYS^®, ABAQUS or COMSOL using the finite element analyses. In the following investigation ANSYS^® is used as simulation tool.

3.2. Mesh density and singularities

In the thermal electrical calculation the degree of freedom is only the temperature. The mesh of the surrounding material depends in the investigated case of the metallization and due to this of the metallization mesh. A coarse mesh of the metallization can lead to an insufficient determination of the temperature in the isolation materials. Out of this the mesh of the metallization should be investigated with regard of an optimization.

The mesh can be easily refined by varying the count of the nodes. If a convergence concerning the change of the potential voltage is found the mesh is sufficient. A predication concerning the convergence of the current density at inhomogeneities like vias or edges in the metallization cannot be done. At these positions a strong dependence of the maximum current density due to current crowding in relation to the model of this inhomogeneity occurs. The layout from the EDA tools give at such positions rectangular edges or edges with a defined angle. Taking into account such an angle in the model may lead to insufficient calculation of the local current density at this position. Due to this the relation of the current crowding and angle of the edge will be investigated here.

Considering the two-dimensional case, choosing a right angle at the edge of a metallization a singularity will occur. In this case the current density can be calculated analytical with an infinity value [Betz]. The result of the numerical analysis of the current density at this position converges to a defined value due to the fact that the neighbor elements are taken into account for the calculation. If the right angle is replaced by a radius, the calculated current density depends on the selected shape of the radius, the element size and the element shape itself. For a determination of the real current density in the structure and the optimized radius, the edge of a simplified metallization consisting of aluminum, is investigated for the 2-dimensional and 3-dimensional case in sub-models. The boundaries in the calculations were taken of the coarse model.

In the 2-dimensional as well as in the 3-dimensional case the current crowding j_max/j_in increases with increasing element count cubically and the relation of the maximum current density to the applied current density j_max/j_in increases with decreasing element size exponentially. Out of this the numerical solution converges in the case of infinite small elements to the theoretical expected value. In the case of a curvature of the edges, the current crowding j_max/j_in has for increasing element counts as well as decreasing element size for every taken radius a constant value. Due to this the element count used in the simulations depends on the element size as well as the used radius in the model. A decreasing radius down to the range of an rectangular angle with constant element count leads to an increase of j_max/j_in.

Not in every case optical or scanning electron microscopy pictures are available for a determination of the structure shapes after the production process. Due to this for a determination of the optimized radius two dimensional simulations can be compared with results from the calculation of the maximum current density out of conformal mapping. The maximum current density and details about the homogeneity of the current density and resistance behavior at the investigated places can be achieved by the conformal mapping.

The calculation of the maximum current density is done by the following approach. In general for the analytical calculation of the current density distribution at the edge of an interconnection, a bent metallization with different width g and h and a radius r can be converted by two times conformal mapping into a flat conducting band with width π in the Z’ plane. In this flat conducting band a homogenous current density distribution can be assumed (figure 2). With the back transformation of this homogenous current density distribution into the Z plane the current density distribution along the path A-B-C can be achieved and due to this the maximum current density. With the parameters S=g/h as ratio of the metallization width and P=r/h as ratio of the radius r to the smaller metallization width h and the applied current density j_in in the vertical part of the metallization shall apply under the condition the P << 1 [Hagedorn]:

$$j_{\max }=1,{04}^{{\left[\frac{S^2+1}{P{S}^2}\right]}^{\frac{\text{1}}{\text{3}}}}{\text{ j}}_{\text{in}}$$E10

With equation 10 the current crowding j_max/j_in can be determined depending on the metallization width and the radius.

Figure 2.

Original figure of a two dimensional edge of an metallization in the Z-plane (left). Metallization after two times of conformal mapping in the Z’ plane (right).

In the 2-dimensional case a good correlation between the analytical solution and the maximum current density determined by the simulations can be found. In the 3-dimensional case the analytical solution is valid under the assumption of a homogenous current density distribution in the via, which means that the current density distribution can be mapped on the 2-dimensional case.

The maximum current density depending on the radius for a fixed element count is compared to the simulated values. The results are shown in figure 3. For a radius of 0.03µm a good correspondence between simulation and analytical solution is found. For a radius above 30nm the simulated values are caused by the increased element size above the analytical values. A radius below 30nm has an insufficient mesh and due to this the simulated values are beneath the analytical values. To get compliance between simulation and analytical solution the mesh has to be refined. This leads to an increase of the calculation time and is therefore not useful.

Figure 3.

Max. current density depending on the radius for the interconnect/via edge. Analytical solution done by conformal mapping and simulation.

3.3. Mechanical stress and stress gradients

In the finite element method simulation of the mechanical stress, the calculated unknown parameters are the displacement of the nodes. The different forces are calculated by the derivation of the displacements. The form function is continuously but not smoothed between the elements. Due to this the values at the nodes have to be averaged. Due to this and caused by the simplified form function and the limited element expansion the FEM is an approximation method. This can lead to failures in the calculation of the dilatations. When the results of an inaccurate mesh are derived the inaccuracy concerning a second derivation will be intensified. The correlation of the accuracy and mesh refinement is shown in figure 4 [Eichelseder]. From this knowledge it is clear that the distribution of the dilatation of the components S1 to S3 is more accurate than the strain distribution resp. the Von Mises Stress (VMS). This fact also implies that a simplified mesh or a coarse mesh can deliver accurate results concerning the dilatation. The derivate strain shows accurate values only after a refinement. After an additional refinement the gradient of the strain and out of this the hydrostatic stress (HS) can be calculated with a good accuracy. The calculation of the mass flux divergence due to stress migration includes also a derivation of the stress gradients, which will affect the accuracy additional. The accuracy of the finite element mesh has to be investigated before the mass flux divergence calculation is done. And the results have to be handled with care if the mesh is not adequate tested.

Figure 4.

Correlation of the accuracy and mesh refinement [Eichelseder].

4. Metallization

4.1. Metallization variation in the 90nm node

The investigated DD copper structure based on the dimension of the 90nm scaled down to the 65nm process technology. The model consists of SiCN as cap layer, Ta/TaN as barrier layer and different dielectrics like Silk™, Black Diamond II™ and SiCOH. In the mechanical calculations the process temperature using the birth and die algorithm in ANSYS are included, caused by the fact that the use of a reference temperature for the stress free state is not sufficient [Weide-Zaage 2008].

Figure 5.

Mesh of the investigated structure (left) and mass flux divergence vs metallization height for a variation of M1 and M2 (right).

The mesh of the investigated structure is shown in figure 5 (left). For simplification one via with two metal links right and left hand was generated. As one example for this proposal the metallization height here was varied. In the investigated cases for electromigration the current flux direction was calculated downstream from the second metallization M2 to the first metallization M1 and in the opposite direction upstream. The mass flux divergence of the variation is shown in figure 5 (right). It was found that the variation of the first metallization M1 and the variation of the second metallization M2 show a big influence on the mass flux divergence and out of this on the reliability of the structure. The weakest location found here is in and beneath the via.

Figure 6.

Maximum hydrostatic Stress for BDII™, Silk und SiCOH.

For technology nodes in the range of 500nm and below SiO₂ or combinations of USG (undoped-silicate-glass), PSG (phosphosilicated-glass) or FSG (fluorine-doped-silicate-glass) as a customary dielectric (IMD) is used. For products produced in nowadays low-κ dielectrics, with a low dielectric constant such as Silk ™, SiCOH, Black Diamond, or MSQ ™ II are used.

Figure 6 shows the maximum hydrostatic stress as a function of the applied current density for the different dielectric materials. The greatest mechanical stress occurs at a SiCOH dielectric and the lowest at Silk ™.

4.2. Separated migration mechanism with different via shape

A second model is based on dimensions of the 65nm technology node with a wide line and different via bottom geometries. The model is shown in figure 7 based on SEM pictures from the literature [Delsol, Lee]. Tantalum and TaN were chosen as barrier material, SiCN as capping material and SiCOH was chosen as dielectric material. The structures were investigated under test conditions with an applied current density of 1.5MA/cm² taken from [Lin2006a, Lin2006b]. The thermo-mechanical investigations were carried out considering the process-induced stress. The process temperatures are given in table 2.

The stress state in the copper metallization is strongly related to the process temperature of the copper. A high process temperature leads to a higher nearly stress free state of the whole structure compared to structures with lower temperatures [Matsuyama]. The nearly stress free state determined by investigating of the process induced stress is approx. 200°C. The calculation of the stress migration in the structures is described in [Weide-Zaage 2010].

Figure 7.

SEM pictures from [Lin2006, Lin2006(2)] (left) and Finite Element Mesh of the via region of the four different bottom geometries (right).

	Young Modul (GPa)	Poisson	CTE (300K)	Process Temperature (°C)
Cu	125	0.34	16.7 10^-6	200
TaN	185	0.33	6.6 10^-6	40
SiCN	100	0.17	3 10^-6	300
SiCOH	15	0.3	11.6 10^-6	350
Si	98	0.45	2.64 10^-6	25

Table 2.

Mechanical properties and process temperatures used in the simulations

4.2.1. Electromigration test temperature 325°C

The electromigration behavior of the four models was investigated with current (flux in both directions) and an activation energy of 0.9eV. The different via bottom geometries were verified by comparison with investigations from literature. The calculated maximum mass flux divergences as well as the reciprocal values which are related to the MTF are given in table 3. Comparing a Gouging and a Flat via and a Cone Shape and a V-Gouging via the Gouging via is more reliable than the Flat via and the Cone shape via more reliable than the V-Gouging via. The simulated models show that the Cone Shape via has the best results concerning the electromigration behavior.

Via-Shape	div EM	(div EM)^-1
Gouging Via	1.143 10^-3	875
Flat Via	0.917 10^-3	1091
V-Gouging Via	0.939 10^-3	1065
Cone Via	0.672 10^-3	1488

Table 3.

Mass flux divergence and reciprocal mass flux divergence.

The failure locations due to the void formation after the electromigration test in the literature are found downstream at the bottom of the via and upstream in and over the via at the interface metallization cap layer. A schematically overview of the different migration effects like electro-, thermo- and stress migration mass flux electron flux upstream and downstream and supposed void formation (yellow) is drawn in the via structure shown in figure 6. The electromigration mass flux is depending on the current flux direction and temperature gradients, the thermomigration mass flux depend on the temperature gradients and the stress migration on the stress gradients, which occur at tensile regions also under electromigration test conditions. An arbitrary assumed grain distribution is grey colored in this graphic. Thermomigration mass flux (blue) in the investigated structures proceeds independent of the current flux direction from first as well as the second metallization into the via. The electromigration mass flux (green) is at the upper part of the via downstream in the same direction and at the via bottom in opposite direction of the temperature gradients and thermomigration mass flux. In the upstream case it is vice versa. Due to the process induced stress distribution in the metallization, under electromigration test conditions, the metallization is mostly compressive with some tensile regions. The stress gradients in the tensile regions may lead to a current direction independent migration. High stress migration (orange) is found at the bottom and beneath the via above and below the barrier. The stress gradients show above the barrier upwards and below the barrier downwards. At the interface of the cap layer and the metallization they show down. Only the occurrence of a migration pathway leads to voiding. For interface migration the stress gradients have to have a component into the direction of the interface. Also the existence of a grain boundary pathway supports the voiding. Out of this foot voiding as well as voids in the via itself, both indicated in yellow can be explained by this (figure 8).

Figure 8.

Schematically electro-, thermo- and stress migration mass flux electron flux upstream and downstream and supposed void formation (yellow) in the via structure EM>SM>TM.

4.2.2. Thermomigration stress temperature 325°C

The thermomigration (Soret-Effect) was investigated under electromigration test conditions of 325°C. The calculation of the thermomigration shows high, values for the mass flux divergence. With temperature gradients of 50K/µm the gradients are high but not high enough to induce thermomigration. Gradients of 100K/µm are supposed to be an acceleration factor for the thermomigration. In figure 9 the mass flux divergence distribution due to thermomigration is shown. High values in the Gouging model are found at the via bottom in the via and at the corner of the metallization, barrier- and cap-layer. For the Cone Shape via high values are found only at the corner of the corner the metallization, barrier- and cap-layer. The activation energy for the thermomigration calculation is supposed to be too small and should be >1.2eV. Measurements of thermomigration activation energies of copper in copper are not known until now to verify this.

Figure 9.

Thermomigration mass flux divergence distribution for the Gouging and the Cone Shape model.

5. µ-Bump, CuSn-pillar and TSV

5.1. µ-Bump in comparison with BGA-PoP

A variation of the applied current in a Package-on-Package (PoP) bumps [Meinshausen 2010 (a)] and µ-bump [Meinshausen 2011] was carried out and the mass flux divergence distribution was determined. The bumps in the FE model of the PoP device with one bump chain consist of SAC305 (SnAgCu). For the under-bump metallization (UBM) and the surface finishes 6µm thick Ni layers were used. As mold compound (MC) around the upper contact surfaces "Stycast 1090" was chosen. The simulations were carried out with anisotropic and temperature depending material parameters. In the FR-4 substrates five layers of Cu traces and Cu vias between the top and the bottom package are placed to connect the top and the bottom bumps. The height of the upper, the middle and the lower copper layers is 20µm, 18µm and 36µm. The width of all layers is 250µm. The model of a µ-bump between two ICs is shown in figure 10 (left). The dimensions of the µ-bumps are similar to the test structures used in [Labie]. The diameter of the µ- bump is 25µm and the height is 10µm. Over and under the µ-bump a 100μm silicon layer resp. a 50μm thick silicon layer is representing the ICs of a CoC (Chip-on-Chip) structure. The ICs are covered with a 1µm thick Si₃N₄ passivation layer. The copper traces at the upper and the lower contact surface have a height of 0.5μm and a width of 32μm. The pitch is 40μm.

Figure 10.

Model of the µ-bump (left). The different materials are indicated by colors. The mass flux divergence vs. the applied current for SAC bump and µ-bump (right).

In figure 10 (right) the mass flux divergence depending on the applied current for the SAC PoP bump and the µ-bump is shown. Under the same applied current the mass flux divergence of µ-bumps is compared to PoP bumps about four orders of magnitude higher. Due to this fact the reliability of the µ-bump has a high electromigration risk. The weakest point of the µ-bump was found at the top of the bump due to current crowding combined with high temperature gradients.

Figure 11.

Temperature distribution in a BGA PoP Package with ICs (left) and maximum temperature in the bumps vs applied IC’s power loss (right).

[Meinshausen 2010(b), Meinshausen 2012] investigated in previous simulations the heat flux density of IC stacks with a constant surface load. The heat loss of the ICs is one reason for the existence of an inhomogeneous temperature distribution in packages. An applied power of 2-8W for all ICs in the IC stack was investigated as boundary condition for the simulations. The temperature distribution in the package for a power loss of 8W is shown in figure 11, left side. The maximum temperature occurs in the IC stack. The bump temperature is much lower (figure11, right). This leads to a strong inhomogeneous temperature distribution. On the other hand the temperature gradients in the bumps are not as high compared to the EM simulation. Due to this in this investigated case thermomigration will not occur.

5.2. TSV with µ-bump

Through silicon vias (TSV) have a wide range of applications in the modern packaging. Three different kinds of TSV can be identified which are common under the topic of 3D IC packaging, 3D IC integration and 3D Si integration. Due to this the processing as well as the size of the TSV is completely different. The application can be divided into three topics.

1. Package on Package stacking with TSV as connections between them

2. Chip to Chip, Chip to Wafer or Wafer to Wafer stacking with SiO₂ to Sio₂ and Cu to Cu bonding.

3. Wafer to Wafer stacking; Memory chip stacking TSV are used as connection between the chips; active or passive interposer.

In the case of the 3D IC integration the TSVs are necessary. Due to the fact that the data width is limited, the use of TSV with small sizes in the range of 5-10µm with a pitch in the range of 20-40µm a much wider I/O width is possible. The TSV shortens the way of the connections compared to the common wire bonding. Using the wire bonding the chip size has to be increased and due to this the costs increase. TSV as solution will solve this problem in future. For future 3D integration Cu-TSV as well as tungsten W-TSV will be used [Murugesan] caused by the fact that W is a capable material for sub µ-vias but not usable for power and ground applications due to its high specific resistance. Therefore Cu with a Ta barrier will be used.

Figure 12.

Schematic picture of a TSV and µ-bump (left) and FE-Model (right).

As an example in figure 12 (right) the mesh of a TSV with a SAC pillar is shown. The right side of figure 12 shows schematically a TSV with a µ-bump [Leduc]. The geometrical data for the simulations are taken from [Kitada, Lo]. The diameter of the TSV was set to 9µm and the height was set to 80µm. The copper barrier in the TSV consists of SiN and TiN. The thickness of the interconnections is 7µm. The passivation consists of SiN.

On top of the TSV high VMS was found caused by the fact that the copper in the TSV is under pressure and the material is pressed out of the TSV. Concerning the EM high mass fluxes were found at the bottom of the TSV to the metal trace. At this position also high thermal mass flux occurs. High hydrostatic stress occurs at the bottom of the TSV to the metal trace. The stress increases about 10% for an increase of the applied current from 0.1 to 0.3A.

5.3. CuSn-pillar

For a simplified assembly, a smaller pitch between the bumps and a low interconnect inductance CuSn-pillars can be used. Copper pillars are copper bumps with a thin layer of Sn on the top. These layer thicknesses can vary [Huffman, Syed]. The CuSn-pillars can be formed under pressure and a temperature load. This leads to the formation of Cu₃Sn and Cu₆Sn₅ phases. Due to this CuSn pillars with different Sn thickness and location in the bump were investigated. The bumps had a diameter of 24µm and a complete high of 45µm. Above the bumps a Cu trace and below the bumps a TSV (figure 12) were placed in the model. The applied current was set to 175 mA, the substrate temperature was varied from -50 to 150°C and the stress free temperature was set to 135°C. In figure 13 a part of the mesh and the temperature gradient distribution in the different bumps is shown.

Figure 13.

Mesh of the CuSn Pillar with Cu in magenta and Sn in yellow (above) and temperature gradient distribution for the different models (below).

The homogenous temperature distribution the bumps can be achieved by a placement of the Sn in the middle (model A). In the case of model B-C high temperature gradients occur above the Sn near the copper trace. At this position also current crowding occurs. Both can lead to a weak link at this position. Depending on the current flow direction the flux will be increased or decreased. In figure 14 the hydrostatic stress depending of the substrate temperature is shown for model A the highest stress occurs for temperatures about -50°C near the reference temperature of the stress free state the stress is the lowest. The processing temperatures should be included in the model.

	HS [MPa]	Tmax [K]	Div Jmax [a.u.]
B	174	396,86	0,486e12
C	175	396,54	0,459e12
D	172	397,02	0,512e12

Table 4.

Hydrostatic Stress, Maximum Temperature and Maximum Mass flux Divergence.

Figure 14.

Hydrostatic stress depending of the substrate temperature (model A).

6. Conclusion

The finite element analysis of metallization structures or packages and bumps provide an insight into local temperatures, current density and stress gradient distributions with the possibility of mass flux divergence calculation. A prediction of weak links in the structures helps to increase the reliability during the design phase. Also the costs will be decreased and redesigns avoided.

In future phase separation in the CuSn bumps by IMC growth as well as the influence of special placed grain boundaries and concentration gradients should be included in the model. The process induced stress should be also included in the package modeling.