PROGRESS IN THE DEVELOPMENT OF AN ELECTRO-MIGRATION MODEL-

In this work, previous simulation works on electro-migration have been reviewed, and a multi-physics EM simulation method that combines electric, thermal, atomic diffusion and stress analysis has been proposed. The method is able to solve vacancy concentration distribution with coupling electric current, temperature gradient, hydro-static stress gradient simultaneously and to mimic void formation in realistic structure of microelectronics devices.


INTRODUCTION
Electro-migration is the phenomenon of atomic migration under the influence of electric current, and is usually associated with thermal, and stress effects [1].It is widely accepted that high current density is the main cause that increases the probability of the collision between electrons and diffusing atoms which result in void and hillock formation at the cathode and the anode respectively.In recent years, high density packaging (HDP) in electronics manufacturing has been increasingly adopted to meet the needs of miniaturization and increasing functionality in electronic products.Moreover, higher current density and complexity of structure has also generated high temperature gradients and high stress gradients leading to thermo-migration and stress-migration respectively in conductors.As a result, the true causes of EM involve a multi-physical cross coupling relationship and are poorly understood and characterized.In the past few decades, many researchers have studied EM and more recently computer simulation has become an ever more important tool in this research area.In this paper, previous works on EM simulation have been reviewed, and a new multi-physics EM simulation method that combines electric, thermal, atomic diffusion and stress analysis has been described.This new method is based on the multi-physics software PHYSICA and the aim of the work is to simulate the EM evolution process in interconnects in a unified environment and to help engineers achieve EM-aware electronics designs efficiently.

PREVIOUS STUDIES REVIEW
In 1960s, Black proposed a very simple model which first approaches the expression for the mean time to failure (MTTF) of a metal line subjected to electro-migration.He considered the MTTF is inversely proportional to the product of density of conducting electrons (݊ ), the momentum transfer from the electrons to the ions ‫,)∆(‬ and the density of thermally activated ions (ܰ ) as Eq.(1).
Furthermore, he assumed that both the electron density as well as the momentum transferring is proportional to the current density j.
And following the Arrhenius equation, we have Eq.( 3).

‫ܨܶܶܯ‬ =
మ ‫ݔ݁‬ ( where A is a material and geometry dependent constant [2 -4], ‫ܧ‬ is the activation energy, T is the temperature, and k is Boltzmann's constant.However, many observations found that not all experimental observation match Eq. ( 4), but they could be matched by amending the current density exponent.Therefore, Black's equation was modified to Eq. (5).
Many studies have been carried out to determine the value of n.A study on bulk metals suggested that a value of n=1 is more appropriate [5].D'Heurle and Agarwala [6,7], through a series of experimental works, fitted the n between 2 and 3 [8].For a long time, there has been no theoretically trusted value, and some of previous experiments were challenged with regards to their thermal management.When a test is run with relatively low current density stressing and with the thermal effects decoupled, n is reportedly between 1 and 2, but a value of 2 is probably used most widely [9].
Thermo-migration (TM) is almost inextricably linked to EM In addition, the composition of solder materials is moving towards lead-free Sn-based alloys that necessitate more investigation on the EM reliability issue due to poor field data availability [10]; Chih Chen et al, reported that the temperature gradient in a solder bump can reach 900-1000 K/cm with the diameter of a flip-chip solder bump being reduced to 25μm [11,12] and that this temperature gradient is high enough to trigger thermo-migration (TM).A thermal analysis has been carried out on the software ANSYS and the result showed the contact site between Cu line and solder bump has the highest temperature is around 503 K (with ambient temperature 373 K) and the temperature will significantly decrease to approximately 410 K because of the lower current density and the lager volume to dissipate heat in solder bump [12].Chih Chen and et al. concluded that the Al trace generates most of the heat and therefore, the Al trace is the weakest site in the aspect of TM.A verification experiment was also carried out using an infrared microscopy.The results showed a SnAg solder bump with a 5-µm Cu / 3-µm Ni UBM can reach 2237 K/cm with 0.6 A, and the solder bump will likely melt with higher current stressing.After examining current density, Chih Chen et al also concluded that if the current density has a magnitude of 10 4 A/cm 2 the driving force of EM is of the same magnitude as TM, but the driving force of TM will rise significantly beyond the driving force of EM with further current density increase and might cause solder bump melting.
Moreover, Stress-migration (SM) is due to the development of large tensile and compressive stress region at cathode and anode respectively as a combined consequence of the atomic migration and also due to the different thermal expansion coefficients of a narrow metal line and passivation or inter-metal dielectric material [13].Hydro-static stress is proved to be the main driving force for void nucleation; regions with the most concentrated hydrostatic tensile stress are assumed to be the most probable sites for void nucleation [14,15].The hydro-static stress, σ HS , was defined as the average of σ x , σ y , σ z .The hydrostatic stress gradient was defined as the ratio of the difference between stress values of two consecutives nodes to the distance between the nodes.The magnitude of the hydrostatic stress gradient Many researchers have attempted to model the growth of voids based on the SM studies in line structures with varying degrees of success over the years.Sauter et al [16] has developed a model, from which the time to failure, t f , prediction due to void growth was found to follow the simple relation where σ is the stress in the line, d is the grain size, Q is the activation energy for grain boundary diffusion in aluminum, R is the universal gas constant, and T is the absolute temperature at test.The stress σ is given by ߪ = 0.26ሺ3ሻ‫ܶ∆ߙ∆ܤ‬ Where B is the bulk modulus of Aluminum, α ∆ is the difference in thermal expansion coeffi- cient between the line and the substrate and T ∆ is the difference between the dielectric deposition temperature and the test temperature.P. M. Igic [17] proposed a method for calculating the thermal stress in aluminum interconnections during processing of the multilevel structure.
Compared to previous studies [18], this technique allows the residual stresses from one processing step to be used as the initial conditions for a subsequent step by using a so called activating-deactivating strategy.
Since EM is caused by atomic/vacancy migration under the influence of electric current, convective diffusion algorithms can be employed to describe the EM process.Based on Fick's laws of diffusion, the governing equation for EM equation can be written as Eq.(8).
where ‫ܥ‬ ௩ is the vacancy concentration, ‫ܬ‬ ௩ is the total vacancy flux, r is the source/sink term.H owever, in many EM works, atomic concentration ‫ܥ‬ , rather than vacancy concentration ‫ܥ‬ ௩ was being used for simulating atomic movement.Because of the vacancy exchange mechanis m, at equilibrium, the relation between the vacancy diffusivity and atomic diffusivity is given by [23]: where ‫ܦ‬ is the diffusivity of atom.By using the Eq. ( 8), Kirsten and et al. [19] carried out a thermo-electrical-mechanical simulation of EM using ANSYS.The atomic flux from mechanical effect was also investigated but the atomic self-diffusion was not taken into account.The mass flux as well as mass flux divergence calculated.A steady state analysis was also implemented in this work to predict highest vacancy concentration site (weakest part of structure).
Another similar computer simulation on EM for solder joint was proposed by Yong Liu et al [20][21][22].In the simulation, elements with high vacancy concentration were 'killed' automatically to simulate void formation.F. Cacho et al [23] used the software package COMSOL to simulate atomic self-diffusion and stress evolution in 3D and 2D structures and predicted their MTTF for pre-set failure criterion in the terms of vacancy concentration.The failure criterion is arbitrary and moreover, the models could not select "dead elements" automatically.

MODELLING
In most computer simulation research work that has been done to date, not all of the physical processes which affect EM have been factored and this makes it difficult to compare simulation and experimental results.In this work, a closely coupled multi-physics modelling method has been proposed.It can be used to predict atomic concentration and void formation in metals where electrical, thermal, stress, and geometry factors are all taken into account.The schematic diagram Fig. 1 describes the methodology.The method has been implemented using the multi-physics software package PHYSICA [24], which is capable of solving fluid flow, heat transfer, electric field and general diffusion equations simultaneously.

GOVERNING EQUATIONS
The governing diffusion-convection equation for vacancy evolution is Eq. ( 8).The total atomic flux can be described as Eq.(10). where and ௗ is the vacancy flux caused by self-diffusion, is the flux caused by electric current effects, ௧ is the vacancy flux caused by thermal effects, ௦ is the vacancy flux caused by stress effects, ‫ܦ‬ ୴ is the diffusivity of atoms, ܼ * is the effective charge, e is the elementary charge, k is the Boltzmann's constant, T is the temperature, E is the electrical field, ܳ * is the heat of transport, f is the vacancy relaxation factor, ߗ is the atomic volume, ߪ ுௌ is the hydrostatic stress, ‫ܦ‬ ୴ is the pre-exponential factor and ‫ܧ‬ is the activation energy.The atomic flux ୟ = − ୴ .

MODEL VALIDATION
To test the modeling framework, a known 1D analytic solution of Eq. ( 8) and Eq. ( 10) (without the thermal and stress terms) as detailed by R.L. de Orio and his colleagues [25] was compared with the simulation results.The boundary conditions are shown in Eq. ( 12) which states that perfect blocking boundary conditions are imposed at both ends of the 1D domain.
where L is the length of the line. where and and In the numerical simulation of this 1D problem, the domain has a length of 100 µm and it was divided into 20 elements with constant cross-section area.An electric potential difference was applied at the two ends of the domain generating a uniform electric filed and a constant current.The parameters in Table I are used.The EM drift velocity is 0.214 m/s if temperature or stress gradients are not taken into account.Fig. 2 shows that the analytical and simulation results for the vacancy concentration (C v /C v0 ) distribution are almost identical.

BLECH STRUCTURE
The Blech structure (Fig. 4) has been modelled using the methods described above.The parameters in Table II were used in the modelling and a voltage load of 0.3 V was ap-    From Table III, we found that the maximum flux divergence of EM (J em ) is significantly greater TM (J th ) and SM (J s ) and all maximum vacancy flux divergences are located around the interfaces of materials.We also noted the maximum temperature gradient is just about 1570 K/cm (Fig. 7).This value just reaches the criterion value of thermo-migration [10,11,33].Therefore, at an ambient temperature of 295 K and in well ventilated conditions, the TM contributes the least to atomic migration compared to other effects on this structure.Future work will be focused on the modeling of the combined effects of atomic migration in conductors and to model the void nucleation with the help of experiment to determine accurate activation energy E a and diffusion coefficient D v , which significantly affect the modeling as discussed in our previous study [34].

CONCLUSION
In this work, we reviewed the previous modelling work on EM and proposed a multiphysics modelling framework using the software package PHYSICA.The processes of diffusion, electric current, temperature gradient and stress gradient in electro-migration can be solved simultaneously.The computational model has been tested by a known analytical solution for a one dimensional problem and results are highly consistent with theoretical values.The Blech structure also has been implemented for comparing the contribution of individual effects.

Figure 1 :
Figure 1: Flow Chart of the EM model

Figure 2 :
Figure 2: The analytical and simulation results for the vacancy concentration (C v /C v0 ) along the X-axis at different times
ambient temperature was assumed to be 295 K.The resulting current density, temperature and hydro-static stress distributions are shown Fig.5.The highest current density is about 6.18x105 A/cm 2 .The maximum temperature is 456K and it's located at the centres of tungsten parts.The highest hydrostatic stress magnitude is about 101 MPa which appears at the W/Al interfaces.The vacancy flux divergence ߘ • ‫ܬ‬ ௩ due to electro-migration, thermomigration and stress-migration are shown in TableIIIand Fig.6.

Figure 4
Figure 4 the schematic diagram of a Blech model and boundary conditions in the numerical simulation.

Figure 5 (
Figure 5 (a) Current density distribution, (b) Temperature and (c) Hydro-static stress distribution of Blech model

Figure 6 :
Figure 6: Snapshot of the distribution of divergence at t=2000s for the (a) electrical, (b) thermal, and (c) stress vacancy flux.

Table I :
Parameters used in the calculations

Table II :
Parameters used in Blech Structure

Table III .
Maximum atomic flux divergences due to EM, TM and SM with initial vacancy concentration C v0 =1