A Comprehensive Microstructure-Aware Electromigration Modelling Framework; Investigation of the Impact of Trench Dimensions in Damascene Copper Interconnects

1. Introduction

Electromigration (EM) refers to the phenomenon where the movement of atoms within a conductor is induced by a high-density electric current. This process can lead to the gradual degradation of the conductor, causing, material displacement, void formation and ultimately, device failure. As electronic devices continue to shrink in size and increase in complexity, the current densities in interconnects drastically increase, intensifying the effects of electromigration. This renders the understanding of EM in modern electronics crucial due to its significant implications for device reliability and longevity. Variability of EM time to failure (TTF) has always been a key challenge for back-end-of-line (BEOL) reliability predictions. Experiments have shown that time to failure has a large range of values for similar interconnects fabricated by the same technology and tested under same working conditions. Many factors contribute to TTF variability including the variation in the actual cross-sectional area of the interconnect [1], the variation in quality of the interfaces between the interconnect’s metal and the surrounding materials and other process related factors such as via misalignments. However, one of the most important intrinsic sources of variation is the microstructure of the interconnect. Microstructure has always been identified as a key contributor to the mechanical, thermal, transport and electrical properties of materials [2], and not surprisingly it has been shown to be strongly implicated in EM [3] [4]. Microstructure has been identified to be the main factor causing shortened lifetime of interconnects with down-scaling of the interconnect cross-sectional area for a given current density [3], making the inclusion of microstructure in the modelling of EM in nanointerconnects a very important aspect of a comprehensive and a predictive model.

Numerous scientific contributions have delved into the intricate relationship between microstructural features and the EM behaviour of nanoscale interconnects. Zhang et. al. [4] conducted EM experiments on two sets of Cu-based interconnects with different average grain sizes, that were realized by manipulating the parameters used during electroplating. A 2-fold increase in lifetime was observed for an average grain size of 198 nm compared to a smaller average grain size of 125 nm. A similar work studying the impact of average grain size was done by Choi et. al. [3], where a wider range of variation in average grain size was achieved by varying the trench width while fixing all other fabrication parameters resulting in a significant enhancement in EM lifetime as the metal line width increases from the minimum to three times the minimum width. However, beyond this limit, further widening of the metal line did not yield any additional improvement in EM lifetime. Oates et al. derived the diffusivity coefficients in different diffusion paths (i.e., grain boundaries and interfaces) for estimation of the drift velocity by conducting EM experiments and texture analysis on damascene interconnects fabricated with different dimensions and technological schemes [5].

Only a handful of studies on numerical modelling of EM have attempted to explicitly simulate microstructure by considering heterogeneous diffusivity domains. In contrast, to reduce model complexity and the entailed computational costs, the widely adopted approach has been homogenization using the effective diffusivity concept [6] [7], where a homogenous atomic/vacancy flux through the entire interconnect’s cross-section is employed to approximate the distributed flux through the grain boundaries and interfaces. One of the few such studies was the seminal work by Gleixner [8] where the atomic flux in each type of diffusivity paths was given in terms of the current density and the angle between the path and the vector of the electric current, giving a stress distribution varying between compression and tensile, and concentrated within the diffusivity paths. Later, the impact of microstructure on void dynamics was introduced by Bower [9] where the diffusivity value along the void surface was considered to vary depending on its interaction with grain boundaries and other interfaces causing a variation in the atomic flux. Subsequently, the normal velocity of the void’s surface was estimated using the normal component of the volumetric atomic flux and the divergence of its tangential component. An attempt to include microstructure in both stress evolution and void dynamics was shown separately by Kteyan and Sukharev [10] using the phase field method. An innovative approach to avail of the benefits of the effective model, being its simplicity and low computational cost, while the impact of microstructure could still be approximated, was implemented in the works of Zahedmanesh [11] and Ceric [12], where the interconnect was divided into segments of polycrystalline and single crystal regions, each having its own effective diffusivity value based on the dominating diffusivity path, that is grain boundaries for polycrystalline regions and the cap interface for the single crystals. In essence, the incorporation of microstructure into simulations of EM is indispensable for navigating the challenges posed by the relentless march towards miniaturization in the ever-evolving landscape of nanoelectronics.

The model presented in this investigation constitutes a major enhancement of our antecedent model expounded comprehensively in [7], where detailed expositions and derivations for the pertinent equations and their implementation can be found. Notably, the antecedent model relied upon a singular effective diffusivity coefficient necessitating calibration for each distinct interconnect dimension and technology. In the current study, this aspect has undergone refinement, wherein the solitary effective diffusivity has been supplanted by a more precise distributed diffusivity scheme. The interconnect cross-section is treated as a heterogeneous diffusive medium with distinct paths within a generated microstructure. This approach provides a more nuanced representation of microstructural influences on diffusive processes capturing the impact of dimensional interconnect scaling and technological variations, seamlessly. A visual comparison between the two approaches is shown in Figure 1. The presented model complements existing literature models due to its comprehensive simulation encompassing all stages of EM, including stress evolution, void nucleation, stress relaxation, and void dynamics, along with their intricate interactions with high diffusivity pathways. Additionally, its realistic microstructure generation module considers interconnect aspect ratios, demonstrating commendable alignment with experimental observations.

The subsequent sections of this study are structured as follows; in Section 2, the implemented methodologies are described, encompassing the algorithms employed for the generation of microstructures and the methodology employed for their integration into EM simulations. Subsequently, in section 3, the comprehensive ability of the presented model to capture all stages of EM failure is demonstrated. The simulation framework was then used to study the impact of trench width (and thereby linewidth) on microstructure and EM lifetime, where the simulation was compared with experimental findings. This validation is then followed by an application of the model beyond the experimental domain to reveal the full scaling and aspect ratio impact on EM reliability. Finaly, in section 4, a conclusive summary of the findings is provided.

2. Methods

As a first step, a realistic grain structure has to be implemented. We propose a 2-dimensional model, where grains are spanning the 3rd dimension (i.e., the linewidth). The synthesis of artificial but realistic microstructures for simulation is implemented by employing a novel methodical process. This process is based on the packing of circles of various sizes in the interconnect domain (Figure 2A), where the radii of the circles follow a prescribed statistical distribution, as shown in Figure 2C. This distribution is obtained from experimental grain size analyses on Cu interconnects of different dimensions [13] [3]. The best-fit analysis yielded a lognormal distribution of grain sizes, with an average size approximating the minimum cross-sectional dimension (i.e., the minimum of line width and height) and a standard deviation of ~ 0.45 × average grain size. Initially, the simulation incorporates dynamic packing of the generated circles, which utilizes a random perturbation force and ensures contact resolution amongst the packed circles to minimize their overlap and the inter-circular empty spaces. An example of the dynamic packing process is provided as a video in the supplementary materials. Our algorithm stands out among other algorithms due to its unique feature of ensuring that all circular objects exist as specified from the beginning, thereby preserving the original experimental size distribution. In contrast to other algorithms, such as the random placement [14], which dynamically adds circles based on available spaces left by earlier placements, the utilized approach guarantees the inclusion of all pre-specified circle sizes. This prevents deviation of the size distribution of the circular objects during the packing process from the pre-scribed distribution. Subsequently, circles are used to create actual grains by employing a weighted Voronoi tessellation, also known as the power diagram [15] (Figure 2B). The circles’ centres are utilized as seeds for the Voronoi diagram, while their radii act as weights, influencing the spatial extent of each Voronoi cell to match as closely as possible with the corresponding circle [15]. This weighted Voronoi tessellation creates a microstructure that faithfully represents the specified grain size probability distribution. An example of the described steps to create a microstructure is illustrated in Figure 2.

The microstructure generation module creates the computational domain devised for the resolution of atomic flux, which is a network comprising high-diffusivity paths, specifically the interconnect interface with the cap layer and the grain boundaries. Concomitantly, flux within the bulk and across metal barrier interfaces is disregarded, as substantiated by markedly lower diffusivity values (~10 order of magnitude lower for bulk) [5] [16]. This geometry domain is then used in our finite element (FE) solver, employed within MATLAB^®, to model the stress evolution with time. The numerical mesh configuration is accordingly tailored to reflect a higher density at triple points, where the divergence of the atomic flux is most pronounced, and a coarser mesh at the midpoints of the paths, as illustrated in Figure 3. For the calculation of the current density, a separate FE solver is used, following the same principles outlined in our preceding publication [7].

The stress evolution with time is determined by the divergence in atomic flux (i.e., depletion and accumulation of atoms) and is given by [7],

$${\displaystyle \frac{{\partial }_{\sigma }}{\partial t}}=B\Omega f\nabla .J_a$$

(1)

where $\sigma$ is the stress, $t$ is time, $B$ is the bulk modulus, $\mathsf{\Omega }$ is the atomic volume, $f$ is the vacancy relaxation volume, and $J_a$ is the atomic flux. The atomic flux is defined as the summation of the flux driven by the electron wind force and the flux driven by the stress gradient [8],

$$J_a={\displaystyle \frac{D_a{F}_e}{\Omega K_{B}T}}+{\displaystyle \frac{D_{a}f\nabla \sigma }{K_{B}T}}$$

(2)

where $D_a$ is the value of atomic diffusivity in the path where the atomic flux is calculated. For instance, $D_a=D_{gb}$ inside grain boundaries and $D_a=D_{cap}$ for atomic flux of the cap interface (see Table 1). $F_e$ is the electron wind force, $K_B$ is the Boltzmann constant, and $T$ is temperature. The electron wind force on each path is given by the amount of momentum exchange between electrons and atoms in the direction tangential to the path [8],

$$F_e=Z^{\mathrm{*}}e\rho j_e\cos \left(\theta \right)$$

(3)

where $Z^{*}$ is the effective charge, $e$ is the electronic charge, $\rho$ is the resistivity, $j_e$ is the current density, and $\theta$ is the angle that the diffusivity path of interest makes with the direction of the electric current flow. In this computational framework, the determination of the atomic flux within each path is initiated by employing equations (2) and (3), with a representative calculation illustrated in Figure 4 for the case of electrons flowing from left to right. At the interconnect boundaries, such as at the end terminals isolated by diffusion barriers, a blocked boundary condition that ensures a zero normal component of the total flux was applied.

The atomic flux takes place along the fast diffusion paths, being higher when the path aligns more closely with the electron flow direction, in contrast to paths perpendicular to it. Subsequently, the flux divergence is computed, elucidating regions where atoms accumulate or deplete, notably at triple points. An illustrative calculation of flux divergence is presented in Figure 5. Finally, the evolution of stress is evaluated for each discrete time step, employing the formulation provided in equation (1). The subsequent phases of the simulation, encompassing void nucleation and void dynamics adhere to the methodologies expounded in our prior publication [7].

For void nucleation, it is imperative to surpass a critical tensile stress threshold, ${\sigma }_{crit}$, initiating void growth. Triple points within the microstructure are conducive to nucleation due to their role as focal points for stress evolution and flux divergence. Moreover, the lower activation energy of atoms at triple points reduces the critical stress required for voiding [17]. The stress at the void surface is introduced as a boundary condition post-nucleation. Utilizing the principle of continuity of chemical potential, the stress on the void’s surface, ${\sigma }_s$, is related to surface energy, $\gamma$, and curvature, $k$, by equation (4) [7].

$${\sigma }_s=\gamma k$$

(4)

Assuming a circular initial void geometry, the initial void radius, $r_0$, is estimated as outlined in equation (5) [18]. This equation was deduced from minimizing the free energy of a stressed continuum containing a void by variating its size.

$$r_0={\displaystyle \frac{2\gamma }{{\sigma }_{crit}}}$$

(5)

The nucleated void generates a new high diffusivity path represented by its free surface. The normal velocity of the evolving void surface under mass transport, $V_{\perp }$, is determined by equation (6) using the calculated flux on the void surface [9].

$$V_{\perp }=-{\delta }_s\nabla .j_{{vol}_{//}}+j_{{vol}_{\perp }}$$

(6)

where $j_{{vol}_{//}}$ is the volumetric flux tangential to the void surface, $j_{{vol}_{\perp }}$is the volumetric flux normal to the void surface, and ${\delta }_s$ is the void’s surface phase thickness.

3. Results and Discussion

The simulation results on stress evolution, void nucleation, and void dynamics for a 100 $\mathsf{\mu }\mathrm{m}$long dual damascene Cu interconnect with a linewidth of 40 nm and an aspect ratio (AR) of 2 with a conventional SiCN cap layer, are shown in Figure 6 zoomed on the region of the cathode end for clarity. See Table 1 for the physical parameters and working conditions employed in the simulations.

Because of the substantial difference in stress ranges during the initial phase of stress evolution and subsequent stages of void nucleation and dynamics, distinct color bars are employed to ensure clarity and enhance contrast in stress visualization across all phases. A comprehensive animation illustrating the complete simulation is provided in the supplementary materials, while selected snapshots of interest are highlighted in Figure 6. At the onset (time 0), denoted as the initiation of current flow through the interconnect, a uniform initial stress of 185 MPa is assumed, as estimated in our dedicated thermal-mechanical study [25]. This initial tensile stress originates during the cooling phase from the stress-free temperature down to the working condition temperature, attributed to the differing thermal expansion coefficients of copper (Cu) and the surrounding materials [25]. After current turn-on, momentum exchange between electron wind and the Cu atoms induces an atomic flux through the high-diffusivity paths. Depending on the directions and magnitudes of these fluxes, flux divergence results in either atomic depletion, augmenting the initial tensile stress, or atomic accumulation, shifting the stress towards compression. Large single crystals serve as flux divergence segments, with atomic flux locally occurring solely through the cap interface. In contrast, polycrystalline segments exhibit higher flux due to the presence of both cap interface and grain boundaries as high-diffusivity paths. Subsequently, an increasing number of atoms diffuse toward the anode end, leading to the highest tensile stress near the cathode end, see Figure 6. This stress evolution continues until reaching the minimum critical stress for voiding, assumed to occur at a triple point on the cap interface [26], see Figure 6F. At this juncture, a void nucleates, and tensile stress locally relaxes around the nucleated void [27]. The void undergoes growth primarily under the influence of atomic flux driven by the stress gradient, with regions intersecting high-diffusivity paths expanding more rapidly. As the void continues to grow, the stored mechanical energy is released, leading to stress relaxation and a subsequent weakening of the initially strong gradient, as depicted in Figure 7A, where the onset (time 0) is denoted in this figure as the time of void nucleation. Simultaneously, the reduction of the interconnect’s cross-sectional area due to void growth results in an increase of interconnect resistance (see Figure 7B), generating about 3x increase in current density around the void surface and propelling accelerated void migration toward the cathode end.

Using the developed modelling framework, a comprehensive case study investigating the impact of trench width (and thereby the damascene interconnect’s linewidth) on electromigration was investigated and subsequently compared with experimental data [3] obtained from interconnects with an identical height of 80 nm. Specifically, the experimental trials in [3] involved interconnects with linewidths of 32 nm, 64 nm, 96 nm, and 156 nm. In the simulations, linewidth variations were systematically introduced, ranging from 20 nm to 150 nm. For each simulation data point, a short segment of the generated microstructure at the corresponding scale is exemplified in conjunction with the average grain size (d) in Figure 8.

It is evident that for interconnects with AR<1, the microstructure is predominantly bamboo, and transitions towards increased polycrystallinity and smaller average grain size at narrower linewidths. Thereby, the increased grain boundary density increases the atomic flux and accelerates stress evolution, consequently leading to a reduced time to nucleation. The elongation in time required for nucleation grows proportionally with the expansion of line width, exhibiting a linear trend until reaching an aspect ratio of 0.8, below which it stabilizes. This saturation emerges because, within this range, the constant line height acts as the limiting factor for the average grain size, leading to a saturation point in total diffusivity value and consequently the time needed for nucleation. In Figure 8, the simulation results are compared to experimental findings in [3] demonstrating a correlation that proves the validity of the presented model.

Expanding the scope of this investigation to encompass various line heights provides a comprehensive understanding of the scaling effects on EM lifetime, as depicted in Figure 9. Utilizing simulated data points (see Figure 9 left hand side plot) and using 2D interpolation, the time to nucleation is estimated for each combination of line height and width, extending up to and including 250 nm (see Figure 9 right hand side plot).

The consistent trend of pronounced reliance on linewidth for high aspect ratio lines, followed by a plateauing effect, is evident across all line heights. However, due to the variation in the value of line width at which aspect ratio 1 is reached for each line height, the point of saturation varies accordingly, where doubling the line height for a fixed large line width (i.e., lying in the saturation region) was found to give about 2.5x increase in the value of nucleation time saturation. However, for a fixed small line width (i.e., lying in the aspect ratio 1 and above region) a reduction in time to nucleation is noted as line height increases, even if the minimum cross-sectional dimension dictating the average grain size is the fixed line width. This influence gradually reaches a plateau with further increase in line height. This phenomenon can be attributed to the probability of encountering a single crystal region or a polycrystalline region along the length of the line, quantified by the “p” factor, which is the probability of encountering a polycrystalline region [28] and defined by,

$$P={\displaystyle \frac{L_p}{L_p+L_s}}$$

(7)

where $L_p$ and $L_s$ are the mean length of the polycrystalline and single crystal regions, respectively. In a single crystal region, the atomic flux relies solely on the cap interface as the primary high diffusivity pathway. Consequently, even with a fixed average grain size dictated by the minimum cross-sectional dimension, elevating the line height diminishes the likelihood of single crystal regions to persist (indicated by a higher “p” factor), until such regions become absent for high aspect ratios (resulting in the “p” factor reaching a maximum of 1). Consequently, subsequent increases in line height exert minimal influence on the time to nucleation. The “p” factor as a function of aspect ratio (AR) is shown in Figure 10, assuming the quality of interfaces is the same across different scales. Results from our simulations in Figure 10, are in perfect agreement with reported experimental studies in [11] for the “p” factor measured at aspect ratios 2 and 3.

In addition to previously discussed cases of varying line height while keeping width constant and vice versa, represented by vertical and horizontal lines on the full scaling colormap (Figure 9, right), there are two additional lines of interest. These include iso-volumetric lines, indicating a fixed cross-sectional area with different aspect ratios (Figure 11, left), and fixed aspect ratio scaling lines, indicating different cross-sectional areas (Figure 11, right). Results indicate that an aspect ratio of 1 maximizes EM lifetime for each cross-sectional area, with an estimated 2-fold increase when moving from an aspect ratio of 2 to 1. This can be explained by the fact that at an aspect ratio of 1, both line width and height are the minimum dimension, allowing microstructure grains to fill the entire interconnect trench cross-section. This fosters a bamboo-like structure, reduces the probability of polycrystalline regions, creates a low diffusivity medium for atomic flux, and extends EM lifetime.

Examining the escalation of resistance post-void nucleation across various dimensions provides insights into the influence of scaling on void growth, as depicted in Figure 11. The figure illustrates three distinct combinations of line height and width. Notably, time 0 is defined as the moment of void nucleation for each case, acknowledging the varying times to nucleation for each case as highlighted earlier.

Comparing the scenario where both line height and width are 80 nm to the case where the line height is 80 nm and the width is 40 nm, illustrates the consequence of altering line width while maintaining a constant line height. Doubling the line width results in approximately a 200-fold increase in the time required to achieve the same rise in resistance. This disparity arises from the heightened presence of grain boundaries interacting with the void surface, attributable to a smaller average grain size in the case of a smaller width. Consequently, a larger flux of atoms migrates away from the void surface, accelerating its growth. Likewise, comparing the scenario where both line height and width are 80 nm to the case where the line height is 160 nm, and the width is 80 nm elucidates the effect of varying line height while keeping line width constant. Doubling the line height leads to approximately a 500-fold increase in the time required to attain the same rise in resistance. This discrepancy is attributed to the larger cross-sectional area, necessitating a greater distance for the void surface to traverse to generate an equivalent resistive impact in the case of a larger height. Finally, comparing the scenario where the line height is 80 nm, and the width is 40 nm with the case where the line height is 160 nm, and the width is 80 nm shows the impact of scaling both dimensions while maintaining a constant aspect ratio. This comparison highlights the most pronounced variation in the time needed for the void to generate the same resistive impact, attributable to the combined influence of a larger average grain size and a larger cross-sectional area.

4. Conclusion

In conclusion, our study presents a comprehensive simulation framework for assessing electromigration failure in nanointerconnects, wherein the nuanced impacts of microstructure, and linewidth scaling, are systematically examined. Our simulation framework captures the entire spectrum of EM stages, including stress evolution, void nucleation, and void dynamics. The modelling of stress evolution unveils the intricate interplay of atomic flux, stress gradients, and microstructural features, shedding light on the dynamics leading to void nucleation. The subsequent void growth phase is characterized by complex interactions, including the release of stored mechanical energy and changes in stress gradients, influencing void growth rates and interconnect resistance. According to our model’s predictions, under the assumption of consistent interface quality across various scales and for the same cross-sectional area, the proportion of polycrystalline segments within the interconnect substantially rises from approximately 10% to around 90% of the total length causing a reduction in lifetime as the aspect ratio escalates from 1 to 2. Notably, for a given cross-sectional area, an aspect ratio of 1 optimizes EM lifetime by minimizing grain boundary density and fostering a bamboo-like structure.