- Email: [email protected]
Continuum and atomistic modeling of electromechanically-induced failure of ductile metallic thin films
Computational Materials Science 23 (2002) 242–249 www.elsevier.com/locate/commatsci
Continuum and atomistic modeling of electromechanicallyinduced failure of ductile metallic thin films Dimitrios Maroudas *, M. Rauf Gungor Department of Chemical Engineering, University of California, Santa Barbara, CA 93106-5080, USA Accepted 1 June 2001
Abstract A multiscale modeling approach is presented for the analysis of electromechanically-induced void morphological evolution and failure in ductile metallic thin films, which are used for device interconnections in integrated circuits. Selfconsistent mesoscopic simulations of surface morphological evolution are combined with atomistic calculations of surface properties and molecular-dynamics (MD) simulations of plastic deformation mechanisms in the vicinity of void surfaces. Results are presented that demonstrate a coupled mode of surface instability that leads to formation of electromigration-induced slit-like features and stress-induced crack-like features on void surfaces. In addition, MD results are presented for the dislocation-mediated mechanism and the kinetics of void growth in ductile metallic systems subject to hydrostatic and biaxial tensile strains. The incorporation of MD-derived constitutive information for plastic deformation into the mesoscopic analysis and simulation framework also is discussed. Ó 2002 Elsevier Science B.V. All rights reserved. PACS: 66.30.Qa; 68.35.Ja; 05.70.Ln; 85.40.Qx Keywords: Models of non-equilibrium and non-linear phenomena; Thin-film failure; Surface diffusion; Electromigration; Void growth; Dislocation emission; Boundary-element method; Molecular-dynamics simulations
1. Introduction A major materials reliability concern in microelectronics is the failure of Al or Cu thin films which are used for device interconnections in integrated circuits [1–3]. Failure mechanisms are driven mainly by the mass transport phenomenon of electromigration [1–3], as well as thermome-
*
Corresponding author. Tel.: +1-805-893-7346; fax: +1-805893-4731. E-mail address: [email protected] (D. Maroudas).
chanical stresses which develop in metallization layers during thermal processing due to the thermal expansivity mismatch between the metallic films and the dielectric material layers that encapsulate them [4,5]. Mechanisms of stress relief in these metallic films include void nucleation, crack initiation, and delamination phenomena at metal/ dielectric interfaces. In particular, void dynamics subsequent to void nucleation mediates most of the commonly observed interconnect failure cases [2,3]: such dynamical phenomena include void growth, migration, and surface morphological evolution under surface electromigration conditions.
0927-0256/02/$ - see front matter Ó 2002 Elsevier Science B.V. All rights reserved. PII: S 0 9 2 7 - 0 2 5 6 ( 0 1 ) 0 0 2 2 9 - 4
D. Maroudas, M. Rauf Gungor / Computational Materials Science 23 (2002) 242–249
Recent mesoscopic-scale theoretical studies have elucidated a broad range of electromigrationinduced non-linear dynamics of transgranular void surfaces [6–19], especially the crucial role of surface mobility anisotropy in governing surface morphological instabilities [16–19]. Aided also by our understanding of surface instabilities in mechanically stressed solids [20,21], such theoretical studies have shed new light on the complex void non-linear dynamics under the combined action of electric fields and mechanical stresses [22–24]; most notably, it has been suggested based on self-consistent numerical simulations that electromigration-induced void morphological evolution can be stabilized by appropriate tailoring of the residual stress level in the metallic film [22–24]. In spite of these recent theoretical advances, however, the effects of plastic deformation mechanisms on the void evolution in such ductile thin films have not been addressed rigorously and their role in electromechanically-induced thin-film failure is not understood. Linking mesoscopic-scale dynamical simulations of void morphological evolution with atomistic simulations of plastic deformation phenomena provides an ideal framework for predictive modeling of ductile thin-film response to electric and stress fields. This is an important step toward developing predictive simulation tools for electromechanically-driven evolution of realistic thin-film microstructures and computational testing of interconnect reliability. In this paper, such a multiscale modeling strategy is presented that aims at predicting electromechanically-induced void dynamical phenomena. The relevant length scales range from 1 nm, typical of crack-like feature widths at void tips, to 0:1–1 lm, typical of thinfilm widths and grain sizes for the systems of interest here. The corresponding time scales range from a few ps, characteristic of dislocation emission from crack tips under high strain, to several 100 h, characteristic of thin-film failure due to electromigration-induced void propagation under high current densities. Special emphasis is placed on obtaining atomic-scale information for rigorous development of plasticity constitutive equations for mesoscopic-scale modeling of void evolution. Results are presented for the dynamics
243
of void growth in Cu under hydrostatic and biaxial straining and the incorporation of these results into the mesoscopic modeling is discussed.
2. Multiscale modeling strategy Our multiscale modeling strategy for simulating thin-film failure mechanisms mediated by void dynamics combines atomistic simulations with mesoscopic modeling and it is outlined diagrammatically in Fig. 1. At the heart of our modeling approach lies a mesoscopic formulation of void morphological evolution due to surface mass transport under the action of electric fields, as described in detail in [16–19] and [22–24]. The predictive capabilities of this mesoscopic framework are enhanced by atomistic calculations of surface and interface properties, where the interatomic interactions are described according to the embedded-atom method (EAM) [25]; of particular importance to the study of void morphological instabilities is the dependence on the void surface orientation of the diffusivity of surface atoms [16–19,22–24]. In addition, EAM-based moleculardynamics (MD) simulations are implemented using large simulation supercells in order to probe the atomic-scale phenomena that govern plastic deformation in the vicinity of voids in strained ductile metallic systems. Together with the mechanistic understanding, the MD simulations provide constitutive information for the local plastic displacement field, upl , at the void surface. Within this framework, void dynamics can be simulated by computing the evolution of the local displacement normal to the void surface, un ðtÞ, according to the equation oun ¼ Xrs Js þ u_ n;pl ; ot
ð1Þ
where X is the atomic volume, Js is the total mass flux on the void surface, rs is the surface gradient operator, and u_ n;pl is the time derivative of the component of upl normal to the void surface. Note that Eq. (1) is reduced into the usual continuity equation for u_ n;pl ¼ 0. Js is expressed by a Nernst–Einstein equation, where the driving force includes contributions from capillarity, stress, and
244
D. Maroudas, M. Rauf Gungor / Computational Materials Science 23 (2002) 242–249
Fig. 1. Diagrammatic outline of multiscale strategy for modeling of void-mediated failure mechanisms in ductile metallic thin films.
the applied electric field [16–19,22–24], i.e., Js consists of curvature- and stress-induced surface diffusion and surface electromigration. The corresponding anisotropic surface diffusivity as a function of surface orientation can be quantified based on atomistic simulations [16–19,22–24]. Clearly, Eq. (1) is coupled strongly with the distribution of the electric and deformation fields in the metallic film and must be solved self-consistently with Maxwell’s equations and Cauchy’s mechanical equilibrium equations; the former equations can be reduced into Laplace’s equation for the electrostatic potential, while the latter ones require an appropriate set of constitutive and kinematic relations for closure. Dimensional analysis of the governing equations, Eq. (1) and the coupled electrostatic and mechanical problems, determines a multidimensional parameter space for the systematic study of thin-film mechanical response and failure [16– 19,22–24]: such parametric studies are based on self-consistent dynamical simulations and can contribute significantly to the design of experiments toward improving interconnect reliability. Important dimensionless parameters express the
applied electric field strength and stress with respect to capillarity, the initial void size, surface diffusional anisotropy, and the applied stress tensor anisotropy [16–19,22–24]. For dynamical simulations at the mesoscopic scale, domain discretization methods, such as boundary-integral or finite-element methods, are implemented for computing electric and strain fields along the evolving voids coupled with moving-boundary propagation techniques for the integration of Eq. (1). Our numerical simulation methods in the limit u_ n;pl ¼ 0 have been described in detail in [16–19] and [22–24]. In summary, the electrostatic potential and the displacement vector on the domain boundary are computed using a symmetric Galerkin boundary-integral formulation [26]. The corresponding electric and stress fields are computed using the surface-derivative calculation method of [27]. rs Js is calculated using a centered finitedifference scheme and Eq. (1) is integrated by implementing an Adams–Bashforth algorithm and using appropriately fine time steps [16–19,22–24]. Our MD simulation methodology also has been described in detail in [28]. In brief, an EAM
D. Maroudas, M. Rauf Gungor / Computational Materials Science 23 (2002) 242–249
parametrization for Cu is used [25] to carry out isothermal–isostrain simulations in supercells that contain typically around 106 atoms. Here, we focus on the response of cylindrical voids that are placed in the middle of the supercell and extend throughout the supercell thickness along the [1 1 1] crystallographic direction; this choice of crystallographic orientation is motivated by the preferred h1 1 1i grain orientation in metallic thin-film interconnects. The initial configuration is taken to be a perfectly cylindrical hole with circular crosssection on the (1 1 1) plane. Two types of straining are examined: hydrostatic tensile straining employing supercells with periodic boundary conditions (PBCs) in all three directions and biaxial tensile straining using slab supercells with tractionfree surfaces normal to the axis of the cylindrical void and PBCs in the two directions parallel to the surface.
3. Results and discussion Mesoscopic void dynamics simulations by integrating Eq. (1) in the limit u_ n;pl ¼ 0 are valuable in analyzing void evolution at constant void volume and have explained various failure phenomena in metallic thin-film interconnects. Such a catastrophic failure mechanism, captured by selfconsistent mesoscopic simulation [22–24], is shown in Fig. 2 through the evolution of a transgranular
245
void in a h1 1 1i-oriented grain of a metallic thin film under hydrostatic tension and an electric field directed from left to right. The simulation parameters are characteristic of electromigration testing conditions in passivated interconnect lines [22–24]. Failure occurs as a result of a coupled mode of electromigration-induced and stressinduced surface morphological instability; by failure here we mean the propagation of the void all the way across the film. Specifically, under the action of the electric field, an initially semi-circular void becomes faceted and obtains a semi-hexagonal shape due to the six-fold symmetry of the grain orientation [16–19]. Subsequently, a facet selection process becomes operative and leads to formation of a wedge-like void shape through elimination of one facet; this is followed by destabilization of a remaining facet and formation of a fast propagating faceted slit. After the slit propagates a certain distance, a finer-scale crack-like feature forms at the slit tip, it propagates at a speed higher by about two orders of magnitude than the average void migration speed, and causes failure of the film. This type of response constitutes a novel mode of instability that was introduced in the literature in [22–24]. In addition, the formation of cusp-like features on void surfaces, assuming isotropic surface mobility, was first mentioned in the literature by Xia et al. [14]. The void morphological features of Fig. 2 resemble qualitatively experimentally observed void configurations [29];
Fig. 2. Evolution of a transgranular void in a h1 1 1i-oriented grain of a metallic thin film under hydrostatic tension and an electric field directed from left to right. See [22–24] for the values of the dimensionless parameters used in this simulation, which are typical of electromigration testing conditions in passivated interconnect lines.
246
D. Maroudas, M. Rauf Gungor / Computational Materials Science 23 (2002) 242–249
however, the experimentally observed crack-like features appear to be blunted, as one would expect for cracks propagating in a ductile material. In addition, the simulated crack-like features are narrower by at least one order of magnitude than the original slits, i.e., their widths are of nanometer scale. Thus, understanding atomic-scale mechanisms of plastic deformation and calculating the corresponding plastic displacement rates become particularly important for the quantitative analysis of thin-film failure. MD simulations provide ideal means for probing the dislocation-mediated mechanisms of plastic deformation around voids in the ductile systems of interest. A representative example of capturing such a mechanism of void growth is highlighted in
Fig. 3 for the biaxial straining of a thin Cu film at 100 K and an applied strain of 5%; for the simulation, a slab supercell is used that contains 631,260 atoms and has sizes of 411.5, 429.5 and along the ½1 1 0 , ½1 1 2 and [1 1 1] crystal43.8 A lographic directions, respectively, at its unstrained state with an initial void fraction of 3.8%. In Fig. 3, the traces of crystallographic slip are evident on the film surface after 6.17 ps of straining; during this period, the void has become faceted with the morphology of an almost hexagonal prism, consistent with the six-fold symmetry of the chosen grain orientation for an fcc lattice, and has grown to a volume fraction of 5.6%. A pair of perfect screw dislocations with opposite Burgers vectors is emitted at each corner of the faceted void surface;
Fig. 3. Atomic configuration at time t ¼ 6:17 ps from an MD simulation of void evolution in a Cu thin film at 100 K and an applied biaxial strain of 5%; at t ¼ 0 the void is a perfectly cylindrical hole at the center of simulation cell. The atoms are shaded according to the gray scale to show the distribution in the film of the instantaneous potential energy per atom, U/N.
D. Maroudas, M. Rauf Gungor / Computational Materials Science 23 (2002) 242–249
detailed characterization of all the dislocation pairs is given in [28]. The dislocations glide along h1 1 0i directions on f1 1 1g planes. The screw dislocation lines intersect the void and form steps on the void surface. Dislocation glide causes motion of the step pairs on the void surface, which contributes to the radially outward evolution of the void surface. Dislocation pair generation and propagation occurs repeatedly and is responsible for the step pattern on the film surface, where the dislocation lines end. The same mechanistic elements of strain relaxation due to void growth are observed in Cu systems under both hydrostatic and biaxial straining: over a range of temperature and applied strain, the MD simulations reveal that void growth is mediated by screw dislocation pair emission from the void surface and subsequent propagation. As a result, another important feature common to both modes of straining is the presence of a plastic zone that surrounds the void: outside of this zone, the strain is uniform corresponding to the applied strain level; this becomes evident, e.g., by moni-
247
toring the distribution of the atomic-level pressure [28]. Interestingly, such a plastic zone formation is consistent with the theoretical analysis by Langer and Lobkovsky [30] of plastic deformation near a circular hole based on a rate-and-state theory of plasticity. Another common feature also is revealed by monitoring the atomic squared displacement distribution: stress-enhanced diffusional activity is high on the void surface and highest at the corner regions of the faceted surface. In addition, the atomic mobility distribution is strongly correlated with the dislocation distribution in the plastic zone, which indicates that ‘‘pipe’’ diffusion also is operative along the dislocation cores. During the MD simulations, the computed trajectories of the void surface atoms can be used to derive the evolution of the local plastic displacements according to the dislocation-mediated mechanism of void growth. In addition, the kinetics of void growth can be derived directly by monitoring the void volume evolution. Such representative kinetic results for the void fraction evolution, f ðtÞ, are shown in Fig. 4 for both
Fig. 4. MD simulation results for void fraction evolution, f ðtÞ, in Cu at 100 K under (a) hydrostatic and (b) biaxial tensile strain of 4% (}) and 5% (). In both cases (a) and (b), the solid curves represent fits to the MD results according to Eq. (2) for strains of (1) 4% and (2) 5%.
248
D. Maroudas, M. Rauf Gungor / Computational Materials Science 23 (2002) 242–249
hydrostatic and biaxial straining, Figs. 4(a) and (b), respectively, for two levels of applied strain, 4% and 5%, at 100 K. In both cases, Figs. 4(a) and (b), the solid curves represent least-squares fits to the MD simulation results according to the same phenomenological equation: f ðtÞ ¼ f0 þ Að1 exp ½ ðt=sÞn Þ;
ð2Þ
where f0 ¼ f ðt ¼ 0Þ and A ¼ limt!1 ½f ðtÞ f0 , which represents the void fraction increase when the void volume reaches each steady state. Eq. (2) implies that limt!0 ½f ðtÞ f0 / tn , which is an important scaling relation for the initial stage of void growth. For each mode of straining, the exponent n in Eq. (2) has been found to be invariant over the entire range of applied tensile strain that was examined: from 1% to 8%. Throughout this strain range, n 4=3 for hydrostatic straining and n 3=2 for biaxial straining; for example, n ¼ 1:33 for both curves of Fig. 4(a), while n ¼ 1:51 for both curves of Fig. 4(b). In conjunction with the invariance of this exponent, n, for given mode of straining, it should be mentioned that the void growth mechanism of strain relaxation remains unchanged for each mode over the entire strain range examined. In Eq. (2), the relaxation time, s, depends on the applied strain level: for the results of Fig. 4(a), s ¼ 15:29 and 13.25 ps for strains of 4% and 5%, respectively, while for the results of Fig. 4(b), s ¼ 5:73 and 7.06 ps for strains of 4% and 5%, respectively. In the case of biaxial straining, the shorter relaxation time is understood because of the thin-film finite thickness and presence of the free surfaces. The analysis of MD trajectories to derive void growth kinetics, such as Eq. (2), is particularly promising for the role of MD simulations, as computer experiments, for development of constitutive plasticity relationships that can be incorporated directly into the mesoscopic dynamical framework, Eq. (1). Specifically, the void growth rate, that can be expressed analytically by the time derivative of the right-hand side of Eq. (2), also is expressed as the integral over the surface of the void of the local plastic displacement rate normal to the void surface. In the case of biaxial straining on the (1 1 1) plane of a thin film of thickness h, for example, this relationship is given by
df ðtÞ 6h ¼ dt Vcell
Z
p=3
Rðh; tÞu_ n;pl ðh; tÞ dh;
ð3Þ
0
where Vcell is the simulation supercell volume, R is the local radius of curvature of the cylindrical void surface, h is the corresponding polar angle, and the six-fold symmetric geometry of a cylindrical void with its axis along h1 1 1i has been taken into account to simplify the integral. Most elegantly, instead of directly using the MD simulation results, analytical constitutive theories for the plastic displacement rate can be constructed and serve as parametrization schemes; the corresponding parameters can be quantified by fitting the MD results. Such analytical theories are currently under development and will be reported in a forthcoming publication.
4. Conclusions and future directions A systematic multiscale modeling strategy has been presented for the analysis of electromechanically-induced void dynamics in ductile metallic thin films and void-mediated failure mechanisms in such metallic systems. This multiscale strategy links self-consistent dynamical simulations of void morphological evolution due to surface transport in a deforming solid with atomistic calculations of surface and interfacial properties and MD simulations of plastic deformation mechanisms in the vicinity of the void surface. An important theoretical conclusion of our study is that MD simulations of strain-induced void dynamics can be used to quantify rigorously parametrization schemes provided by mesoscopic or phenomenological theories of void growth and, thus, develop constitutive theories for plastic displacement rates; such constitutive equations can be incorporated directly into the mesoscopic void evolution equations. An additional important conclusion regarding thin-film failure is that, in the most general case, a coupled mode of failure is operative involving the formation of electromigration-induced slits and other morphological features and stress-induced crack-like features at the slit tips. However, the most challenging task ahead is
D. Maroudas, M. Rauf Gungor / Computational Materials Science 23 (2002) 242–249
the identification of possible regimes of electromechanical conditions where stress and electromigration compete to give rise to stable void behavior: a morphologically stable void surface that is either immobile or simply migrates along the film and is driven out of the film. Addressing this challenging problem motivates exciting future research directions for multiscale modeling in ductile metallic thin films, such as the Al and Cu films used in modern interconnect technologies. Such future studies include self-consistent dynamical simulations over a broad range of electromechanical conditions incorporating into the mesoscopic formulation MD-derived plastic response for the evolving void surface. This is a particularly demanding task that requires not only extensive mesoscopic parametric studies and nonlinear analyses but also systematic MD simulation schedules for fundamental mechanistic investigations of stress-induced void dynamics over a wide range of thermomechanical conditions. Multiscale theoretical studies within this framework are currently underway and will be reported in future publications. Acknowledgements This work was supported by the National Science Foundation through a CAREER Award (No. ECS-9501111) to one of the authors (D.M.), the Frontiers of Materials Science Program of the UCSB Materials Research Laboratory and Los Alamos National Laboratory (Award No. STBUC:97-63), and the Oak Ridge Institute for Science and Education through a Faculty Research Fellowship to one of the authors (D.M.). Fruitful discussions with L.J. Gray, S.J. Zhou, G.E. Beltz, and J.S. Langer also are gratefully acknowledged.
249
References [1] P.S. Ho, T. Kwok, Rep. Prog. Phys. 52 (1989) 301. [2] C.V. Thompson, J.R. Lloyd, MRS Bull. 18 (12) (1993) 19, and references therein. [3] C.-K. Hu, J.M.E. Harper, Mater. Chem. Phys. 52 (1998) 5. [4] A.I. Sauter, W.D. Nix, J. Mater. Res. 7 (1992) 1133. [5] A.F. Bower, L.B. Freund, J. Appl. Phys. 74 (1993) 3855. [6] Z. Suo, W. Wang, M. Yang, Appl. Phys. Lett. 64 (1994) 1944. [7] W.Q. Wang, Z. Suo, T.-H. Hao, J. Appl. Phys. 79 (1996) 2394. [8] D. Maroudas, Appl. Phys. Lett. 67 (1995) 798. [9] D. Maroudas, M.N. Enmark, C.M. Leibig, S.T. Pantelides, J. Comp.-Aided Mater. Des. 2 (1995) 231. [10] O. Kraft, E. Arzt, Appl. Phys. Lett. 66 (1995) 2063. [11] O. Kraft, E. Arzt, Acta Mater. 45 (1997) 1599. [12] M. Mahadevan, R.M. Bradley, J. Appl. Phys. 79 (1996) 6840. [13] M. Mahadevan, R.M. Bradley, Phys. Rev. B 59 (1999) 11037. [14] L. Xia, A.F. Bower, Z. Suo, C.F. Shih, J. Mech. Phys. Solids 45 (1997) 1473. [15] M. Schimschak, J. Krug, Phys. Rev. Lett. 80 (1998) 1674. [16] M.R. Gungor, D. Maroudas, Appl. Phys. Lett. 72 (1998) 3452. [17] M.R. Gungor, D. Maroudas, Surf. Sci. 415 (1998) L1055. [18] M.R. Gungor, D. Maroudas, J. Appl. Phys. 85 (1999) 2233. [19] M.R. Gungor, D. Maroudas, Surf. Sci. 461 (2000) L550. [20] see, e.g., W.H. Yang, D.J. Srolovitz, Phys. Rev. Lett. 71 (1993) 1593. [21] W.H. Yang, D.J. Srolovitz, J. Mech. Phys. Solids 42 (1994) 1551. [22] M.R. Gungor, D. Maroudas, L.J. Gray, Appl. Phys. Lett. 73 (1998) 3848. [23] M.R. Gungor, D. Maroudas, Surf. Sci. 432 (1999) L604. [24] M.R. Gungor, D. Maroudas, Int. J. Fract. 109 (2001) 47. [25] S.M. Foiles, M.I. Baskes, M.S. Daw, Phys. Rev. B 33 (1986) 7983. [26] L.J. Gray, D. Maroudas, M.N. Enmark, E.F. D’Azevedo, Eng. Anal. Boundary Elements 23 (1999) 267. [27] L.J. Gray, D. Maroudas, M.N. Enmark, Comp. Mech. 22 (1998) 187. [28] M.R. Gungor, D. Maroudas, S.J. Zhou, Appl. Phys. Lett. 77 (2000) 343. [29] Y.-C. Joo, C.V. Thompson, J. Appl. Phys. 81 (1997) 6062. [30] J.S. Langer, A.E. Lobkovsky, Phys. Rev. E 60 (1999) 6978.
Continuum and atomistic modeling of electromechanicallyinduced failure of ductile metallic thin films Dimitrios Maroudas *, M. Rauf Gungor Department of Chemical Engineering, University of California, Santa Barbara, CA 93106-5080, USA Accepted 1 June 2001
Abstract A multiscale modeling approach is presented for the analysis of electromechanically-induced void morphological evolution and failure in ductile metallic thin films, which are used for device interconnections in integrated circuits. Selfconsistent mesoscopic simulations of surface morphological evolution are combined with atomistic calculations of surface properties and molecular-dynamics (MD) simulations of plastic deformation mechanisms in the vicinity of void surfaces. Results are presented that demonstrate a coupled mode of surface instability that leads to formation of electromigration-induced slit-like features and stress-induced crack-like features on void surfaces. In addition, MD results are presented for the dislocation-mediated mechanism and the kinetics of void growth in ductile metallic systems subject to hydrostatic and biaxial tensile strains. The incorporation of MD-derived constitutive information for plastic deformation into the mesoscopic analysis and simulation framework also is discussed. Ó 2002 Elsevier Science B.V. All rights reserved. PACS: 66.30.Qa; 68.35.Ja; 05.70.Ln; 85.40.Qx Keywords: Models of non-equilibrium and non-linear phenomena; Thin-film failure; Surface diffusion; Electromigration; Void growth; Dislocation emission; Boundary-element method; Molecular-dynamics simulations
1. Introduction A major materials reliability concern in microelectronics is the failure of Al or Cu thin films which are used for device interconnections in integrated circuits [1–3]. Failure mechanisms are driven mainly by the mass transport phenomenon of electromigration [1–3], as well as thermome-
*
Corresponding author. Tel.: +1-805-893-7346; fax: +1-805893-4731. E-mail address: [email protected] (D. Maroudas).
chanical stresses which develop in metallization layers during thermal processing due to the thermal expansivity mismatch between the metallic films and the dielectric material layers that encapsulate them [4,5]. Mechanisms of stress relief in these metallic films include void nucleation, crack initiation, and delamination phenomena at metal/ dielectric interfaces. In particular, void dynamics subsequent to void nucleation mediates most of the commonly observed interconnect failure cases [2,3]: such dynamical phenomena include void growth, migration, and surface morphological evolution under surface electromigration conditions.
0927-0256/02/$ - see front matter Ó 2002 Elsevier Science B.V. All rights reserved. PII: S 0 9 2 7 - 0 2 5 6 ( 0 1 ) 0 0 2 2 9 - 4
D. Maroudas, M. Rauf Gungor / Computational Materials Science 23 (2002) 242–249
Recent mesoscopic-scale theoretical studies have elucidated a broad range of electromigrationinduced non-linear dynamics of transgranular void surfaces [6–19], especially the crucial role of surface mobility anisotropy in governing surface morphological instabilities [16–19]. Aided also by our understanding of surface instabilities in mechanically stressed solids [20,21], such theoretical studies have shed new light on the complex void non-linear dynamics under the combined action of electric fields and mechanical stresses [22–24]; most notably, it has been suggested based on self-consistent numerical simulations that electromigration-induced void morphological evolution can be stabilized by appropriate tailoring of the residual stress level in the metallic film [22–24]. In spite of these recent theoretical advances, however, the effects of plastic deformation mechanisms on the void evolution in such ductile thin films have not been addressed rigorously and their role in electromechanically-induced thin-film failure is not understood. Linking mesoscopic-scale dynamical simulations of void morphological evolution with atomistic simulations of plastic deformation phenomena provides an ideal framework for predictive modeling of ductile thin-film response to electric and stress fields. This is an important step toward developing predictive simulation tools for electromechanically-driven evolution of realistic thin-film microstructures and computational testing of interconnect reliability. In this paper, such a multiscale modeling strategy is presented that aims at predicting electromechanically-induced void dynamical phenomena. The relevant length scales range from 1 nm, typical of crack-like feature widths at void tips, to 0:1–1 lm, typical of thinfilm widths and grain sizes for the systems of interest here. The corresponding time scales range from a few ps, characteristic of dislocation emission from crack tips under high strain, to several 100 h, characteristic of thin-film failure due to electromigration-induced void propagation under high current densities. Special emphasis is placed on obtaining atomic-scale information for rigorous development of plasticity constitutive equations for mesoscopic-scale modeling of void evolution. Results are presented for the dynamics
243
of void growth in Cu under hydrostatic and biaxial straining and the incorporation of these results into the mesoscopic modeling is discussed.
2. Multiscale modeling strategy Our multiscale modeling strategy for simulating thin-film failure mechanisms mediated by void dynamics combines atomistic simulations with mesoscopic modeling and it is outlined diagrammatically in Fig. 1. At the heart of our modeling approach lies a mesoscopic formulation of void morphological evolution due to surface mass transport under the action of electric fields, as described in detail in [16–19] and [22–24]. The predictive capabilities of this mesoscopic framework are enhanced by atomistic calculations of surface and interface properties, where the interatomic interactions are described according to the embedded-atom method (EAM) [25]; of particular importance to the study of void morphological instabilities is the dependence on the void surface orientation of the diffusivity of surface atoms [16–19,22–24]. In addition, EAM-based moleculardynamics (MD) simulations are implemented using large simulation supercells in order to probe the atomic-scale phenomena that govern plastic deformation in the vicinity of voids in strained ductile metallic systems. Together with the mechanistic understanding, the MD simulations provide constitutive information for the local plastic displacement field, upl , at the void surface. Within this framework, void dynamics can be simulated by computing the evolution of the local displacement normal to the void surface, un ðtÞ, according to the equation oun ¼ Xrs Js þ u_ n;pl ; ot
ð1Þ
where X is the atomic volume, Js is the total mass flux on the void surface, rs is the surface gradient operator, and u_ n;pl is the time derivative of the component of upl normal to the void surface. Note that Eq. (1) is reduced into the usual continuity equation for u_ n;pl ¼ 0. Js is expressed by a Nernst–Einstein equation, where the driving force includes contributions from capillarity, stress, and
244
D. Maroudas, M. Rauf Gungor / Computational Materials Science 23 (2002) 242–249
Fig. 1. Diagrammatic outline of multiscale strategy for modeling of void-mediated failure mechanisms in ductile metallic thin films.
the applied electric field [16–19,22–24], i.e., Js consists of curvature- and stress-induced surface diffusion and surface electromigration. The corresponding anisotropic surface diffusivity as a function of surface orientation can be quantified based on atomistic simulations [16–19,22–24]. Clearly, Eq. (1) is coupled strongly with the distribution of the electric and deformation fields in the metallic film and must be solved self-consistently with Maxwell’s equations and Cauchy’s mechanical equilibrium equations; the former equations can be reduced into Laplace’s equation for the electrostatic potential, while the latter ones require an appropriate set of constitutive and kinematic relations for closure. Dimensional analysis of the governing equations, Eq. (1) and the coupled electrostatic and mechanical problems, determines a multidimensional parameter space for the systematic study of thin-film mechanical response and failure [16– 19,22–24]: such parametric studies are based on self-consistent dynamical simulations and can contribute significantly to the design of experiments toward improving interconnect reliability. Important dimensionless parameters express the
applied electric field strength and stress with respect to capillarity, the initial void size, surface diffusional anisotropy, and the applied stress tensor anisotropy [16–19,22–24]. For dynamical simulations at the mesoscopic scale, domain discretization methods, such as boundary-integral or finite-element methods, are implemented for computing electric and strain fields along the evolving voids coupled with moving-boundary propagation techniques for the integration of Eq. (1). Our numerical simulation methods in the limit u_ n;pl ¼ 0 have been described in detail in [16–19] and [22–24]. In summary, the electrostatic potential and the displacement vector on the domain boundary are computed using a symmetric Galerkin boundary-integral formulation [26]. The corresponding electric and stress fields are computed using the surface-derivative calculation method of [27]. rs Js is calculated using a centered finitedifference scheme and Eq. (1) is integrated by implementing an Adams–Bashforth algorithm and using appropriately fine time steps [16–19,22–24]. Our MD simulation methodology also has been described in detail in [28]. In brief, an EAM
D. Maroudas, M. Rauf Gungor / Computational Materials Science 23 (2002) 242–249
parametrization for Cu is used [25] to carry out isothermal–isostrain simulations in supercells that contain typically around 106 atoms. Here, we focus on the response of cylindrical voids that are placed in the middle of the supercell and extend throughout the supercell thickness along the [1 1 1] crystallographic direction; this choice of crystallographic orientation is motivated by the preferred h1 1 1i grain orientation in metallic thin-film interconnects. The initial configuration is taken to be a perfectly cylindrical hole with circular crosssection on the (1 1 1) plane. Two types of straining are examined: hydrostatic tensile straining employing supercells with periodic boundary conditions (PBCs) in all three directions and biaxial tensile straining using slab supercells with tractionfree surfaces normal to the axis of the cylindrical void and PBCs in the two directions parallel to the surface.
3. Results and discussion Mesoscopic void dynamics simulations by integrating Eq. (1) in the limit u_ n;pl ¼ 0 are valuable in analyzing void evolution at constant void volume and have explained various failure phenomena in metallic thin-film interconnects. Such a catastrophic failure mechanism, captured by selfconsistent mesoscopic simulation [22–24], is shown in Fig. 2 through the evolution of a transgranular
245
void in a h1 1 1i-oriented grain of a metallic thin film under hydrostatic tension and an electric field directed from left to right. The simulation parameters are characteristic of electromigration testing conditions in passivated interconnect lines [22–24]. Failure occurs as a result of a coupled mode of electromigration-induced and stressinduced surface morphological instability; by failure here we mean the propagation of the void all the way across the film. Specifically, under the action of the electric field, an initially semi-circular void becomes faceted and obtains a semi-hexagonal shape due to the six-fold symmetry of the grain orientation [16–19]. Subsequently, a facet selection process becomes operative and leads to formation of a wedge-like void shape through elimination of one facet; this is followed by destabilization of a remaining facet and formation of a fast propagating faceted slit. After the slit propagates a certain distance, a finer-scale crack-like feature forms at the slit tip, it propagates at a speed higher by about two orders of magnitude than the average void migration speed, and causes failure of the film. This type of response constitutes a novel mode of instability that was introduced in the literature in [22–24]. In addition, the formation of cusp-like features on void surfaces, assuming isotropic surface mobility, was first mentioned in the literature by Xia et al. [14]. The void morphological features of Fig. 2 resemble qualitatively experimentally observed void configurations [29];
Fig. 2. Evolution of a transgranular void in a h1 1 1i-oriented grain of a metallic thin film under hydrostatic tension and an electric field directed from left to right. See [22–24] for the values of the dimensionless parameters used in this simulation, which are typical of electromigration testing conditions in passivated interconnect lines.
246
D. Maroudas, M. Rauf Gungor / Computational Materials Science 23 (2002) 242–249
however, the experimentally observed crack-like features appear to be blunted, as one would expect for cracks propagating in a ductile material. In addition, the simulated crack-like features are narrower by at least one order of magnitude than the original slits, i.e., their widths are of nanometer scale. Thus, understanding atomic-scale mechanisms of plastic deformation and calculating the corresponding plastic displacement rates become particularly important for the quantitative analysis of thin-film failure. MD simulations provide ideal means for probing the dislocation-mediated mechanisms of plastic deformation around voids in the ductile systems of interest. A representative example of capturing such a mechanism of void growth is highlighted in
Fig. 3 for the biaxial straining of a thin Cu film at 100 K and an applied strain of 5%; for the simulation, a slab supercell is used that contains 631,260 atoms and has sizes of 411.5, 429.5 and along the ½1 1 0 , ½1 1 2 and [1 1 1] crystal43.8 A lographic directions, respectively, at its unstrained state with an initial void fraction of 3.8%. In Fig. 3, the traces of crystallographic slip are evident on the film surface after 6.17 ps of straining; during this period, the void has become faceted with the morphology of an almost hexagonal prism, consistent with the six-fold symmetry of the chosen grain orientation for an fcc lattice, and has grown to a volume fraction of 5.6%. A pair of perfect screw dislocations with opposite Burgers vectors is emitted at each corner of the faceted void surface;
Fig. 3. Atomic configuration at time t ¼ 6:17 ps from an MD simulation of void evolution in a Cu thin film at 100 K and an applied biaxial strain of 5%; at t ¼ 0 the void is a perfectly cylindrical hole at the center of simulation cell. The atoms are shaded according to the gray scale to show the distribution in the film of the instantaneous potential energy per atom, U/N.
D. Maroudas, M. Rauf Gungor / Computational Materials Science 23 (2002) 242–249
detailed characterization of all the dislocation pairs is given in [28]. The dislocations glide along h1 1 0i directions on f1 1 1g planes. The screw dislocation lines intersect the void and form steps on the void surface. Dislocation glide causes motion of the step pairs on the void surface, which contributes to the radially outward evolution of the void surface. Dislocation pair generation and propagation occurs repeatedly and is responsible for the step pattern on the film surface, where the dislocation lines end. The same mechanistic elements of strain relaxation due to void growth are observed in Cu systems under both hydrostatic and biaxial straining: over a range of temperature and applied strain, the MD simulations reveal that void growth is mediated by screw dislocation pair emission from the void surface and subsequent propagation. As a result, another important feature common to both modes of straining is the presence of a plastic zone that surrounds the void: outside of this zone, the strain is uniform corresponding to the applied strain level; this becomes evident, e.g., by moni-
247
toring the distribution of the atomic-level pressure [28]. Interestingly, such a plastic zone formation is consistent with the theoretical analysis by Langer and Lobkovsky [30] of plastic deformation near a circular hole based on a rate-and-state theory of plasticity. Another common feature also is revealed by monitoring the atomic squared displacement distribution: stress-enhanced diffusional activity is high on the void surface and highest at the corner regions of the faceted surface. In addition, the atomic mobility distribution is strongly correlated with the dislocation distribution in the plastic zone, which indicates that ‘‘pipe’’ diffusion also is operative along the dislocation cores. During the MD simulations, the computed trajectories of the void surface atoms can be used to derive the evolution of the local plastic displacements according to the dislocation-mediated mechanism of void growth. In addition, the kinetics of void growth can be derived directly by monitoring the void volume evolution. Such representative kinetic results for the void fraction evolution, f ðtÞ, are shown in Fig. 4 for both
Fig. 4. MD simulation results for void fraction evolution, f ðtÞ, in Cu at 100 K under (a) hydrostatic and (b) biaxial tensile strain of 4% (}) and 5% (). In both cases (a) and (b), the solid curves represent fits to the MD results according to Eq. (2) for strains of (1) 4% and (2) 5%.
248
D. Maroudas, M. Rauf Gungor / Computational Materials Science 23 (2002) 242–249
hydrostatic and biaxial straining, Figs. 4(a) and (b), respectively, for two levels of applied strain, 4% and 5%, at 100 K. In both cases, Figs. 4(a) and (b), the solid curves represent least-squares fits to the MD simulation results according to the same phenomenological equation: f ðtÞ ¼ f0 þ Að1 exp ½ ðt=sÞn Þ;
ð2Þ
where f0 ¼ f ðt ¼ 0Þ and A ¼ limt!1 ½f ðtÞ f0 , which represents the void fraction increase when the void volume reaches each steady state. Eq. (2) implies that limt!0 ½f ðtÞ f0 / tn , which is an important scaling relation for the initial stage of void growth. For each mode of straining, the exponent n in Eq. (2) has been found to be invariant over the entire range of applied tensile strain that was examined: from 1% to 8%. Throughout this strain range, n 4=3 for hydrostatic straining and n 3=2 for biaxial straining; for example, n ¼ 1:33 for both curves of Fig. 4(a), while n ¼ 1:51 for both curves of Fig. 4(b). In conjunction with the invariance of this exponent, n, for given mode of straining, it should be mentioned that the void growth mechanism of strain relaxation remains unchanged for each mode over the entire strain range examined. In Eq. (2), the relaxation time, s, depends on the applied strain level: for the results of Fig. 4(a), s ¼ 15:29 and 13.25 ps for strains of 4% and 5%, respectively, while for the results of Fig. 4(b), s ¼ 5:73 and 7.06 ps for strains of 4% and 5%, respectively. In the case of biaxial straining, the shorter relaxation time is understood because of the thin-film finite thickness and presence of the free surfaces. The analysis of MD trajectories to derive void growth kinetics, such as Eq. (2), is particularly promising for the role of MD simulations, as computer experiments, for development of constitutive plasticity relationships that can be incorporated directly into the mesoscopic dynamical framework, Eq. (1). Specifically, the void growth rate, that can be expressed analytically by the time derivative of the right-hand side of Eq. (2), also is expressed as the integral over the surface of the void of the local plastic displacement rate normal to the void surface. In the case of biaxial straining on the (1 1 1) plane of a thin film of thickness h, for example, this relationship is given by
df ðtÞ 6h ¼ dt Vcell
Z
p=3
Rðh; tÞu_ n;pl ðh; tÞ dh;
ð3Þ
0
where Vcell is the simulation supercell volume, R is the local radius of curvature of the cylindrical void surface, h is the corresponding polar angle, and the six-fold symmetric geometry of a cylindrical void with its axis along h1 1 1i has been taken into account to simplify the integral. Most elegantly, instead of directly using the MD simulation results, analytical constitutive theories for the plastic displacement rate can be constructed and serve as parametrization schemes; the corresponding parameters can be quantified by fitting the MD results. Such analytical theories are currently under development and will be reported in a forthcoming publication.
4. Conclusions and future directions A systematic multiscale modeling strategy has been presented for the analysis of electromechanically-induced void dynamics in ductile metallic thin films and void-mediated failure mechanisms in such metallic systems. This multiscale strategy links self-consistent dynamical simulations of void morphological evolution due to surface transport in a deforming solid with atomistic calculations of surface and interfacial properties and MD simulations of plastic deformation mechanisms in the vicinity of the void surface. An important theoretical conclusion of our study is that MD simulations of strain-induced void dynamics can be used to quantify rigorously parametrization schemes provided by mesoscopic or phenomenological theories of void growth and, thus, develop constitutive theories for plastic displacement rates; such constitutive equations can be incorporated directly into the mesoscopic void evolution equations. An additional important conclusion regarding thin-film failure is that, in the most general case, a coupled mode of failure is operative involving the formation of electromigration-induced slits and other morphological features and stress-induced crack-like features at the slit tips. However, the most challenging task ahead is
D. Maroudas, M. Rauf Gungor / Computational Materials Science 23 (2002) 242–249
the identification of possible regimes of electromechanical conditions where stress and electromigration compete to give rise to stable void behavior: a morphologically stable void surface that is either immobile or simply migrates along the film and is driven out of the film. Addressing this challenging problem motivates exciting future research directions for multiscale modeling in ductile metallic thin films, such as the Al and Cu films used in modern interconnect technologies. Such future studies include self-consistent dynamical simulations over a broad range of electromechanical conditions incorporating into the mesoscopic formulation MD-derived plastic response for the evolving void surface. This is a particularly demanding task that requires not only extensive mesoscopic parametric studies and nonlinear analyses but also systematic MD simulation schedules for fundamental mechanistic investigations of stress-induced void dynamics over a wide range of thermomechanical conditions. Multiscale theoretical studies within this framework are currently underway and will be reported in future publications. Acknowledgements This work was supported by the National Science Foundation through a CAREER Award (No. ECS-9501111) to one of the authors (D.M.), the Frontiers of Materials Science Program of the UCSB Materials Research Laboratory and Los Alamos National Laboratory (Award No. STBUC:97-63), and the Oak Ridge Institute for Science and Education through a Faculty Research Fellowship to one of the authors (D.M.). Fruitful discussions with L.J. Gray, S.J. Zhou, G.E. Beltz, and J.S. Langer also are gratefully acknowledged.
249
References [1] P.S. Ho, T. Kwok, Rep. Prog. Phys. 52 (1989) 301. [2] C.V. Thompson, J.R. Lloyd, MRS Bull. 18 (12) (1993) 19, and references therein. [3] C.-K. Hu, J.M.E. Harper, Mater. Chem. Phys. 52 (1998) 5. [4] A.I. Sauter, W.D. Nix, J. Mater. Res. 7 (1992) 1133. [5] A.F. Bower, L.B. Freund, J. Appl. Phys. 74 (1993) 3855. [6] Z. Suo, W. Wang, M. Yang, Appl. Phys. Lett. 64 (1994) 1944. [7] W.Q. Wang, Z. Suo, T.-H. Hao, J. Appl. Phys. 79 (1996) 2394. [8] D. Maroudas, Appl. Phys. Lett. 67 (1995) 798. [9] D. Maroudas, M.N. Enmark, C.M. Leibig, S.T. Pantelides, J. Comp.-Aided Mater. Des. 2 (1995) 231. [10] O. Kraft, E. Arzt, Appl. Phys. Lett. 66 (1995) 2063. [11] O. Kraft, E. Arzt, Acta Mater. 45 (1997) 1599. [12] M. Mahadevan, R.M. Bradley, J. Appl. Phys. 79 (1996) 6840. [13] M. Mahadevan, R.M. Bradley, Phys. Rev. B 59 (1999) 11037. [14] L. Xia, A.F. Bower, Z. Suo, C.F. Shih, J. Mech. Phys. Solids 45 (1997) 1473. [15] M. Schimschak, J. Krug, Phys. Rev. Lett. 80 (1998) 1674. [16] M.R. Gungor, D. Maroudas, Appl. Phys. Lett. 72 (1998) 3452. [17] M.R. Gungor, D. Maroudas, Surf. Sci. 415 (1998) L1055. [18] M.R. Gungor, D. Maroudas, J. Appl. Phys. 85 (1999) 2233. [19] M.R. Gungor, D. Maroudas, Surf. Sci. 461 (2000) L550. [20] see, e.g., W.H. Yang, D.J. Srolovitz, Phys. Rev. Lett. 71 (1993) 1593. [21] W.H. Yang, D.J. Srolovitz, J. Mech. Phys. Solids 42 (1994) 1551. [22] M.R. Gungor, D. Maroudas, L.J. Gray, Appl. Phys. Lett. 73 (1998) 3848. [23] M.R. Gungor, D. Maroudas, Surf. Sci. 432 (1999) L604. [24] M.R. Gungor, D. Maroudas, Int. J. Fract. 109 (2001) 47. [25] S.M. Foiles, M.I. Baskes, M.S. Daw, Phys. Rev. B 33 (1986) 7983. [26] L.J. Gray, D. Maroudas, M.N. Enmark, E.F. D’Azevedo, Eng. Anal. Boundary Elements 23 (1999) 267. [27] L.J. Gray, D. Maroudas, M.N. Enmark, Comp. Mech. 22 (1998) 187. [28] M.R. Gungor, D. Maroudas, S.J. Zhou, Appl. Phys. Lett. 77 (2000) 343. [29] Y.-C. Joo, C.V. Thompson, J. Appl. Phys. 81 (1997) 6062. [30] J.S. Langer, A.E. Lobkovsky, Phys. Rev. E 60 (1999) 6978.