In this work, we utilize a phase-field model to investigate electromigration-mediated defects in non-columnar polycrystalline interconnects. We find that the misalignment of the grain boundary with respect to an externally applied electric field governs the non-linear kinetics of electromigration-induced slit evolution. We explore the mechanisms by which electromigration-induced defects propagate in interconnects comprising equiaxed and randomly distributed grains. We deduce that when atomic mobility in grain boundaries (Mgb) is two orders of magnitude larger than along the surface (Ms), the defect manifests as grain boundary slits, while a smaller Mgb/Ms promotes surface drift. By the aid of an extensive parametric study, the presence of a mixed mode at intermittent values of Mgb/Ms is established. Our simulations of slit formation in a network of randomly distributed grains validate our hypothesis that grain boundary alignment and the grain size distribution determine failure rates. Finally, we found that the failure rates in 3D are several times faster than in 2D, which indicates the importance of accounting the physics of three-dimensional capillarity in the present modeling approach.

Electromigration (EM) is the transport of materials caused due to the momentum transfer between conducting electrons and the diffusing metal atoms.1 Technological downscaling has resulted in increased current densities in micro- and nanoelectronic components exposing them to enhanced electromigration, which often leads to failure.1–4 Challenges in predicting and controlling the proliferation of EM defects that are borne out of technology scaling have been highlighted in the ITRS roadmap of 2011.5 Although electromigration can never be completely mitigated, preventive strategies can help in prolonging the life of electronic circuits. Devising preventive strategies and their efficient implementation warrant a detailed understanding of the failure mechanisms, which are directly linked to interconnect microstructures.

EM failure in interconnects manifests in the form of voids and hillocks’ formation, grain-boundary (GB) slit propagation, and metallization, all of which significantly hamper the microelectronics’ reliability limiting our capabilities to downscale electronic circuits and chips. Complex interactions of electrical, thermal, and stress fields obscure our understanding of the diffusion mechanisms by which EM-mediated failures in interconnects and solder materials occur while making it impossible to predict the onset and propagation of defects. Therefore, a comprehensive understanding of diffusional mechanisms in electronic materials is required for mitigating the deleterious effects of electromigration on interconnect reliability.

In interconnects that comprise bamboo-like grain structures, GBs are aligned nearly perpendicular to current flow as opposed to a typical polycrystalline microstructure consisting of a random distribution of grains. Bamboo microstructures are more resistant to slit formation due to a 90° misalignment of grain boundaries with respect to current flow. Numerous studies have investigated the defect pathways anticipated in bamboo microstructures, specifically in the context of interconnect reliability.6–8 On the contrary, studies reporting EM-mediated defects in a polycrystalline interconnect consisting of a random distribution of grains have been rare and limited to microstructures comprising columnar grains.3,4,9

Previous finite element studies have shown that the interface curvature influences the local distribution of an electric field.10–12 An increased current density facilitates electromigration more readily while also increasing Joule heating locally. The increase in the temperature also increases the diffusion rates, which accelerates electromigration, while the temperature gradients induced by localized Joule heating can cause thermomigration.13 Combined action can result in void and slit formation. Once such defects form, current flows can occur around them, diverting the flow of electrons and creating areas of increased current density. An understanding of how these events influence void and slit propagation is important from a fundamental standpoint.

In an inert environment, atomic diffusion along the surface of the interconnect is typically much larger than along the grain boundaries.14,15 However, as interconnects age, they begin to develop surface defects, which cannot only hamper reliability16,17 but also enhance the diffusion rates in the grain boundary relative to the surface.18 More comparable atomic mobility coefficients along the surface and the grain boundaries can result in distinct defect evolution modes. Moreover, diffusion along grain boundaries and its relative misalignment with the electric current direction can significantly impact the morphological evolution of electromigration-mediated defects. Therefore, the role of grain boundary misalignment and the competing diffusion pathways on the evolution of electromigration-mediated defects need to be explored in order to devise techniques for mitigating the interconnect failure. Naturally, this begs the question: how can one predict the morphological evolution of EM-induced defects in polycrystalline interconnects that may comprise a network of grain boundaries, which are misaligned with respect to the direction of electric current flow? The objective of the present work is to explore this question through the influence of electric current and diffusion pathways on the electromigration-induced grain boundary grooving in polycrystalline interconnects. Although isolated attempts for analyzing the evolution of EM-mediated growth of intermetallic compounds (IMCs),19,20 slits,3,4,21,22 hillocks,23 and voids24 have been made in the past via computational or mathematical techniques, the derived models are limited to 2D and marred with over-simplifying physical assumptions of the physics that govern EM damage. In our recent study, where a three-dimensional phase-field model to simulate EM-mediated defects has been introduced, the study is limited to columnar polycrystals where the defect pathways are predetermined.4 However, predicting the damage pathway is not as trivial in interconnects with non-columnar microstructures, which typically comprise a random distribution of grains. Therefore, previous phase-field studies need to be extended in order to gain a complete understanding of how electric field intensity, grain size, and grain boundary misalignment with respect to the applied electric field influence the failure mechanisms.

Building on our earlier work,3,4 we begin with isolating the effect of grain boundary misalignment with respect to the applied electric field in simulations of slit formation. Then, we simulate the EM-mediated defect evolution in an equiaxed polycrystalline microstructure and analyze the mechanisms by which slit formation and surface drift occur. Next, we examine both these defect evolution modes in a realistic setting comprising the Voronoi distribution of grains. Finally, the influence of three-dimensional capillarity on the defect evolution kinetics is explored. The overall goal is to understand the influence of the electric field and the grain size on the morphological evolution of EM-mediated defects.

This paper is organized as follows. In Sec. II, we revisit the phase-field model for simulating EM-mediated evolution of defects in a polycrystalline interconnect, which comprises an arbitrary number of grains. In Sec. III, we explore the influence of GB misalignment with respect to the current flow, grain size, electric field, and three-dimensional capillarity on the EM-mediated damage in polcrystalline interconnects comprising non-columnar grain microstructures. Section IV concludes the article.

In this study, we investigate a periodic array of equiaxed hexagonal grains representing a polycrystalline interconnect, shown in Fig. 1. The region labeled as the underlayer represents the separation between the interconnect and the dielectric layer. The phase-field model to simulate EM-mediated defects’ growth in interconnects comprising equiaxed and random distributions of grains leverages the model already reported in our previous works except that both those studies were limited to analyzing the damage evolution in columnar grains.3,4 In this context, we would like to clarify that the following model, whose applicability is limited to the meso-length-scale, assumes the interconnect to be free of any pre-existing point defects of lattice dimensions.

FIG. 1.

Schematic diagram illustrating the initial configuration of the interconnect, which comprises an equiaxed network of grains, used for phase-field simulations. To distinguish individual grains within the interconnect region, it is sub-divided using the non-conserved order parameter ηi where i{1,2,3,4}. A conserved order parameter, ρ, which represents the density, transitions from 0 to 1 across the underlayer/interconnect surface. Grain triple points (TPs) and the boundaries (GBs) are indicated. The grain size is determined by the distance between circled vertices of hexagonal grains. While no-flux boundary conditions are assigned at the top and bottom edges of the computation domain, the left and right boundaries are periodic. The externally applied electric field is directed normally into the underlayer/interconnect surface.

FIG. 1.

Schematic diagram illustrating the initial configuration of the interconnect, which comprises an equiaxed network of grains, used for phase-field simulations. To distinguish individual grains within the interconnect region, it is sub-divided using the non-conserved order parameter ηi where i{1,2,3,4}. A conserved order parameter, ρ, which represents the density, transitions from 0 to 1 across the underlayer/interconnect surface. Grain triple points (TPs) and the boundaries (GBs) are indicated. The grain size is determined by the distance between circled vertices of hexagonal grains. While no-flux boundary conditions are assigned at the top and bottom edges of the computation domain, the left and right boundaries are periodic. The externally applied electric field is directed normally into the underlayer/interconnect surface.

Close modal

The conserved order parameter, ρ, maps the density field such that ρ=1.0 corresponds to the interconnect, whereas ρ=0.0 represents the underlayer. The non-conserved parameter, ηi, distinguishes between individual grains such that ηi equals 1.0 within the corresponding grain and 0.0 elsewhere. The density, ρ, is defined such that at the grain boundary, there is a small decrease to reflect the decreased density of the grain boundary. The free energy density functional is given by

F=V[f(ρ,ηi)+κρρ2+inκηηi2+fem]dV,
(1)

where f(ρ,ηi) represents the bulk free energy, the middle two terms dependent on ρ and ηi correspond to the diffuse surface and grain boundary energies, respectively, while the final term accounts for the electrostatic free energy associated with the externally applied electric field. κρ and κη are the gradient energy coefficients for ρ and η, respectively. The bulk free energy is given by

f(ρ,ηi)=Aρ2(1ρ)2+Bρ2ξ(ηi)+C(1ρ)2inηi2,
(2)

where A, B, and C are bulk energy coefficients that define the well height in the bulk free energy-order parameter space, while ξ(ηi) is defined as

ξ(ηi)=in[η44η22+2j>inηi2ηj2]+0.25.
(3)

The bulk free energy expression is chosen to ensure free energy minima at (ρ=0,ηi=0,ηj=0), (ρ=1,ηi=1,ηj=0), and (ρ=1,ηi=0,ηj=1) where “i” and “j” represent grains in the interconnect. The electrostatic free energy is given by

fem=eNAzρϕ,
(4)

where e represents the coulombic charge on an electron, NA is Avogadro’s number, z is the effective valence of the species undergoing EM, and ϕ denotes the potential field as defined in Eq. (10). Temporal evolution of the density field, ρ, is governed by a Cahn–Hillard equation as a Cahn–Hillard modeled approach,25,26

ρt=Mμ,
(5)

where μ is defined as

μ=δFδρ,
(6)

such that the temporal evolution is given by

ρt=M(ρ,ηi)[f(ρ,ηi)ρ2κρ2ρ+zeNAϕ].
(7)

The atomic mobility, M(ρ,ηi), is chosen as

M(ρ,ηi)=Mb+4Mgbj>i[ηi2ηj2]1/2+16Msρ2(1ρ)2
(8)

to account for the variability in atomic mobilities along the surface of the interconnect (Ms), grain boundaries (Mgb), as well as in the bulk (Mb). The kinetic equation for the evolution of the grain field, η, is described by a system of Allen–Cahn equations,27 

ηit=LδFδηi=L[f(ρ,ηi)ηi2κη2ηi],
(9)

where L is the grain relaxation parameter, which determines the coarsening rate. The electrostatic field distribution is calculated by solving a Laplace equation,

[σ(ρ)ϕ]=0,
(10)

where σ, which represents the electrical conductivity, is linearly interpolated in terms of ρ as

σ(ρ)=σintρ+σdie(1ρ),
(11)

with σint and σdie representing the conductivity of the metallic interconnect and the underlayer region, which is composed of a dielectric material, respectively.

Non-dimensionalization of all the above quantities is done by selecting a characteristic energy scale E=A, timescale t=L2(ρintρdie)2MbE, and length scale L=(κρA)1/2.

The coupled system of Cahn–Hillard and Allen–Cahn is solved iteratively using the finite difference method. Grid spacing along both the x and y directions (Δx=Δy=1.0) is assumed to be equal, and the time step is represented by Δt. A Neumann (no-flux) boundary condition is imposed along the y-direction while periodic boundaries in the x-direction for effectively allowing a periodic array of interconnects. To investigate the influence of the grain size on the morphological evolution of defects, different domains of sizes 107Δx×100Δx, 143Δx×141Δx, 171Δx×173Δx, and 194Δx×200Δx are considered. The equiaxed grain size is defined as the length to join the opposite corners of the hexagon, encircled in Fig. 1. For studies involving random grain distributions, the interconnect domain is composed of 171Δx×173Δx with 20Δx assigned to the underlayer domain, which is in contact with the polycrystalline interconnect. This implicitly assigns an average grain size of 102Δx.

Among other external parameters such as the current density, our previous phase-field simulations of EM in columnar grains3,4 suggest that the defects’ evolution is governed by a ratio of atomic mobilities along the surface and the grain boundaries.3,4 Therefore, one method to mitigate defect proliferation is to minimize grain boundaries by increasing the grain size. Since interconnects are typically polycrystalline, they offer several pathways for the diffusion to occur. With this prior knowledge, we are able to identify a parameter space of interest, which is described below.

Simulation parameters are selected to isolate different defect modes as well as to inhibit grain coarsening when applicable, as listed in Table I. The three defect modes investigated in this study are slit-formation, surface drift, and mixed modes. Isolated evolution of slits with little to no grain coarsening is necessary to develop a baseline for future studies involving more complex environments such as slit formation with grain coarsening, anisotropic interface energies, or anisotropic grain boundary atomic mobilities.

TABLE I.

Parameters used for simulations reported in Sec. III.

SimulationGrain size (Δx)LMbMgb/MsΔϕκηκρσintσdie
Influence of GB misalignment N/A 0.1 10−6 107 0.4 0.33 1.0 3.0 0.3 
Influence of grain sizes 59, 83, 102, 117 0.1 10−6 107 0.4 0.33 1.0 3.0 0.3 
Influence of electric fields 102 0.1 10−6 107 0.4, 0.8, 1.2, 1.6, and 2.0 0.33 1.0 3.0 0.3 
Mixed modes 59 0.1 10−6 103, 10, and 4 0.4 0.33 1.0 3.0 0.3 
Randomly distributed grains N/A 0.001 and 0.01 10−6 107, 106, and 10 0.4 0.33 1.0 3.0 0.3 
SimulationGrain size (Δx)LMbMgb/MsΔϕκηκρσintσdie
Influence of GB misalignment N/A 0.1 10−6 107 0.4 0.33 1.0 3.0 0.3 
Influence of grain sizes 59, 83, 102, 117 0.1 10−6 107 0.4 0.33 1.0 3.0 0.3 
Influence of electric fields 102 0.1 10−6 107 0.4, 0.8, 1.2, 1.6, and 2.0 0.33 1.0 3.0 0.3 
Mixed modes 59 0.1 10−6 103, 10, and 4 0.4 0.33 1.0 3.0 0.3 
Randomly distributed grains N/A 0.001 and 0.01 10−6 107, 106, and 10 0.4 0.33 1.0 3.0 0.3 

The ratio of atomic mobilities along the surface (Ms), grain boundaries (MGB), and bulk (Mb) were varied to explore the role of preferential diffusion pathways. Grain coarsening was controlled by assuming a small relaxation parameter, L. The ratio, κρ/κη=3.0, was chosen to assign a dihedral angle of 139.6°. Although a change in the dihedral angle can alter the surface curvature near the groove region, thereby influencing the morphological evolution, as seen in Fig. 2, it was left unaltered to keep the parameter space of interest manageable.

FIG. 2.

A plot depicting the influence of the dihedral angle on the kinetics of slit formation. Both simulations were run with a Mgb/Ms of 107, L=0.1, and Δϕ=0.4. The dihedral angle of 139.6° corresponds to a κρ/κη ratio of 3 and is utilized in this present work. A κρ/κη ratio of 1 generates a dihedral angle of 121.8°, and the disparate kinetics can be visualized above. The inset represents the initial condition of a bicrystalline interconnect with current flow perpendicular to the surface.

FIG. 2.

A plot depicting the influence of the dihedral angle on the kinetics of slit formation. Both simulations were run with a Mgb/Ms of 107, L=0.1, and Δϕ=0.4. The dihedral angle of 139.6° corresponds to a κρ/κη ratio of 3 and is utilized in this present work. A κρ/κη ratio of 1 generates a dihedral angle of 121.8°, and the disparate kinetics can be visualized above. The inset represents the initial condition of a bicrystalline interconnect with current flow perpendicular to the surface.

Close modal

Starting from a two-dimensional symmetric grain configuration shown in the inset of Fig. 3, we leverage the phase-field model described in Sec. II A to simulate and compare the relative displacements of groove roots at distinct GB misalignment ranging from 0° to 70° with respect to the direction of the electric field. One can trivially predict that the groove propagation velocity at finite GB misalignment (θ0) will scale as sec(θ). However, we observe that the average slope (ΔxΔt) differs by a factor of 3.606 when comparing 0° and 60° of misalignment with respect to current flow. This outcome demonstrates that the groove velocity scales non-linearly and can be influenced by a number of factors, including grain boundary misalignment with respect to the direction of current flow.

FIG. 3.

A plot comparing the temporal displacement of groove roots for distinct misalignment between the grain boundary and the applied current, θ (shown in the inset), from 0° of misalignment to 70°. The inset shows the initial configuration wherein surface and grain boundaries are represented by thin- and thick-solid lines, respectively. The dotted line represents the direction along which an external electric current is applied, i.e., normal with respect to the interconnect surface.

FIG. 3.

A plot comparing the temporal displacement of groove roots for distinct misalignment between the grain boundary and the applied current, θ (shown in the inset), from 0° of misalignment to 70°. The inset shows the initial configuration wherein surface and grain boundaries are represented by thin- and thick-solid lines, respectively. The dotted line represents the direction along which an external electric current is applied, i.e., normal with respect to the interconnect surface.

Close modal

Here, we investigate the influence of grain sizes on the morphological evolution of grain boundary slits under EM. Simulated evolution of the slit can be viewed in Figs. 4(a)4(d). With grain boundary atomic diffusivity 7 orders of magnitude above that of the surface and the bulk, the surface drift observed is negligible.4 Due to the presence of equiaxed grains, the propagation of the slit is symmetrical along the GBs.

FIG. 4.

Progression of slit formation at intermittent time steps: (a) t=5000 , (b) t=9000, (c) t=15000, and (d) t=21000. A mobility ratio (Mgb/Ms) of 107 was used to simulate slit formation. The relaxation parameter, L, was set to 0.1 in order to suppress grain coarsening, and a potential difference, Δϕ=0.4, is applied such that electric field lines are pointing normally into the interconnect surface. The grain size, as per the definition in Fig. 1, is 102Δx.

FIG. 4.

Progression of slit formation at intermittent time steps: (a) t=5000 , (b) t=9000, (c) t=15000, and (d) t=21000. A mobility ratio (Mgb/Ms) of 107 was used to simulate slit formation. The relaxation parameter, L, was set to 0.1 in order to suppress grain coarsening, and a potential difference, Δϕ=0.4, is applied such that electric field lines are pointing normally into the interconnect surface. The grain size, as per the definition in Fig. 1, is 102Δx.

Close modal

The temporal evolution of the GB slit, corresponding to a grain size of 102Δx, is displayed in Figs. 4(a)4(d). Slit formation also occurs in simulations with a lower grain boundary to surface atomic mobility ratios albeit with comparatively slower growth rates. However, surface drift starts to occur at Mgb/Ms ratios that are lesser than 100, as discussed in Sec. III D.

We observe that the average slit formation velocity decreases with an increase in the grain size. To characterize the kinetics of slit formation wherein the GB alignment with respect to the electric field varies within the polycrystalline interconnect, we plot the temporal evolution of the path fraction in Fig. 5(a). The path fraction at an instant is calculated as the corresponding length of the GB slit, which is measured with respect to the initial position of the interconnect/underlayer surface, divided by the total GB length within the interconnect. The grains here being perfect hexagons, the relative misalignment of the GB can either be 0° or 60° with respect to the direction of the electric field, which results in variable slit growth rates.

FIG. 5.

(a) Slit path fraction plotted as a function of time highlights the influence of the grain size on the velocity of slit propagation. The path fraction is used to normalize the slit length by the maximum displacement of the slit tip before it touches the bottom edge of the computational domain. (b) Temporal variation in the average velocity of the slit tip depicts the influence of the grain size on defect kinetics. Simulations reported in (a) and (b) correspond to a mobility ratio (Mgb/Ms) of 107, which facilitates slit formation; relaxation parameter, L, is set to 0.1 to suppress grain coarsening; and a potential difference, Δϕ=0.4, is applied such that electric field lines are pointing normally into the interconnect surface. Grain sizes are chosen as per the definition in Fig. 1. The plot in (c) depicts the influence of the increasing electric field magnitude on the temporal displacement of the slit tip in a network of equiaxed grains each of size, 102Δx. For the sake of comparison, other material-specific parameters are kept consistent with respect to (a) and (b).

FIG. 5.

(a) Slit path fraction plotted as a function of time highlights the influence of the grain size on the velocity of slit propagation. The path fraction is used to normalize the slit length by the maximum displacement of the slit tip before it touches the bottom edge of the computational domain. (b) Temporal variation in the average velocity of the slit tip depicts the influence of the grain size on defect kinetics. Simulations reported in (a) and (b) correspond to a mobility ratio (Mgb/Ms) of 107, which facilitates slit formation; relaxation parameter, L, is set to 0.1 to suppress grain coarsening; and a potential difference, Δϕ=0.4, is applied such that electric field lines are pointing normally into the interconnect surface. Grain sizes are chosen as per the definition in Fig. 1. The plot in (c) depicts the influence of the increasing electric field magnitude on the temporal displacement of the slit tip in a network of equiaxed grains each of size, 102Δx. For the sake of comparison, other material-specific parameters are kept consistent with respect to (a) and (b).

Close modal

A decrease in the slit propagation velocity in the region of the angled grain boundary is inversely proportional to the grain size, as shown in Fig. 5(b). The presence of larger grains prolongs failure, i.e., the time required for the slit to penetrate completely into the network of grains.

However, in the context of interconnects, instances when grain boundaries align perfectly with current flow are rare. In Fig. 5(a), the slit propagation rate, represented by the slope, is slightly slower for the first vertical GB when compared to the two sets of vertical grain boundaries following a triple point. The difference between the slopes is attributed to temporally increasing the distance between the groove root and the surface exposed to the underlayer region, which contributes to a back-fill of the groove roots. As the groove root descends down inside the interconnect, the surface atoms need to progressively diffuse over a larger distance to replenish the groove root. In the first GB section, this distance, however, is smaller, which leads to a reduced growth rate. As the slit descends down, the growth rate increases due to reduction in the back-fill flux, which depends on the relative distance between the groove root and surface through the chemical potential gradient dependence.

In this section, we explore the role of the electric field intensity on the kinetics of slit propagation. It is known that the current density scales with voltage. In the present model, conductivity, σ, is interpolated in terms of the density, as shown in Eq. (11). Thus, we can investigate the influence of the current density by increasing the potential difference, Δϕ, across the interconnect while keeping the grain size unchanged.

In Fig. 6, we visualize the electric field gradient vectors as the slit approaches and traverses past a triple point. The electric field gradient vector is calculated by multiplying the electrostatic free energy density [Eq. (4)] and the conductivity at every point on the grid, finding its gradient, and then using vector addition to obtain a net gradient vector.

FIG. 6.

Electric field gradient vectors plotted at intermittent time steps: (a) t=5000, (b) 9000, (c) 15 000, and (d) 21 000. Sky blue is used to represent the vacuum phase (ρ=0), while the red lines correspond to grain boundaries between the equiaxed grains. For the simulation, the following parameters are used: Mgb/Ms=107, Δϕ = 0.4, and L=0.1.

FIG. 6.

Electric field gradient vectors plotted at intermittent time steps: (a) t=5000, (b) 9000, (c) 15 000, and (d) 21 000. Sky blue is used to represent the vacuum phase (ρ=0), while the red lines correspond to grain boundaries between the equiaxed grains. For the simulation, the following parameters are used: Mgb/Ms=107, Δϕ = 0.4, and L=0.1.

Close modal

In Fig. 7, we illustrate the temporal changes in the gradient electric field vector, which correspond to Fig. 6, plotted at time steps 2000, 3000, and 6000. The dynamics of slit formation can be visualized in the vicinity of the triple point. Current flow drives the mass along the surface and toward the grain boundary, as seen in Fig. 7(a), which facilitates slit formation along the grain boundary. It is worth noting that subsequent transport of atoms down the grain boundary is largely limited by the value of Mgb. In Fig. 7(b), it can be observed that the current flow, which was initially pushing the atoms along the grain boundary, is assisting the transport into the bulk portion of the grains. As the slit tip progresses, transport into the bulk continues, as seen in Fig. 7(c). This justifies a decrease in the slit formation propagation velocity as plotted in Figs. 5(a)5(c). Comparing the current flow in Figs. 7(a) and 7(c) where the grains are equiaxed, it can be concluded that the wrapping of electric field gradient vectors around the slit influences the growth kinetics. Therefore, in interconnects with non-columnar grain microstructures, non-homogeneous current flow around the slit tip can cause variable growth rates, which will be discussed further in Sec. III E.

FIG. 7.

Magnified view of electric field gradient vectors displayed in Fig. 6, showing behavior of the electric field as the slit approaches the triple point. (a) The slit is approaching the triple point at t=2000. (b) At t=3000, the groove root reaches the triple point, and (c) at t=6000, the slit traverses past the triple point where it splits. Sky blue is used to represent the vacuum phase (ρ=0), while the red lines correspond to grain boundaries between the equiaxed grains.

FIG. 7.

Magnified view of electric field gradient vectors displayed in Fig. 6, showing behavior of the electric field as the slit approaches the triple point. (a) The slit is approaching the triple point at t=2000. (b) At t=3000, the groove root reaches the triple point, and (c) at t=6000, the slit traverses past the triple point where it splits. Sky blue is used to represent the vacuum phase (ρ=0), while the red lines correspond to grain boundaries between the equiaxed grains.

Close modal

In Fig. 5(c), we plot the temporal displacement of the slit tip at different electric field strengths. It can be observed that the tip velocity increases with the strength of the electric field. Sudden changes in the slope can be attributed to misalignment of GB with respect to the applied electric field, similar to the plot shown in Fig. 5(a).

Mixed mode failure caused by electromigration involves slit formation as well as surface drift. Simulations with a Mgb/Ms ratio less than 100 depict mixed mode failure as the grain boundary atomic flux is replenished through mass transport from the surface of the interconnect into the groove root. The diffusional mechanism by which such a failure occurs has previously been reported3 for columnar grains. Here, we analyze

the defect propagation rate in a network of non-columnar grains, by varying the ratio,Mgb/Ms.

We simulate the morphological evolution corresponding to Mgb/Ms ratios of 4, 10, and 1000 to compare the failure rate observed in three cases, as shown in Figs. 8(a), 8(b), and 8(c), respectively. As the Mgb/Ms ratio decreases, the slit evolution rate becomes more sluggish; however, a relative increase in Ms results in smoother corners along the slit, which intensifies surface diffusion along the slit edges, resulting in widening of slits. Surface drift is more prominent in simulations with comparable Mgb and Ms. This is due to an occurrence of enhanced surface flux into the groove walls, resulting in more mass transport into the grain boundary. While diffusion along grain boundaries is not sufficiently high for a slit to form, a lower Mgb also results in slower defect propagation.

FIG. 8.

Comparison of defect modes at time steps 1500, 5000, and 10 000 for Mgb/Ms ratios of (a) 4, (b) 10, and (c) 1000. Parameters Δϕ=0.4, L=0.1, while the grain size equals 59Δx.

FIG. 8.

Comparison of defect modes at time steps 1500, 5000, and 10 000 for Mgb/Ms ratios of (a) 4, (b) 10, and (c) 1000. Parameters Δϕ=0.4, L=0.1, while the grain size equals 59Δx.

Close modal

In Figs. 9(a)9(d), we report the morphological evolution simulated at a Mgb/Ms ratio of 10, where a characteristic surface drift, in the direction of the applied electric field, is observed. The mechanism by which this phenomenon occurs has been reported in our previous work.3,4 The phenomenon of surface drift under current stressing was shown to converge toward a global steady state28,29 using the level-set method, which is faithfully reproduced in the present simulations, as indicated by the linear portions of the plot in Fig. 10(a) that have a non-zero slope. Furthermore, the grain size can have a profound influence on the rate of surface drift, as seen in Fig. 10(a). From these plots that correspond to different sizes, two distinct steady states are identified. Surface drift occurs at a nearly constant rate until the top two grains are depleted [Fig. 9(b)]. The direction of diffusion flux that originates at the grain surface and points toward the groove root and subsequently down the grain boundaries is depicted in Fig. 9(a). Figures 9(b) and 9(d) correspond to the attainment of local steady states following the consumption of grains exposed to the surface. In both these instances, the surface atomic fluxes are able to replenish the atoms within the triple-point regions, which were otherwise transported away by the GB fluxes. This results in little to no movement of the groove root relative to the translating surface, a phenomenon that has previously been classified as surface drift in the experimental studies.28,29 The progression of the defect when the number of grains reduces to two [Fig. 9(d)] obeys the same kinetics as reported by Mukherjee et al.3 Failure caused by surface drift is much slower as compared to slit formation, which requires significantly greater mass transport, as depicted in Fig. 10(b).

FIG. 9.

Defect propagation at Mgb/Ms=10, shown at intermittent time steps: (a) t=5250, (b) t=20250, (c) t=60250, and (d) t=100250. The grain size equals 83Δx, assigned as per the definition in Fig. 1. For these simulations, L=0.1 and Δϕ=0.4.

FIG. 9.

Defect propagation at Mgb/Ms=10, shown at intermittent time steps: (a) t=5250, (b) t=20250, (c) t=60250, and (d) t=100250. The grain size equals 83Δx, assigned as per the definition in Fig. 1. For these simulations, L=0.1 and Δϕ=0.4.

Close modal
FIG. 10.

(a) Temporal evolution of the slit path fraction at an Mgb/Ms ratio of 10 plotted at different grain sizes. Grain sizes are chosen as per the definition in Fig. 1. (b) Temporal evolution of the slit path fraction at Mgb/Ms ratios of 107 and 10 to compare the two defect mode slit formation and surface drift, respectively. The grain size equals 102Δx, L=0.1 and Δϕ=0.4.

FIG. 10.

(a) Temporal evolution of the slit path fraction at an Mgb/Ms ratio of 10 plotted at different grain sizes. Grain sizes are chosen as per the definition in Fig. 1. (b) Temporal evolution of the slit path fraction at Mgb/Ms ratios of 107 and 10 to compare the two defect mode slit formation and surface drift, respectively. The grain size equals 102Δx, L=0.1 and Δϕ=0.4.

Close modal

In this section, we examine slit formation and surface drift in a more realistic polycrystalline setting, which comprises a random distribution of grains. Starting from the Voronoi distribution of grains, we simulate the morphological evolution of defects assuming the same set of simulation parameters described in Sec. II. The relaxation parameter, L, from Eq. (9) is assumed to be small (L=0.001) to suppress grain coarsening.

Simulation snapshots plotted at intermittent times in Fig. 11 differentiate the morphological evolution at distinct Mgb/Ms ratios. The overall failure rate is sluggish at Mgb/Ms=10, as opposed to when the ratio is 106 times larger due to a transition from surface drift to a slit formation regime. The nature of the slit propagation on the right panel is found to be strongly dependent on GB alignment. In Fig. 11(b), it is observed that the growth of slits after branching at the triple points, designated by the arrows, along the misaligned GB becomes sluggish or comes to a standstill. On the contrary, slit branches evolving along GBs, which are more favorably aligned with respect to the applied electric field, continue to progress at the rate determined by relative misalignment, as discussed in Fig. 3. In Fig. 12, the gradient vector plot composed of arrows allows us to visualize the temporal flow of current around the groove roots as the slits proliferate into the polycrystalline interconnect. Since the GBs of grain designated by η1 are almost equally aligned with respect to current flow, the adjacent groove displacements are nearly equal. On the contrary, slit progression between the GB of η4 and η5 is slow. This slowing down of slit progression is characterized by the wrapping of gradient vectors around the slit tip, which inhibits the growth. EM-mediated defects’ evolution kinetics is thus heavily dependent on grain boundary misalignment with respect to current flow. Therefore, complementary future experiments that can validate the mechanisms proposed through computational studies are desirable.

FIG. 11.

Phase-field simulations of (a) surface drift corresponding to Mgb/Ms=10 and (b) slit formation corresponding to Mgb/Ms=107 in an interconnect comprising the Voronoi distribution of grains. For both the simulations, the relaxation parameter, L, is set to 0.001. The potential difference across the interconnect, Δϕ, is set to 0.4.

FIG. 11.

Phase-field simulations of (a) surface drift corresponding to Mgb/Ms=10 and (b) slit formation corresponding to Mgb/Ms=107 in an interconnect comprising the Voronoi distribution of grains. For both the simulations, the relaxation parameter, L, is set to 0.001. The potential difference across the interconnect, Δϕ, is set to 0.4.

Close modal
FIG. 12.

Electric field gradient vectors corresponding to the microstructural evolution shown in Fig. 11 plotted at intermittent time steps: (a) t=11000, (b) t=34000, and (c) t=54000.

FIG. 12.

Electric field gradient vectors corresponding to the microstructural evolution shown in Fig. 11 plotted at intermittent time steps: (a) t=11000, (b) t=34000, and (c) t=54000.

Close modal

To investigate the influence of three-dimensional capillarity on EM-mediated defect evolution in interconnects comprising non-columnar grains, we extend our two-dimensional calculations reported above to 3D. For this purpose, we generated a random distribution of three grains using the Voronoi algorithm in a three-dimensional domain with periodic boundary conditions along the x and y axes and no-flux along the z-direction. For comparison with 2D, several equivalent two-dimensional domains comprising distinct distributions of three grains are generated such that the average grain size is retained and is represented by a gray region in Fig. 13(e). Preserving the grain size in 2D and 3D facilitates a direct comparison of the growth kinetics while mitigating the possibility of grain size effects that can influence defects’ evolution (see Sec. III B). The number of two-dimensional domains analyzed was kept sufficiently large to minimize statistical errors. Since our aim is to compare the morphological evolution of slits, Δϕ is set to 0.4, L=0.01, while the mobility ratio (Mgb/Ms) is assumed to be 106.

FIG. 13.

Three-dimensional evolution of slits in a network of (a) three randomly distributed grains. Snapshots comparing the partially failed interconnect microstructures comprising a network of (b) 3, (c) 4, and (d) 5 grains at t=1620 are also shown. (e) Plot of temporal evolution of the path fraction corresponding to simulations reported in (a)–(d) wherein the shaded region represents the mean rate of failure in two-dimensional simulations while accounting for the standard deviation. Differences in the evolution kinetics highlight the need for four-dimensional characterization of EM-mediated defects. For all the simulations, Δϕ=0.4, L=0.01, and Mgb/Ms=106.

FIG. 13.

Three-dimensional evolution of slits in a network of (a) three randomly distributed grains. Snapshots comparing the partially failed interconnect microstructures comprising a network of (b) 3, (c) 4, and (d) 5 grains at t=1620 are also shown. (e) Plot of temporal evolution of the path fraction corresponding to simulations reported in (a)–(d) wherein the shaded region represents the mean rate of failure in two-dimensional simulations while accounting for the standard deviation. Differences in the evolution kinetics highlight the need for four-dimensional characterization of EM-mediated defects. For all the simulations, Δϕ=0.4, L=0.01, and Mgb/Ms=106.

Close modal

The temporal evolution of slits in a three-dimensional periodic array of randomly oriented grains can be seen in Fig. 13(a). The shaded region within Fig. 13(e) is the mean slit evolution rate for the two-dimensional simulations, which accounts for the standard deviation. While the slit evolution rates in 3D as well as 2D, on an average, follow a linear trend, the slope of the former is found to approximately three times larger than the latter when considering grains of similar sizes. This can be attributed to a multitude of diffusional pathways that are available in 3D as opposed to 2D where the evolution is constrained along a line. The plot in Fig. 13(e) shows that the three-dimensional failure rate increases with an increase in the GB density or with grain refinement with a direct comparison displayed in Figs. 13(b)13(d). As the number of grains are increased without altering the size of the three-dimensional computational domain, the discrepancies in predicted failure rates is found to decrease, which indicates the importance of large scale studies. Finally, while an accurate prediction of defects’ evolution rate warrants large scale studies comprising thousands of randomly oriented grains, in all certainties, our findings highlight the pressing need for future characterization studies in 4D.

Growth of electromigration-mediated defects in polycrystalline interconnects comprising non-columnar grains was examined by a phase-field model. Defect modes studied were slit formation, surface drift, and mixed modes. Below are the main conclusions from this work:

  1. In a network of equiaxed grains, the grain size, as expected, follows an inverse proportionality with the slit propagation rate. GB alignment has a significant influence on the defect kinetics; when the GB is perfectly aligned along the direction of current flow, the gradient vectors wrap around the slit tip, resulting in a mass transport flux directing down the GB. However, increasing the magnitude of GB misalignment with respect to an externally applied electric field results in more electromigration into the bulk, rather than along the GB, thereby retarding the slit propagation velocity.

  2. When the mobility ratio, MGB/Ms, is less than 100, the mixed mode, which imbibes the characteristics of both surface drift as well as slit growth, occurs, significantly lowering the defect propagation rates. This is primarily caused by the mass transport along the surface of the interconnect, which replenishes the groove root. Curving of the interconnect surface occurs when the MGB/Ms ratios are small, which further increases the replenishment of the groove root.

  3. Phase-field simulations of defect evolution in a network of randomly distributed grains reveal that GBs that are equally misaligned with respect to the external electric field result in an identical slit propagation rate, while a severe misalignment can halt the defect progression altogether. While this finding is corroborated through our three-dimensional simulations, the slit formation propagation is found to be several times faster when compared with 2D simulations. This increase in the defect evolution velocity is due to the availability of numerous diffusional pathways in 3D as opposed to 2D where the defect is constrained to propagate along a line. Such a difference in growth kinetics emphasizes the need for four-dimensional characterization of EM-mediated defects for their efficient mitigation.

The authors acknowledge funding from the National Science Foundation (NSF) under Grant No. NSF CMMI-1763128 (Dr. Alexis Lewis, Program Manager). Fruitful discussions with Dr. Nikhilesh Chawla’s group at the Arizona State University are gratefully acknowledged.

1.
H.
Ceric
and
S.
Selberherr
,
Mater. Sci. Eng. Rep.
71
,
53
(
2011
).
2.
P. S.
Ho
and
T.
Kwok
,
Rep. Prog. Phys.
52
,
301
(
1989
).
3.
A.
Mukherjee
,
K.
Ankit
,
R.
Mukherjee
, and
B.
Nestler
,
J. Electron. Mater.
45
,
6233
(
2016
).
4.
A.
Mukherjee
,
K.
Ankit
,
M.
Selzer
, and
B.
Nestler
,
Phys. Rev. Appl.
9
,
044004
(
2018
).
5.
ITRS
, “International Technology Roadmap for Semiconductors,” Technical Report, 2011.
6.
T. O.
Ogurtani
and
O.
Akyildiz
,
J. Appl. Phys.
97
,
093520
(
2005
).
7.
O.
Akyildiz
and
T.
Ogurtani
,
J. Appl. Phys.
110
,
043521
(
2011
).
8.
T. O.
Ogurtani
and
E. E.
Oren
,
J. Appl. Phys.
96
,
7246
(
2004
).
9.
W.
Li
,
C. M.
Tan
, and
Y.
Hou
,
J. Appl. Phys.
101
,
104314
(
2007
).
10.
C.
Basaran
and
M.
Lin
,
Int. J. Mater. Struct. Integr.
1
,
016039
(
2007
).
11.
K.
Croes
,
Y.
Li
,
M.
Lofrano
,
C. J.
Wilson
, and
Z.
Tokei
, in
2013 IEEE International Reliability Physics Symposium (IRPS)
(IEEE, 2013), pp. 2C.3.1–2C.3.4.
12.
F.
He
and
C. M.
Tan
,
World J. Model. Simul.
8
,
271284
(
2012
).
13.
S.
Chakraborty
,
P.
Kumar
, and
A.
Choudhury
,
Acta Mater.
153
,
377
(
2018
).
14.
C. K.
Hu
,
R.
Rosenberg
, and
K. Y.
Lee
,
Appl. Phys. Lett.
74
,
2945
(
1999
).
15.
D. A.
Porter
,
K. E.
Easterling
, and
M.
Youssef Abdelraouf Sherif
,
Phase Transformations in Metals and Alloys
(
CRC Press
,
2009
).
16.
E.
Baibuz
,
S.
Vigonski
,
J.
Lahtinen
,
J.
Zhao
,
V.
Jansson
,
V.
Zadin
, and
F.
Djurabekova
,
Comp. Mater. Sci.
146
,
287
(
2018
).
17.
D.
Contestable-Gilkes
,
D.
Ramappa
,
M.
Oh
, and
S. M.
Merchant
,
J. Electron. Mater.
31
,
1047
(
2002
).
18.
E.
Glickman
and
M.
Nathan
,
J. Appl. Phys.
80
,
3782
(
1996
).
19.
M. S.
Park
and
R.
Arróyave
,
Acta Mater.
58
,
4900
(
2010
).
20.
M. S.
Park
,
S. L.
Gibbons
, and
R.
Arróyave
,
Acta Mater.
61
,
7142
(
2013
).
21.
M.
Mahadevan
and
R.
Bradley
,
Phys. Rev. B
59
,
11037
(
1999
).
22.
J.
Santoki
,
A.
Mukherjee
,
D.
Schneider
,
M.
Selzer
, and
B.
Nestler
,
J. Electron. Mater.
48
,
182
(
2019
).
23.
L.
Klinger
and
L.
Levin
,
J. Appl. Phys.
78
,
1669
(
1995
).
24.
A. F.
Bower
and
S.
Shankar
,
Mod. Sim. Mater. Sci. Eng.
15
,
923
(
2007
).
25.
J. W.
Cahn
and
J. E.
Hilliard
,
J. Chem. Phys.
28
,
258
(
1958
).
26.
J. W.
Cahn
,
Acta Metall. Mater.
9
,
795
(
1961
).
27.
S. M.
Allen
and
J. W.
Cahn
,
Acta Metall. Mater.
27
,
1085
(
1979
).
28.
M.
Nathan
,
E.
Glickman
,
M.
Khenner
,
A.
Averbuch
, and
M.
Israeli
,
Appl. Phys. Lett.
77
,
3355
(
2000
).
29.
M.
Khenner
,
A.
Averbuch
,
M.
Israeli
,
M.
Nathan
, and
E.
Glickman
,
Comput. Mater. Sci.
20
,
235
(
2001
).