Subsections

• 4.5.1 Theoretical Formulation
• 4.5.2 The Diffuse Interface Model
• 4.5.3 Numerical Implementation
• 4.5.4 Setting of the Initial Order Parameter Profile and Initial Grid Refinement $\Lambda _{h}(0)$
• 4.5.5 Finite Element Scheme
• 4.5.6 Maintaining the Grid during Simulation
• 4.5.7 Grid Adaptation
• 4.5.8 Solving Procedure

4.5 Void evolution and Interconnect Resistance Change Models

The description of the time period up to void nucleation must be directly followed by a model which handles void evolution. Once a void is nucleated, material transport on the void metal interface dominates [54,64,60,61,63]. The driving forces are electromigration proportional to the local tangential component of the electric field and the chemical potential. Since the most critical voids are nucleated at the copper/barrier and copper/cap-layer interfaces, which also represent the fastest diffusion paths, available models [64,60] have to be extended in order to include these phenomena.
Void evolution was originally handled applying explicit void surface tracking methods [64]. A further development was introducing numerically efficient diffuse interface models [56,57,58,59]. However, not until introducing an adaptive mesh approach into the diffuse interface method [62], realistic interconnect structures could be handled.
Applying a diffuse interface model we are able to predict the intrinsic void behavior from the nucleation until failure. The simulation of the evolving void/metal interface is accompanied with electrical field calculation and resistance extraction.

Once the void is nucleated, the dominating factor for the void evolution is the self-diffusion on the void surface and in the vicinity of this surface. The driving forces are again the chemical potential and the tangential component of the electrical field gradient. The chemical potential in this case is comprised of components due to the free energy of the surface and the local elastic strain energy. The material is transported along the void surface and into the metal bulk. This effect gives rise to the movement of the whole void and it's shape change. In this context the crystal orientation must be considered because of the anisotropic characteristics of self-diffusion [54].

Since the problem is a moving boundary problem, the first attempts to handle it was the application of explicit tracking of the void/metal interface by means of finite element methods [64,60]. These methods are quite complicated to implement and also have rather poor numerical stability. To overcome the problems mentioned above, we apply a diffuse interface model. In a diffuse interface model void and metal area are presented through an order parameter $\phi$ which takes values $+1$ in the metal area, $-1$ in the void area and $-1 < \phi < +1$ in the void-metal interface area. Demanding explicit tracking of the void-metal interface is not needed and models do not require boundary conditions to be enforced at the moving boundary.

The evolving void/metal interface is described through the time dependent distribution of the order parameter $\phi$. The actual diffuse interface is a belt area around the iso-surface $\phi=0$. This belt area has some finite width. With shrinking width of this belt area the accuracy of the simulation increases. The diffuse interface governing equation (DIGE) is a central part of each diffuse interface based modeling. The equation itself comprises electromigration and stress driven material transport on the void/metal interface and on this basis regulates the distribution of the order parameter. This equation is an extension of the well-known Chan-Hilliard equation [56,57,58,85,86]. Concerning the solution of Chan-Hilliard type equations several approach are proposed in the literature [85,86]. The gap between fast and reliable DIGE solving approaches which are already known from the literature and a solution approach which is also applicable for arbitrary interconnect geometries was bridged by introducing an adaptive grid technique [62].

The main idea of the adaptive grid method was to use local grid refinement and a coarsening algorithm, a fine grid is set dynamically adjusted only in the vicinity of the void/metal interface. The geometry of the interconnect can be arbitrary and, due to the economical distribution of the grid density, computation can be performed efficiently.

Additional equations to be solved such as the Laplace equation for the electric field and mechanical stress-strain model equations are adapted to the diffuse interface model through an order parameter dependence of the electrical conductivity and the elastic compliance [59]. These equations are solved using common finite element schemes assuming a finite width of the diffuse interface area. Generally the accuracy of the calculation will increase by reducing this width. There is a rigorous asymptotic theory which proves that the diffuse interface model converges to the sharp void/interface physical situation both in the equations' solution and in the compliance with boundary conditions [56,57,58,59].

The initial void is placed on the spot predicted by the first stage module and the second stage module is invoked, which simulates evolution of the void under the impact of driving forces, until it reaches a critical size which changes the initial interconnect resistivity to such an extent that it can be considered as actual interconnect failure.
Thereby the resistance change up to the interconnect failure and, the time needed for this failure to occur are provided. We call this time period void evolution time and added to the void nucleation time calculated in the first stage module, it defines the interconnect time to failure.

4.5.1 Theoretical Formulation

We assumed unpassivated monocrystal isotropic interconnects where stress phenomena can be neglected. An interconnect is idealized as two-dimensional electrically conducting via which contains initially a circular void. For simplicity we also neglect the effects of grain boundaries and lattice diffusion. In this case there are two main forces which influence the shape of the evolving void interface: the chemical potential gradient and electron wind. The first force causes self-diffusion of metal atoms on the void interface and tends to minimize energy which results in circular void shapes. The electron wind force produces asymmetry in the void shape depending on the electrical field gradients.
Including contributions from both, electromigration and chemical potential-driven surface diffusion, gives the total surface atomic flux, $\bf {J_{A}}$ = $J_{A}\bf {t}$, where $\bf {t}$ is the unit vector tangent to the void surface [60,87]

$$J_A = D_s (-\vert e\vert Z^{*}E_s+\gamma_{s}\Omega \bf {t}\cdot\nabla_{s}\kappa)$$ (4.41)

$Z^{*}$ is the effective valence (Section 4.2), $e$ is the charge of an electron, $E_s\equiv \bf {E_s}\cdot{t}$ is the local component of the electric field tangential to the void surface, $\kappa$ is the local surface curvature, and $\nabla_s$ is the surface gradient operator; $\kappa \equiv \nabla \cdot \bf {n}$, where $\bf {n}$ is the local unit vector normal to the void surface. Further, $\gamma_{s}$ is the surface energy, $\Omega$ is the volume of an atom, and $D_s$ is given through an Arrhenius law:

$$D_s=\frac{D_0\delta_s}{k T} \Bigl(\frac{-Q_s}{k T}\Bigr)$$ (4.42)

Here, $\delta_s$ is the thickness of the diffusion layer, $T$ is the temperature, $Q_s$ is the activation energy for the surface diffusion, and $D_0$ is the pre-exponential factor for mass diffusion. Equation (4.41) is the Nernst-Einstein equation, where the sum in the parentheses on the right side expresses the driving force. Mass conservation gives the void propagation velocity normal to the void surface, $v_n$, through the continuity equation [61,63,87],

$$v_n = -\Omega \cdot \nabla_s \bf {J_A}.$$ (4.43)

The electric field $\bf {E}$ can be written as the gradient of the electric potential, i.e.,

$$\bf {E} \equiv -\nabla \varphi.$$ (4.44)

The field is also solenoidal, i.e.,

$$\nabla \cdot \bf {E} = 0.$$ (4.45)

Equations (4.44) and (4.45) imply that the potential $\varphi$ obeys Laplace's equation,

$$\Delta \varphi = 0.$$ (4.46)

The void surface and the interconnects edges are modeled as insulating boundaries, i.e.

$$\nabla \varphi \cdot \bf {n}=0,$$ (4.47)

at these boundaries.

4.5.2 The Diffuse Interface Model

The diffuse interface model equations for the void evolving in an unpassivated interconnect line are the balance equation for the order parameter $\phi$ [56,57,58,59],

$$\frac{\partial \phi} {\partial t} = \frac{2D_{s}} {\epsilon \pi} \nabla \cdot (\nabla \mu - \vert e\vert Z^{*} \nabla \varphi),$$ (4.48)

$$\mu = \frac{4\Omega \gamma_{s}}{\epsilon \pi}(f'(\phi) - \epsilon^{2} \Delta \phi),$$ (4.49)

and for the electrical field,

$$\nabla \cdot (\gamma_E(\phi) \nabla \varphi)= 0.$$ (4.50)

where $\mu$ is the chemical potential, $f(\phi)$ is the double obstacle potential as defined in [85,86], and $\epsilon$ is a parameter controlling the void-metal interface width. The electrical conductivity was taken to vary linearly from the metal ( $\gamma_E = \gamma_{E,metal}$) to the void area ( $\gamma_E = 0$) by setting $\gamma_E = \gamma_{E,metal} (1+\phi)/2$. The equations (4.48) and (4.50) are solved on the two-dimensional polygonal interconnect area $T$.
It has been proven [57,59] that in the asymptotic limit for $\epsilon\rightarrow 0$ the diffuse interface model defined by equations (4.48)-(4.50) describes the same voids-metal interface behavior like equations (4.41)-(4.46). The width of the void-metal diffuse interface is approximately $\epsilon\pi/2$, and in order to numerically handle sufficiently thin interfaces one needs a very fine locally placed grid around it.

4.5.3 Numerical Implementation

The equations (4.48)-(4.50) are solved by means of the finite element method on a sequence of the grids $\Lambda_{h}(t_{0}=0),\Lambda_{h}(t_{1}),\Lambda_{h}(t_{2})$ each one adapted to the position of the void-metal interface from the previous time step. The initial grid $\Lambda _{h}(0)$ is produced by refinement of some basic triangulation $T_h$ of area $T$ according to the initial profile of order parameter $\phi$. The motivation of grid adaptation is to construct and maintain a fine triangulated belt of width $\epsilon\pi$ in the interconnect area where $-1 < \phi < +1$, respectively, where the void-metal interface area is placed.

4.5.4 Setting of the Initial Order Parameter Profile and Initial Grid Refinement $\Lambda _{h}(0)$

The initial order parameter profile depends on the initial shape of the void $\Gamma(0)$ and can be expressed as

$$\phi(x,y)= \begin{cases}+1& \text{if} \quad d > \frac{\epsilon\pi... ...{\epsilon\pi}{2}, -1& \text{if} \quad d < -\frac{\epsilon\pi}{2} \end{cases}$$ (4.51)

Where $d=dist(P(x,y),\Gamma(0))$ is the signed normal distance of the point from the initial interface $\Gamma(0)$. To obtain sufficient resolution of this initial profile, the basic grid $T_{h}$ is transformed into grid $\Lambda _{h}(0)$ obeying the following initial grid refinement criterion (IGRC) for the circular void with center $O$ and radius $r$:

$$\vert dist(C_{K},O)-r\vert\leq \frac{\epsilon\pi}{2} \quad \wedge \quad h_{K} > \frac{\epsilon\pi}{n}$$ (4.52)

$n$ is the chosen number of grid elements across the void-metal interface width, $h_{K}$ is the longest vertex of the triangle $K$, and $C_{K}$ is its center of gravity. Now an adaptive algorithm defined in Section 4.5.7 transforms the basic grid $T_{h}$ into an initial grid $\Lambda _{h}(0)$ according to $IGRC$.

4.5.5 Finite Element Scheme

$\Lambda_{h}(t_{n})$ is a triangulation of the area $T$ at discrete time $t_{n}$, and $\{\phi_{i}^{n-1}\}_{i=0}^{N-1}$ are discrete nodal values of the order parameter on this triangulation. A finite element based iteration for solving (4.48) on grid $\Lambda_{h}(t_{n})$ and evaluating the order parameter at the time $t_{n+1}=t_{n}+\Delta t$ consists of two steps [88]:

Step 1. For the $k^{th}$ iteration of the $n+1^{th}$ time step the linear system of equations has to be solved:

$$\epsilon \frac{\pi}{2} M_{ii}\phi_{i}^{n+1,k}+\Delta t D_{s}K_{ii}\mu_{i}^{n+1,k}=\alpha_{i}$$ (4.53)

$$M_{ii}\mu_{i}^{n+1,k}-\tau\biggl(\frac{1}{\epsilon}M_{ii}+\epsilon K_{ii}\biggr)\phi_{i}^{n+1,k}=\beta_{i},$$ (4.54)

where

$$\alpha_{i}=\epsilon \frac{\pi}{2}M_{ii}\phi_{i}^{n}-\Delta t D_{s}\sum_{i\neq j}K_{ij}\mu_{j}^{n+1,k-1}$$ (4.55)

$$\beta_{i}=\tau\epsilon \sum_{i\neq j} K_{ij}\phi_{j}^{n+1,k-1}-\vert e\vert Z^{*}M_{ii}V_{i}^{n}$$ (4.56)

for each $i=0,N-1$ of the nodal values $(\phi^{n+1}_{i},\mu^{n+1}_{i})$ of the triangulation $\Lambda_{h}(t_{n})$. $[M]_{ij}$ and $[K]_{ij}$ are the lumped mass and stiffness matrix, respectively and $\tau = \frac{4 \Omega \gamma_{s}}{\pi}$.

Step 2. All nodal values $\{\phi_{i}^{n+1}\}_{i=0}^{N-1}$ are projected on $[-1, 1]$ by a function

$$\rho(x)=max(-1,min(1,x)).$$ (4.57)

For solving (4.50) a conventional finite element scheme is applied [76].

4.5.6 Maintaining the Grid during Simulation

After an order parameter was evaluated on the $\Lambda_{h}(t_{n})$ a grid needs to be adapted according to the new void-metal interface position. Therefore it is necessary to extract all elements which are cut by the void-metal interface in grid $\Lambda_{h}(t_{n})$. The following condition is used: We take a triangle $K \in \Lambda_{h}(t_{n})$ and denote its vertices as $P_{0},P_{1},P_{2}$. The triangle $K$ belongs to the interfacial elements if for the values of the order parameter $\phi$ at the triangle's vertices holds $\phi(P_{0})\phi(P_{1})<0$ or $\phi(P_{1})\phi(P_{2})<0$. We assume that an interface intersects each edge of the element only once. The set of all interfacial elements at the time $t_{n}$ is denoted as $E(t_{n})$. The centers of gravity of each triangle from the $E(t_{n})$ build the interface point list $L(t_{n})$. The distance of the arbitrary point $Q$ from $L(t_{n})$ is defined as,

$$dist(Q,L(t_{n}))= \underset{P \in L(t_{n})}{min} dist(Q,P).$$ (4.58)

Thus we can define the transitional grid refinement criterion TGRC

$$dist(C_{K},L(t_{n})) \leq \frac{\epsilon\pi}{2}\quad \wedge \quad h_{K} > \frac{\epsilon\pi}{n}$$ (4.59)

The grid adaption for the next time step evaluation of the order parameter $\phi$ is now completed with respect to $TGRC$ as defined in the next section.

Figure 4.6: Atomic refinement operation.

4.5.7 Grid Adaptation

The grid adaption algorithm used in this work is a version of the algorithm introduced in [89] and consists of a grid refinement algorithm and a grid coarsening algorithm. The refinement algorithm is based on recursive bisecting of triangles. A triangle is marked for refinement if it complies with some specific refinement criterion $COND$. In our application we used $COND=IGRC$ for the initial refinement and $COND=TGRC$ for the grid maintaining refinement. For every triangle of the grid, the longest one of its edges is marked as refinement edge. The element and its neighbor element which also contains the same refinement edge are refined by bisecting this edge (Fig. 4.6). We can define refinement of the element in the following way:

Now, the overall refinement algorithm can be formulated as follows:

The parameter max_refinement_depth limits the number of bisecting of a triangle. The grid is refined until there is no more element marked for refinement or the maximal refinement depth is reached.
The coarsening algorithm is more or less the inverse of the refinement algorithm. Each element that does not fullfil criterion $COND$ is marked for coarsening. The basic idea is to find the father element whose refinement produced the element in consideration.

The following routine coarsens as many elements as possible, even more than once if allowed:

The complete adaption of the grid is reached by sequential invoking of the grid refinement and of the grid coarsening algorithm.

4.5.8 Solving Procedure

Combining all the components described in the Section 4.5.4 - Section 4.5.7 we obtain the complete void evolution solving procedure (Fig. 4.7).

Figure 4.7: Void evolution solving procedure.