Microstructure and Granularity Effects in Electromigration

The persistent advancements made in the scaling and vertical implementation of front-end-of-line transistors has reached a point where the back-end-of-line metallization has become the bottleneck to circuit speed and performance. The continued scaling of metal interconnects at the nanometer scale has shown that their behavior is far from that expected from bulk films, primarily due to the increased influence that the microstructure and granularity plays on the conductive and electromigration behavior. The impact of microstructure is noted by a sharp shift in the changing crystal orientation at the grain boundaries and the roughness introduced at the interfaces between metal films and the surrounding dielectric or insulating layers. These locations are primary scattering centers of conducting electrons, impacting a film’s electrical conductivity, but they also impact the diffusion of atoms through the film during electromigration. Therefore, being able to fully understand and model the impact of the microstructure on these phenomena has become increasingly important and challenging, because the boundaries and interfaces must be treated independently from the grain bulk, where continuum simulations become insufficient. In light of this, recent advances in modeling electromigration in nanometer sized copper interconnects are described, which use spatial material parameters to identify the locations of the grain boundaries and material interfaces. This method reproduces the vacancy concentration in thin copper interconnects, while allowing to study the impact of grain size and microstructure on copper interconnect lifetimes.


I. INTRODUCTION
The continued trend in transistor scaling along Moore's Law [1] has, for the most part, been accompanied by simultaneous miniaturization of the interconnect lines. The miniaturization of metal lines down to the size of several nanometers results in an increased impact of the material interfaces (MIs) and grain boundaries (GBs) on the conductivity and reliability of the thin film, which was mainly composed of copper. The impact of the microstructure includes the impact of line and via sidewall roughness, the intersection of porous low-κ voids with the sidewall, copper (Cu) surface and copper/barrier interface roughness, and the presence of GBs [2]. Understanding and mitigating the impact of granularity in interconnects is essential for scaling to continue. Alternatively, new materials will have to be used, whose lifetime and resistivity to electromigration will also need a closer investigation. Of particular interest in this regard is cobalt, which is applied for use in combination with copper for M0 and M1 metallization [3], [4], shown in Fig. 1.Another interesting structure being investigated to replace the copper interconnect is a carbon nanotube [5], [6]. However, integrating a completely novel structure and material is complex and the accepted reality is that we will live with copper for the foreseeable future, at least down to the 3 nm node [7]. Ultimately, copper is expected to continue to be used for the next several technology nodes [4], [8], possibly in combination with cobalt at the M0 and M1 layers and on its own in higher metal stacks. The layers critical for electromigration (EM) are those where the lines have the smallest cross-sectional area, or pitch, resulting in the highest current density. Most often, these are lines closest to the front-end-of-line (FEOL) transistors, nearest to M0 from Fig. 1 which shows the back-end-of-line (BEOL) metallization for Intel's 10 nm technology node [3]. Here, the first five layers (M0-M5) require metal lines with a pitch below 45 nm [3], making them susceptible to EM failure and granularity effects. Therefore, it must be ensured that EM is properly modeled in order to be able to accurately estimate interconnect lifetimes. In this review a framework developed to accelerate EM simulations of grained copper interconnect lines is presented. It should be noted that this framework is not limited to copper only, but can easily be extended to deal with any number of materials which are prone to EM failure and require the inclusion of granularity effects for proper understanding of their conductive and reliability behavior [9].
With increased scaling and reduction in the interconnect half-pitch with every new generation, it was noted that the lifetime of copper interconnects has approximately halved for every new technology node, even when the current density is kept the same [10]. As the thickness is reduced, the average crystal grain size decreases almost linearly [11], as shown in Fig. 2 [9]. The decreased grain size, combined with the overall reduction in metal thickness, means that GBs and MIs play an increasingly important role in determining the film behavior. The influence of these properties on electron scattering, and thus on conductivity, has been explored for many decades, starting with Fuchs [12], Mayadas and Schatzkes [13], and Sondheimer [14]. In addition to the changes in its conductive behavior, the reliability of nanometer sized copper films is significantly influenced by their microstructure. Degradation due to EM is the primary form of failure in metal films, which occurs due to the transport and accumulation of vacancies which ultimately nucleate to form a void. Under a high current density, this void grows to increase the line resistance and finally to FIGURE 2. Relationship between film thickness and resulting average grain diameter. The symbols are measurements from [11], while the lines represent the best-fit linear regression. As the film thickness is decreased, the average grain size decreases, meaning that more grain boundaries are present.
cause an open circuit failure [15]. Alternatively, the stress induced by the accumulation of vacancies can sometimes be significant enough to cause cracking and a failure by itself [15], [16].

II. ELECTROMIGRATION PHYSICS
The physics of EM phenomena is described in detail in, e.g., the work of Ceric and Selberherr [17]. There are two main driving forces for electromigration, the direct force F direct , initiated through the direct action of the external field on the charge of the migrating ion, and the wind force F wind , arising due to the scattering of the conduction electrons by impurities or point defects. The total force F is given by the sum of the two forces Each force is described by its defect valence Z direct or Z wind , respectively, according to The effective valence Z * is commonly used to describe the total force and is the sum of the direct Z direct and wind Z wind valences (Z * = Z direct + Z wind ), ultimately allowing to represent the total force using where e is the elementary charge and E is the electric field. In free-electron-like metals (e.g., Cu) F wind is dominant.
To experimentally determine EM failure on a new technology, the failure behavior of all materials making up the interconnect must be characterized. This is commonly performed under accelerated conditions (i.e., high temperature and high current) to find the mean time to failure (MTTF) which is then extracted to operating conditions using a log-normal plot following Black's equation [18], given by where n and E t are the density exponent and activation energy, respectively, A is a constant determined by the material properties and geometry of the interconnect, j is the current density, T is temperature, and k B is Boltzmann's constant. It should be noted that simulations, just like measurements, are carried out using the same accelerated conditions in order to as closely as possible match the experimental conditions.

A. CONDUCTIVITY
The current density which builds up in metal films is one of the main drivers of vacancy accumulation and EM failure. Therefore, it is important to understand how the conductivity changes with increased impact of microstructure. To assess the conductivity, three main components must be included: 1) Intrinsic resistivity of bulk copper, limited only by the electron mean free path (MFP).
2) The decrease in conductivity due to the electron scattering at copper's surfaces, which includes material interfaces. The surface roughness also plays a role in determining this impact, which is why different adjoining materials (e.g., dielectrics or insulators) will have a different resulting conductive behavior of the copper line.
3) The decrease in conductivity due to the impact of grain boundaries. The conducting electrons will scatter when reaching a grain boundary, limiting their MFP. In order to introduce the influence of GBs and MIs on the copper conductivity in a continuum way, several models have been proposed [12], [13], [14], [19], [20]. The effects of the granular microstructure, surface scattering on GBs and MIs, and cross-sectional area of a copper interconnect on its resistivity ρ f is commonly modeled by applying a continuum equation derived by Clarke et al. [19], based on the works of Mayadas and Shatzkes [13] where ρ i is the bulk resistivity, λ is the MFP, w is the metal width, p is the probability of electron scattering from a MI, D is the average grain diameter, and R is the probability of electron scattering from a GB. The added temperature influence on the final resistivity ρ can be calculated using where T is the temperature, T ref is the reference (room) temperature, and α e = 0.0043K −1 is the temperaturedependent factor for copper resistivity, referred to as the temperature coefficient of resistance (TCR). Several studies suggest that TCR varies according to the microstructure in several metals, when looking at the difference between nanometer sized and micrometer sized grains [21], [22]. However, this was shown not to be the case for nanometer sized copper grains where the size varied between 20 nm and 120 nm [23], a range of interest for the study presented in this manuscript, and in copper lines from different data sets [24]. Most simulations dealing with the conductivity and reliability of copper interconnects use equation (5). This approach provides a new bulk value for the microstructure-dependent resistivity, while the effects of an individual grain boundary cannot be analyzed with this model. For this, we need to make sure the entire line, with its microstructure, is represented in an appropriate simulation environment. This can be done by observing that the proximity to GBs and MIs can be treated as a parameter, which influences conductivity. Therefore, finding the local conductivity based on the distance of each point inside the grain to a GB and MI can help create a spatial representation of the conductivity inside a copper line. Once the distance to the boundaries d b from every point inside the individual metal grains is known, the local resistivities ρ l and conductivities σ l , which depend on the nearest GB and MI, are derived from the intrinsic resistivity ρ i using Conductivity is one of the primary properties which influences the electromigration behavior of copper. A high current density and a high electron wind can lead to the diffusion of metal atoms in the direction of electron motion. The diffusion is governed by the atom diffusivity property of the material, which also varies depending on whether the atom is located in the grain, on the grain boundary, or along the interface between the metal and adjacent material. Because atoms are more strongly bound inside the grain lattice than at the grain boundaries, their migration is more likely to take place along the boundary, meaning that their diffusivity there is increased [9]. It should be noted that self-heating leading to thermo-migration (TM) has an additional effect in vacancy dynamics, which is also included in the presented framework. However, TM is commonly ignored in EM measurements due to the accelerated conditions quickly providing a uniform temperature and minimizing temperature gradients.

B. VACANCY DYNAMICS
The main driver of electromigration is the accumulation of vacancies which then form a void. The diffusion of vacancies D v through a material with pre-exponential diffusivity D v0 is determined by where E a is the activation energy, is the atomic volume, and σ is the hydrostatic stress. The vacancy diffusion determines its flux J v using where C v is the vacancy concentration and F( j, T, σ ) is a function which depends on the current density j, temperature, and hydrostatic stress, with where ρ is the resistivity (EM component), Q * is the heat of transport (TM component), and f is the vacancy relaxation ration (stress-migration). The subsequent accumulation and depletion of vacancies is found using the continuity equation where G is a surface function which describes vacancy generation and annihilation. Furthermore, from (10) it is clear that the resistivity ρ plays a significant role in determining the vacancy flux and thereby in the overall EM behavior. The discussion thus far describes the EM process in a bulk material, which can be modeled assuming a continuum in the material properties. However, granularity can modify this view significantly, as discussed in the next section.

C. EFFECT OF GRANULARITY ON VACANCY DYNAMICS
In addition to the conductivity, the diffusion of vacancies D v is different between atoms located in the GB, MI, or in the grain bulk. From (8) it has been shown that both D v0 and E a depend on the atom's location in the granular structure of copper according to Table 1. Of note is that the atomic diffusivity in MIs is three orders of magnitude larger than the bulk value, which explains why MIs play such an increasing role for electromigration in nanometer sized interconnects. Therefore it is essential that these parameters are properly treated in any EM simulations. Another aspect of EM, which is ignored in continuum models, is that the generation and annihilation of vacancies G, according to (11) only takes place inside the GBs and MIs and not in the grain bulk. The equation which governs this process is given by where C v,T and C v,eq are the trapped and equilibrium vacancy concentrations, respectively, τ is the relaxation time, and ω R and ω T are the vacancy release and trapping rates, respectively. In (12) χ is a step function which is assigned a value of 1 inside the GB and MI, and 0 otherwise. Therefore, a total of four spatial parameters is used to sufficiently include granularity in EM models, those being σ l , D v0 , E a , and χ . A simulation framework designed to implement this model is given here. Solving equations (8) to (12) gives the time dependent change of the vacancy concentration inside the copper film. The accumulation of vacancies at one end of the wire and their bunching on the other end results in an increase in tensile and compressive stresses, respectively. Once a critical stress level is reached, the material can no longer conduct sufficient current for the required application, resulting in failure.

III. ELECTROMIGRATION SIMULATION FRAMEWORK
The simulation framework relies on three main components, namely Voronoi tessellation to generate the grained interconnect line, the assignment of the relevant granularitydependent parameters discussed in the previous section (σ l , D v0 , E a , and χ ), and the solution of the EM equations to find the vacancy accumulation and resulting EM-induced stress, as visualized in Fig. 3. These three components are addressed in further detail in this section.

A. TESSELLATION
The stochastic polycrystalline copper line is generated using a Voronoi tessellation. Assuming spherical grains and knowing the average grain size, the total number of grains which fit into the volume is found. For each grain a seed point is placed at a random location inside the metal line, which then grows isotropically, until the entire volume is filled. When grains hit each other, they merge to form a GB. The Neper tessellation tool [25] is used in order to create the required tesselated structures, allowing for the generation of a Vonoroi tessellation with ideal copper orientations of (1 1 2)[1 11] and (1 1 2)[11 1] [26].
In order to show the key features of the tessellation tool used in this study, the above mentioned technique was applied to copper lines with different average grain sizes. The two-dimensional (2D) lines have a length of 1000 nm and height of 20 nm with an average grain size of 20 nm and 40 nm, shown in the top and bottom sections of Fig. 4, respectively. In the wire with smaller grains, we note a very granular structure. However, when the grain diameter is larger than the dimensions of the metal wire, a bamboolike structure is formed, as depicted in the bottom of Fig. 4. This is consistent with many studies which show that as a wire becomes narrower the grains begin to be more bamboolike and less granular [20], [27], [28]. This means that the GBs are primarily near-perpendicular to the direction of the applied electric field and the direction of the atom diffusion. Therefore, they simultaneously act as fast diffusivity pathways and diffusion barriers, depending on the grain boundary orientation.

B. SPATIAL PARAMETER ASSIGNMENT
The assignment of spatial parameters (σ l , D v0 , E a , and χ ) on a Cartesian grid with spacing d g ensures that the GB and MI locations are explicitly defined and that the EM framework properly treats the granular nature of the interconnect line. Linear interpolation is used in the EM model in order to populate the entire material domain between the defined points. This proceeds according to the flow chart in Fig. 5. A boundary thickness of 1 nm was assumed here as was found to be appropriate from previous studies [29].
A 2D test geometry with dimensions 20 nm × 2000 nm and a grain diameter of 25 nm was used to test the given framework. Due to the 2D nature of the test geometry, the simulation is effectively performed on a Cu sheet and the line width is ignored. The results of the spatial parameter assignment for the diffusivity (D v = D v0 e E a /k B T ) on one section of the structure are shown in Fig. 6. The impact of the GBs and MIs is evident. Noteworthy is the dependence of the diffusivity on the angle between the GB and the current flow, or the direction of the electron wind. When the GB is perpendicular to the flow, the diffusivity is almost zero and the grain boundary acts as a vacancy blocking site. On the other hand, GBs which are parallel to the current flow have  a diffusivity which is higher than that of the bulk material, speeding up the vacancy transport.

C. ELECTROMIGRATION MODEL
The equations which govern EM physics are given in Section II. To model EM requires the solution of three physical phenomena simultaneously, including the electro-thermal problem (current density, self-heating, and temperature), the vacancy dynamics problem, and the solid mechanics problem (induced strain and stress) [30]. Ultimately, finding the induced stress is desired in order to ascertain whether the critical stress is reached, which results in failure or the formation of a void. The flowchart of the electromigration model is given in Fig. 7 [30].
Solving the electro-thermal problem allows to identify the temperature distribution and current density in the interconnect. Joule heating must be considered, as this can lead to higher thermal gradients in the interconnect and an increased proclivity to the thermo-migration component of vacancy diffusion. The vacancy dynamics problem is solved by finding a solution to equations (8) to (12), as given in Section II. Finally, in order to calculate the vacancy-induced strain, the solid mechanics problem must be solved. The change in volume caused by the migration and formation of vacancies is represented by [31] ∂ε where ε is the trace of the strain tensor, while ε m and ε f represent the strain induced due to the migration and formation of vacancies. The second term is multiplied by χ because vacancy formation takes place only at the grain boundaries and material interfaces. The induced strain can then be derived to The above equation (14) connects the vacancy transport ( J v ) and mechanics (ε) models. Given that the strain in a dual-damascene interconnect is anistotropic [17], the induced strain is modeled by applying one third of the strain in (14) in each Cartesian direction. From (9) and (14), it can be observed that the vacancy flux depends on the stress, which in-tern depends on the vacancy flux. In order to solve these mutually-dependent sets of equations, time discretization is required and the time steps must be small enough to ensure that the induced error is minimal. The entire flow sequence shown in Fig. 7 is solved at every time step and segregated solvers are used to calculate each of the electro-thermal, vacancy dynamics, and solid mechanics problems. Newton's method is used to obtain a solution for the entire set of equations at each step.

IV. RESULTS AND DISCUSSION
For the results obtained in this section, a 2D copper line is simulated, with a height of 20 nm. The applied current density if 1MA/cm 2 at an ambient temperature of 300 • C. Furthermore, the results are analyzed primarily for the electromigration component during the relatively early stages of vacancy dynamics up to the point where the electromigration induced vacancy transport balances out the stress-induced transport. After this point, the stress-induced component takes over, further increasing the stress until eventual failure. The simulated times are long enough to ensure that stressmigration is the dominant vacancy dynamics effect. This allows to fully encompass the EM and TM phenomena in the result. Thereafter, the stress-time curve can be extrapolated without the need to solve the complex EM model.

A. IMPACT OF GRAIN SIZE
We performed several simulations on a 1000 nm long copper line while varying the average grain diameter (D g ) from 15 nm to 40 nm. The impact of D g on the maximum vacancy concentration (C v /C v0 − 1) and the maximum induced stress (σ ) are given in Fig. 8 and Fig. 9, respectively. The highest vacancy concentration and EM-induced stress are observed at the end of the Cu line, downwind the electron motion. This is due to the assumption of zero vacancy diffusivity there, causing the highest accumulation of vacancies, independent of the grain size. The first observation from the figures is that C v and σ increase with decreasing D g . This is not surprising, since smaller grains mean that there are less columnar GBs, which act as diffusion barriers, and more GBs in parallel to the electron wind, accelerating the vacancy transport. This was shown previously in Fig. 4.
In order to analyze whether a continuum model could be devised which replicates the behavior of the framework presented here, we analyzed the average values for the spatial parameters of interest, namely conductivity, diffusivity (D v = D v0 e −E v /k B T ), and χ , the ratio of the interconnect volume where vacancy generation and annihilation can take place. The calculated values are given in Table 2. We note  that the average diffusivity does not vary much with increasing grain size, suggesting that a continuum model, which relies on bulk parameter representations of the copper film, might not be easily attainable. Using the values from Table 2 directly resulted in an underestimation of the maximum vacancy accumulation and induced stress, because the continuum model is not able to properly represent the gradients which occur in the film due to the presence of complex interfaces and boundaries. Regardless how we varied the parameters, it was not possible to reproduce both the vacancy concentration and induced stress graphs which were obtained with the microstructure simulations.

B. LOCAL STRESS
Another aspect which cannot be properly treated with a continuum model is the representation of local stresses. For example, triple points (where a GB and MI meet) are known to cause a slight increase in the induced stress, compared to its surrounding. We performed a sample simulation on a 2000 nm × 20 nm copper line with an average grain diameter of 25 nm and the resulting EM-induced stress after 200 s is given in Fig. 10. Here, the influence of the GBs and MIs on EM is evident. The framework is able to reproduce the stress generation at triple points σ TP , shown in the circled regions in Fig. 10, including at (x, y) = (287 nm, 0.5 nm), (x, y) = (287 nm, 19.5 nm), (x, y) = (337 nm, 0.5 nm), and (x, y) = (337 nm, 0.5 nm). With the presented simulation framework, this stress can be accurately modeled, even with very coarse meshes [9], [30]. In fact, when the mesh for the simulations was varied from 0.4 nm to 2 nm (25× speedup in 2D), the variation was under 5% (the parameter grid d g from Fig. 5 was set to 0.1 nm).

C. SIMULATION TIME
The proposed framework allows for a very quick and efficient estimation of EM phenomena while taking the film's granularity into consideration. In Fig. 11 the simulation time is plotted against the chosen grid spacing, when the grid during the electromigration simulation is varied. The spatial parameter grid d g is either set to 0.1 nm ( ) or is varied together with the electromigration grid ( ). A drastic reduction in the simulation time can be achieved by increasing the coarseness of the mesh with relatively little loss of accuracy. For the entire simulation range shown in Fig. 11, the variation in the stress at triple points varied by less than 8%. Therefore, when the goal is to model local stresses, even very coarse meshes will suffice, allowing for simulation times in the order of a few minutes.

V. CONCLUSION
Continuum EM models frequently underestimate the time at which EM effects initiate, due to their inability to properly take into account the granularity of nanometer sized interconnects. The effects of granularity (GBs and MIs) have been known to exacerbate the electromigration phenomena. Therefore, it is essential that they are properly treated. Here, a sophisticated modeling framework is described, which considers granularity by applying spatial material parameters (σ l , D v0 , E a , and χ ) in EM simulation to identify the locations of the GBs and MIs. The framework allows to efficiently model the impact of the average grain size on the resulting EM behavior as well as to study induced local stresses, such as those at triple points, even when very coarse meshes are applied to accelerate the simulation time.