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.

## I. INTRODUCTION

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\xb0$ 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 reliability^{16,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 voids^{24} 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.

## II. METHODS

### A. Phase-field model

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.

The conserved order parameter, $\rho $, maps the density field such that $\rho =1.0$ corresponds to the interconnect, whereas $\rho =0.0$ represents the underlayer. The non-conserved parameter, $\eta i$, distinguishes between individual grains such that $\eta i$ equals 1.0 within the corresponding grain and 0.0 elsewhere. The density, $\rho $, 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

where $f(\rho ,\eta i)$ represents the bulk free energy, the middle two terms dependent on $\u2207\rho $ and $\u2207\eta 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. $\kappa \rho $ and $\kappa \eta $ are the gradient energy coefficients for $\rho $ and $\eta $, respectively. The bulk free energy is given by

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

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

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 $\varphi $ denotes the potential field as defined in Eq. (10). Temporal evolution of the density field, $\rho $, is governed by a Cahn–Hillard equation as a Cahn–Hillard modeled approach,^{25,26}

where $\mu $ is defined as

such that the temporal evolution is given by

The atomic mobility, $M(\rho ,\eta i)$, is chosen as

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, $\eta $, is described by a system of Allen–Cahn equations,^{27}

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

where $\sigma $, which represents the electrical conductivity, is linearly interpolated in terms of $\rho $ as

with $\sigma int$ and $\sigma 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\u2032=A\u2032$, timescale $t\u2032=L\u20322(\rho int\u2032\u2212\rho die\u2032)2Mb\u2032E\u2032$, and length scale $L\u2032=(\kappa \rho \u2032A\u2032)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 ($\Delta x=\Delta y=1.0$) is assumed to be equal, and the time step is represented by $\Delta 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\Delta x\xd7100\Delta x$, $143\Delta x\xd7141\Delta x$, $171\Delta x\xd7173\Delta x$, and $194\Delta x\xd7200\Delta 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\Delta x\xd7173\Delta x$ with $20\Delta x$ assigned to the underlayer domain, which is in contact with the polycrystalline interconnect. This implicitly assigns an average grain size of $102\Delta x$.

### B. Parameter selection

Among other external parameters such as the current density, our previous phase-field simulations of EM in columnar grains^{3,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.

Simulation . | Grain size (Δx)
. | L
. | M_{b}
. | M_{gb}/M_{s}
. | Δϕ
. | κ_{η}
. | κ_{ρ}
. | σ_{int}
. | σ_{die}
. |
---|---|---|---|---|---|---|---|---|---|

Influence of GB misalignment | N/A | 0.1 | 10^{−6} | 10^{7} | 0.4 | 0.33 | 1.0 | 3.0 | 0.3 |

Influence of grain sizes | 59, 83, 102, 117 | 0.1 | 10^{−6} | 10^{7} | 0.4 | 0.33 | 1.0 | 3.0 | 0.3 |

Influence of electric fields | 102 | 0.1 | 10^{−6} | 10^{7} | 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} | 10^{3}, 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} | 10^{7}, 10^{6}, and 10 | 0.4 | 0.33 | 1.0 | 3.0 | 0.3 |

Simulation . | Grain size (Δx)
. | L
. | M_{b}
. | M_{gb}/M_{s}
. | Δϕ
. | κ_{η}
. | κ_{ρ}
. | σ_{int}
. | σ_{die}
. |
---|---|---|---|---|---|---|---|---|---|

Influence of GB misalignment | N/A | 0.1 | 10^{−6} | 10^{7} | 0.4 | 0.33 | 1.0 | 3.0 | 0.3 |

Influence of grain sizes | 59, 83, 102, 117 | 0.1 | 10^{−6} | 10^{7} | 0.4 | 0.33 | 1.0 | 3.0 | 0.3 |

Influence of electric fields | 102 | 0.1 | 10^{−6} | 10^{7} | 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} | 10^{3}, 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} | 10^{7}, 10^{6}, 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, $\kappa \rho /\kappa \eta =3.0$, was chosen to assign a dihedral angle of $139.6\xb0$. 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.

## III. RESULTS AND DISCUSSION

### A. Influence of GB alignment on slit formation

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\xb0$ to $70\xb0$ with respect to the direction of the electric field. One can trivially predict that the groove propagation velocity at finite GB misalignment ($\theta \u22600$) will scale as sec($\theta $). However, we observe that the average slope $(\Delta x\Delta t)$ differs by a factor of 3.606 when comparing $0\xb0$ and $60\xb0$ 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.

### B. Influence of grain sizes on GB slit formation

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.

The temporal evolution of the GB slit, corresponding to a grain size of $102\Delta 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\xb0$ or $60\xb0$ with respect to the direction of the electric field, which results in variable slit growth rates.

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.

### C. Influence of electric fields on GB slit formation

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, $\sigma $, 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, $\Delta \varphi $, 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.

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.

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).

### D. Mixed modes

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 reported^{3} 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.

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 state^{28,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).

### E. Defect evolution in a network of randomly oriented grain boundaries

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 $\eta 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 $\eta 4$ and $\eta 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.

### F. Influence of three-dimensional capillarity

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, $\Delta \varphi $ is set to 0.4, $L=0.01$, while the mobility ratio ($Mgb/Ms$) is assumed to be $106$.

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.

## IV. CONCLUSIONS

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:

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.

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.

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.

## ACKNOWLEDGMENTS

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.