Instruction FACTA UNIVERSITATIS Series:Electronics and Energetics Vol. 27, N o 1, March 2014, pp. 1 - 11 DOI: 10.2298/FUEE1401001C MICROSTRUCTURAL IMPACT ON ELECTROMIGRATION: A TCAD STUDY  Hajdin Ceric 1,2 , Roberto Lacerda de Orio 2 , Wolfhard H. Zisser 1,2 , Siegfried Selberherr 2 1 Christian Doppler Laboratory for Reliability Issues in Microelectronics at the Institute for Microelectronics, TU Wien, Austria 2 Institute for Microelectronics, TU Wien, Gußhausstraße 27–29, A-1040 Wien, Austria Abstract. Current electromigration models used for simulation and analysis of interconnect reliability lack the appropriate description of metal microstructure and consequently have a very limited predictive capability. Therefore, the main objective of our work was obtaining more sophisticated electromigration models. The problem is addressed through a combination of different levels of atomistic modeling and already available continuum level macroscopic models. A novel method for an ab initio calculation of the effective valence for electromigration is presented and its application on the analysis of EM behavior is demonstrated. Additionally, a simple analytical model for the early electromigration lifetime is obtained. We have shown that its application gives a reasonable estimate for the early electromigration failures including the effect of microstructure. Keywords: electromigration, interconnect, reliability, physical modeling, simulation 1. INTRODUCTION Electromigration (EM) experiments indicate that the copper interconnect lifetime decreases with every new interconnect generation. In particular, fast diffusivity paths cause a significant variation in the interconnect performance and EM degradation [1]. In order to produce more reliable interconnects, the fast diffusivity paths must be addressed when introducing new designs and materials. The EM lifetime depends on a variation of material properties at the microscopic and atomistic levels. Microscopic properties are grain boundaries and grains with their crystal orientation [2]. Atomistic properties are configurations of atoms at the grain boundaries, at the interfaces to the surrounding layers, and at the cross-section between grain boundaries and interfaces. Modern Technology Computer-Aided Design (TCAD) tools, in order to meet the challenges of contemporary interconnects, must cover two major  Received December 16, 2013 Corresponding author: Hajdin Ceric Christian Doppler Laboratory for Reliability Issues in Microelectronics at the Institute for Microelectronics, TU Wien, Austria (e-mail: Ceric@iue.tuwien.ac.at) 2 H. CERIC, R. LACERDA DE ORIO, W. ZISSER, S. SELBERHERR areas: physically based continuum-level modeling and first-principle/atomistic-level modeling. We present a computationally efficient ab initio method for calculation of the effective valence for EM and the atomistic EM force. The results of these ab initio calculations are applied for parameterization of a continuum-level model [7] and for simulation of the impact of the copper microstructure on the EM behavior. Additionally, an application of the kinetic Monte Carlo method in combination with the ab initio method for EM analysis is demonstrated. Results of ab initio and atomistic calculations are also used for the derivation of a compact model for early EM failures in copper dual-damascene M1/via structures. The model is based on the combination of a complete void nucleation model together with a simple mechanism of slit void growth under the via. It is demonstrated that the early EM lifetime is well described by a simple analytical expression, from where its statistical distribution can be obtained. Moreover, it is shown that the simulation results provide a reasonable estimate for the EM lifetimes. 2. THEORETICAL BACKGROUND 2.1. Electronic density based calculation of effective valence Generally, the effective valence is a tensor field ( ̅), which defines a linear relationship between the EM force ( ⃗) and an external electric field ( ⃗⃗). ⃗⃗( ⃗⃗⃗) ̅ ⃗⃗⃗ ⃗⃗⃗ (1) For the calculation of the effective valence several methods have been proposed, all of them being based on the computation of electron scattering states [3]. Density functional theory (DFT), in connection with the augmented plane wave (APW) method [4] or the Korringa-Kohn-Rostoker (KKR) method [5], has been established as the most powerful method for the determination of scattering states, however, it requires a demanding computational scheme. The cumbersome representation of scattering wave functions with many parameters is a heavy burden on stability and accuracy of subsequent numerical steps. In this work we introduce a more robust and efficient method to calculate the effective valence, which relies only on the electron density ⃗⃗ ⃗ . The basic idea is given in the following equations for the tensor components: ( ⃗⃗⃗) ∭ ⃗⃗ ( ⃗⃗⃗) ( ⃗⃗)[ ⃗⃗( ⃗⃗) ̂ ] ∭ ⃗ ( ⃗⃗ ⃗)[ ⃗⃗⃗ ( ⃗⃗⃗ ⃗) ̂ ] is the interaction potential between an electron and the migrating atom, ( ⃗⃗) is the relaxationtime due to scattering by phonons, ⃗( ⃗⃗) is the electron group velocity, and is the volume of a unit cell. The first integration is over the k-space and the second over the volume of the crystal. For the calculation of the electron density the DFT tool VASP [6] is used. An example of a VASP calculation is given in Fig. 1. Microstructural Impact on Electromigration: A TCAD Study 3 The electron density alone provides a qualitative explanation for the fact that the effective valence is higher in the bulk than in the grain boundaries. Similar analyses can be performed for atomic structures of different copper/insulator interfaces. Higher electron densities lead to higher effective valences, as can be seen from (2) [7]. For an accurate electron density calculation it is necessary to know the exact positions of the atoms in the structure. Fig. 1 Portion of the bulk copper crystal. The electron density is represented in two orthogonal planes. It varies from higher values (circle regions around atoms) closer to the atomic nucleus to lower in the inter-atomic space 2.2. Kinetic Monte Carlo simulation of electromigration To utilize results of quantum mechanical calculations for kinetic Monte Carlo simulations an average driving force along the diffusion jump path must be calculated. In general, the microscopic force-field depends on the position of the defect along the diffusion jump-path. The average of the microscopic force over the j-th diffusion jump path between locations ⃗ and ⃗ [3] is ⃗⃗ ⃗⃗ ∫ ⃗⃗ ⃗ ⃗⃗ ⃗⃗ ⃗ (3) The change in diffusion barrier height is equal to the net work by the microscopic force as the defect is moved from the initial to final sites over the entire jump path. The rates of defect jumps were calculated using the harmonic approximation to transition state theory (TST) [9]. In this approximation the transition rate is given by (4) is the migration energy (barrier) defined as the difference in energy between the transition state and the initial state, and is an attempt frequency [10]. For each defect site α the residence time is calculated as [11] 4 H. CERIC, R. LACERDA DE ORIO, W. ZISSER, S. SELBERHERR ∑ (5) is the number of possible jump sites from the site α. A single point defect is created at an arbitrary site, the clock is set to zero, and the defect is released to walk through the system. At each step, the jump direction is decided by a random number according to the local jump probabilities (6) The jump is implemented by updating the coordinates of the defect. By repeating the described random walk procedure for millions of defects, their concentration dependence on the effective valence tensor and the external field is calculated. 2.3. Compact model for lifetime estimation In order to calculate the mechanical stress in a three-dimensional copper dual damascene interconnect structure, a complex physically based model including the EM equation, the electro-thermal equation, and the mechanical equations has to be solved [7]. Korhonen et al. [14] proposed a simple one-dimensional model, where the solution for the stress at the cathode of a semi-infinite line is given by √ √ (7) Da is the effective atomic diffusivity and B is the effective modulus, which depends on the metal and the surrounding materials. Void formation occurs as soon as the mechanical stress reaches a critical magnitude at a site of weak adhesion, typically at the copper/capping layer interface [15], [16]. Thus, the void nucleation time is determined by the condition σ(tn)= σc, which applied to (7) yields ( ) (8) Where is the critical stress.The solution given by (8) is a good approximation to the more complete solutionobtained by solving a full physical model [7], [13] numerically, as will be shown later. It should be pointed out that (8) is valid as long as the stress remains significantly smaller than the stress magnitude at the steady state condition, which holds true for the void formation phase. Fig. 2 Early failure mode: slit void growth under the via Microstructural Impact on Electromigration: A TCAD Study 5 2.4. Void growth For a copper dual-damascene M1/via structure with downstream electron flow, EM failure analyses [11] indicate that the early failures are caused by slit voids located under the via, as shown in Fig. 2. Since the void is very thin and does not grow through the line height, void growth can be described by a one-dimensional process, so that the void length is given by (9) where is the drift velocity of the right edge of the void. The atomic flux into the right edge of the void is governed by the diffusivity of the copper/barrier layer interface , while the outgoing flux is governed by the surface diffusivity . Since , using the Nernst-Einstein equation one can write [17] (10) The EM failure occurs, when the void spans the via size, , so that the void growth time contribution to the EM lifetime is given by (11) 3. RESULTS AND DISCUSSION The ab initio method described above is applied for the calculation of the effective valence inside grain boundaries and the calculated value is used to parameterize our continuum-level model [7]. Prior to carrying out the ab initio calculation it is necessary to construct grain boundaries with exact positions of atoms. For this purpose an in-house molecular dynamic (MD) simulator with a many-atom interatomic potential based on effective-medium theory [8] is used. The total energy of the system is expressed as ∑ ∑ ∑ ( ) (12) Fig. 3 Formation of grain boundaries (circled regions) for a N-atom system, where V(rij) describes a pair potential and F(ni) describes the energy due to the electron density. An example of the construction of grain boundaries by means of MD simulation is presented in Fig. 3. 6 H. CERIC, R. LACERDA DE ORIO, W. ZISSER, S. SELBERHERR Ab initio calculations of the effective valence in copper grain boundaries have provided a value 75% lower than in the bulk for 4.3 eV Fermi energy (cf. Fig. 4), which is in good agreement with the results of Sorbello [3]. Along with the determination of the effective valence, ab initio calculations predict a lowering of the energy barrier for atomistic transport. Knowing the influence of the EM force on the diffusional barrier we utilize kinetic Monte Carlo [9] simulations for EM, which provide a closer look into the distribution of atoms in the presence of EM for a specific atomistic configuration. The dependence of the atomic concentration on the angle between the EM force and the jump direction is displayed in Fig. 5. The EM intensity clearly reduces from θ = 0 ◦ , where the EM force acts in the fast diffusivity path direction, to a minimum for θ = 90 ◦ , where the EM force is orthogonal to this direction. Ab inito calculations serve as basis to give a proper consideration of fast diffusivity paths and microstructure in the comprehensive physically based model [7].The solution of such a model is indeed rather complex and a detailed description of the numerical approach can be found in [13]. Fig. 4 Average distribution of the effective valence near a grain boundary. The external electric field is oriented parallel to the grain boundary Fig. 5 Concentration difference at four different angles () between the EM force and the atom migration paths Microstructural Impact on Electromigration: A TCAD Study 7 Fig. 6 shows the mechanical stress close to the via at the cathode end of a simulated line. A high stress develops adjacent to the via, where there is a line of intersection between the copper, the capping layer, and the barrier layer. For a copper dual-damascene M1/via structure with downstream electron flow, this is the typical site for void formation and growth leading to early EM failures. Since EM failure has a statistical character, in order to obtain a distribution of void nucleation times several lines with different microstructures were simulated. In particular, the mechanical stress under the via was monitored for a total of twenty lines, from where the resulting stress build-up for five different structures is shown in Fig. 7. We have observed that the time evolution of the stress curves can be divided into two main parts. In the first one the stress increases linearly with time, while in the second part it increases with the square root of time, as shown in Fig. 8 for a typical stress curve. It should be pointed out that Kirchheim [18] derived a linear stress increase from a one- dimensional version of a full physical model [7] under the condition that the stress is sufficiently low. In turn, Korhonen et al. [14] obtained a square root stress increase, as given by (7), from the solution of a simplified model for EM stress buildup. Thus, the stress build-up obtained from our numerical simulations with a rather complete model and for fully three-dimensional structures can be conveniently described by simple analytical solutions. Since void nucleation is expected to occur at high stress magnitudes, the second part of the stress curve shown in Fig. 8 is fitted by the square root model given in (7), where a is used as fitting parameter. By fitting the stress curves of all simulated structures, the distribution of the parameter a is determined, as shown in Fig. 9. The parameter is well described by lognormal statistics, where the mean and the standard deviation are MPa/s 1/2 and , respectively. Once a is known, the void formation time is obtained from (8). Since the distribution of a is also determined, we are able to obtain the statistical distribution of the void formation times, shown in Fig. 10. Due to the lognormal statistics of a, also follows a lognormal distribution, where the mean and standard deviation are h and . It should be pointed out that Filippi et al. [12] estimated a nucleation time of approximately 5h, which lies within the range predicted by the simulations. Fig. 6 Hydrostatic stress distribution (in MPa). High stress develops at the copper/capping/barrier layer intersection adjacent to the via 8 H. CERIC, R. LACERDA DE ORIO, W. ZISSER, S. SELBERHERR The void growth time is determined by (11), which is a function of the surface diffusivity. Choi et al. [17] obtained activation energy for surface diffusivity of eV on clean copper surfaces. It is expected that their measurement delivers a more precise copper surface diffusivity than the typical ones obtained on oxidized surfaces [17] and, therefore, we have used their estimate in our simulations. Furthermore, we have assumed that the activation energy follows a normal distribution [19]. As a consequence, both the surface diffusivity and the void growth time are lognormally distributed. The mean and the standard deviation of the void growth time distribution are h and , respectively. The void formation and the void growth times are of about the same order of magnitude, as shown in Fig. 10, which highlights the importance of considering both contributions for the early EM lifetime estimation under accelerated test conditions. Fig. 7 Stress build-up at the copper/capping/barrier layer intersection for lines with different microstructures Fig. 8 Fitting of a numerical solution using a linear and a square root model Microstructural Impact on Electromigration: A TCAD Study 9 As the void nucleation and the void growth times are known, the early EM lifetime is given by the combination of (8) and (11), ( ) (13) The distributions of the EM lifetimes are shown in Fig. 10, together with the experimental results obtained from Filippi et al. [12]. The lognormal mean and standard deviation of the simulated lifetimes are ̅ h and , We can see that the simulation results provide a reasonable description for the early EM lifetimes. A major advantage of (13) is that it is a simple analytical formula which is more rigorously related to the physical mechanisms active during the early EM failure development than Black's equation. A critical issue arises, however, with regard to the estimation of the parameter a. This parameter is affected by several factors, like diffusion coefficients, effective valence, mechanical moduli, microstructure, and more, so that it cannot be defined in a closed form in full physical modeling [7], [13]. Nevertheless, we have observed that it can be related to Korhonen's solution. In this way, it can be directly described by an analytical expression and connected to physical parameters according to (7). Fig. 9 Distribution of the square root model fitting parameter. The line represents a lognormal fit The relative difference between the simulated and experimental lifetimes for the same failure percentile varies between 15% and 20%, as shown in Fig. 11. The difference is smaller for shorter lifetimes, since the proposed slit void growth model is more accurate for very early failures, where the void volumes are smaller. Such an error magnitude is reasonable, given the required assumptions for the parameters and considering the simplicity of the model. 10 H. CERIC, R. LACERDA DE ORIO, W. ZISSER, S. SELBERHERR Fig. 10 Early EM lifetime distribution Fig. 11 Error between the simulation and the experimental results 4. CONCLUSION Our work demonstrates a novel approach for the calculation of the EM force on an atomistic level and its application to continuum-level modeling. The consideration of the accurate effective valence in grain boundaries allows a realistic simulation of EM behavior. The presented combination of atomistic force calculations with a kinetic Monte Carlo simulation enables sophisticated analyses of vacancy dynamics. A compact model for estimation of the early EM lifetimes in M1/via structures of copper dual-damascene interconnects was developed. The model was derived through the combination of a complete model for void nucleation together with a simple slit void growth mechanism under the via. Given the simplifications and assumptions made for the simulations, a reasonable approximation to experimental early EM failures has been obtained. Microstructural Impact on Electromigration: A TCAD Study 11 Acknowledgment: This work was partly supported by the Austrian Science Fund FWF, project P23296-N13. REFERENCES [1] Z.-S. Choi, R. Mönig, and C. V. Thompson, "Dependence of the Electromigration Flux on the Crystallographic Orientations of Different Grains in Polycrystalline Copper Interconnects," Appl. Phys. Lett., vol. 90, p. 241913, 2007. [2] E. Zschech and P. R. Besser, "Microstructure Characterization of Metal Interconnects and Barrier Layers: Status and Future," Proc. Interconnect Technol. Conf., pp. 233-330, 2000. [3] R. S. Sorbello, "Microscopic Driving Forces for Electromigration," in Materials Reliability Issues in Microelectronics, edited by J. R. Lloyd, F. G. Yost, and P. S. Ho, vol. 225 pp. 3-10, 1996. [4] R. P. Gupta, "Theory of Electromigration in Noble and Transition Metals," Phys. Rev. B, vol. 25, pp. 5118-5196, 1982. [5] D. N. Bly and P. J. Rous, "Theoretical Study of the Electromigration Wind Force for Adatom Migration at Metal Surfaces," Phys. Rev. B, vol. 53, pp. 13909, 2006. [6] G. Kresse and J. Furthmüller, "Efficient Iterative Schemes for ab initio Total-Energy Calculations Using a Plane-Wave Basis Set," Phys. Rev. B, vol. 54, pp. 11169, 1996. [7] H. Ceric, R. L. de Orio, J. Cervenka, and S. Selberherr, "A Comprehensive TCAD Approach for Assessing Electromigration Reliability of Modern Interconnects," IEEE Trans. Dev. Mat.Rel., vol. 9, pp. 9,2009. [8] K. W. Jacobsen, J. K. Norskov, and M. J. Puska, "Interatomic Interactions in the Effective-Medium Theory," Phys. Rev. B, vol. 35, pp. 7423, 1987. [9] R. Sorensen, Y. Mishin, and A. F. Voter, "Diffusion Mechanisms in Cu Grain Boundaries," Phys. Rev. B, vol. 62, pp. 3658, 2000. [10] M. Gall, C. Capasso, D. Jawarani, R. Hernandez, H. Kawasaki, and P. S. Ho, "Statistical Analysis of Early Failures in Electromigration," J. Appl. Phys., vol. 90, no. 2, pp. 732-740, 2001. [11] A. S. Oates and M. H. Lin, "Electromigration Failure Distribution of Cu/Low-k Dual-Damascene Vias: Impact of the Critical Current Density and a New Reliability Extrapolation Methodology," IEEE Trans. Device Mater. Rel., vol. 9, no. 2, pp. 244-254, 2009. [12] R. G. Filippi, P.-C.Wang, A. Brendler, P. S. McLaughlin, J. Poulin, B. Redder, and J. R. Lloyd, "The Effect of a Threshold Failure Time and Bimodal Behavior on the Electromigration Lifetime of Copper [13] Interconnects," Proc.Intl. Reliability Physics Symp., pp. 444-451, 2009. [14] R. L. de Orio, Dissertation, Technische Universität Wien, (2010). [Online]. Available: http://www.iue. tuwien.ac.at/phd/orio/ [15] M. A. Korhonen, P. Borgesen, K. N. Tu, and C.-Y. Li, J. "Stress Evolution due to Electromigration in Confined Metal Lines," Appl. Phys., vol. 73, no. 8, pp. 3790-3799, 1993. [16] R. J. Gleixner, B. M. Clemens, and W. D. Nix, "Void Nucleati on in Passivated Interconnect Lines: Effects of Site Geometries, Interfaces, and Interface Flaws," J. Mater. Res., vol. 12, pp. 2081-2090, 1997. [17] M. W. Lane, E. G. Liniger, and J. R. Lloyd, "Relationship Between Interfacial Adhesion and Electromigration in Cu Metallization," J. Appl. Phys., vol. 93, no. 3, pp. 1417-1421, 2003. [18] Z. S. Choi, R. Mönig, and C. V. Thompson, "Activation Energy and Prefactor for Surface Electromigration and Void Drift in Cu Interconnects," J. Appl. Phys., vol. 102, p. 083509, 2007. [19] R. Kirchheim, "Stress and Electromigration in Al-Lines of IntegratedCircuits," Acta Metall. Mater., vol. 40, no. 2, pp. 309-323, 1992. [20] L. Doyen, X. Federspiel, L. Arnaud, F. Terrier, Y. Wouters, and V. Girault, "Electromigration Multistress Pattern Technique for Copper Drift Velocity and Black's Parameters Extraction," Proc. Intl. Integrated [21] Reliability Workshop, pp. 74-78, 2007.