Subsections

• 4.6.1 Two-Dimensional Void Evolution
• 4.6.2 Estimating the Void Growth Time and Resistance Change
• 4.6.3 A Study of Electromigration Promoting Factors and Vacancy Dynamics

4.6 Simulation Results

There are three major questions to be answered by TCAD regarding electromigration reliability issues. First, interconnect layout designers need to know, how long a drafted interconnect structure will work properly, e.g. for the given operating conditions, how long does it take for complete electromigration failure to happen. The second question of interest is, how does the interconnect resistance change due to electromigration damage. And, the third interesting question is how does the local picture of all electromigration promoting factors looks like in the case when electromigration damage takes place.

The thorough knowledge of the physical phenomena behind the electromigration failure and the capability to model and simulate these phenomena is the only way to accomplish interconnect layouts which optimally meet technology requirements and at the same time avoid far to restrictive design rules.

In the following sections we will present application of the introduced physical models and numerical concepts, realized in the simulation tools, on the different aspects of the electromigration reliability problem.

4.6.1 Two-Dimensional Void Evolution

In all following simulations a circle was chosen as initial void shape. The resolution of the parameter $\phi$ profile can be manipulated by setting parameter $n$ which is the mean number of triangles across the void-metal interface. On Fig. 4.8 initial grids for $n=1$ and $n=5$ are presented.

Figure 4.8: Initial grid refinements.

$n=1$

$n=5$

We consider a two-dimensional, stress free, electricaly conducting interconnect via.

Figure 4.9: Interconnect via with initial void.

$A$ constant voltage is applied between points $A$ and $B$ (Figure 4.9). At $B$ a refractory layer is assumed. Because of geometrical reasons there is current crowding in the adjacencies of the corners $C$ and $D$. The analytical solution of equation (4.50) has at these points actually a singularity [90,76].
High electrical field gradients in the area around the corner points increase overall error of the finite element scheme for the equation (4.50) which is overcome by applying an additional refinement of the basic mesh $T_h$ according to the local value of the electric field gradient (Figure 4.10).

Figure 4.10: Profile of the current density (in $A/m^2$) at the corners of the interconnect.

A fine triangulated belt area which is attached to the void-metal interface at the initial simulations step follows the interfacial area throughout the simulations whereby the interconnect area outside the interface is coarsened to the level of the basic grid $T_{h}$ (Figure 4.11).

Figure 4.11: Refined grid around the void in the proximity of the interconnect corner.

Figure 4.12: Void evolving through interconnect in the electric current direction

Figure 4.13: Time dependent resistance change during void evolution for the different initial void radius $r$.

In our simulations a void evolving through the straight part of the interconnect geometry exhibits similar shape changes as observed in earlier models [60,59]. There is also no significant fluctuation of the resistance during this period of interconnect evolution. The situation changes when the void evolves in the proximity of the interconnect corner (Figure 4.12). Due to current crowding in this area the influence of the electromigration force on the material transport on the void surface is more pronounced than the chemical potential gradient. This unbalance leads to higher asymmetry in the void shape then observed in the straight part of the interconnect (Figure 4.12). A void evolving in the proximity of the interconnect corner causes significant fluctuations in the interconnect resistance due to void asymmetry and position. The resistance change shows a charasteristic profile with two peaks and a valley (Figure 4.13). The extremes are more pronounced for the larger initial voids.
The capability of the applied adaptation scheme is also presented in the simulation of void collision with the interconnect refractory layer (Figure 4.14).
The time step $\Delta t$ for the numerical scheme (4.53)-(4.56) is fitted at the simulations begin taking into account inverse proportionality of the speed of the evolving void-metal interface to the initial void radius [60]:

$$\Delta t = \frac{\epsilon\pi r l}{2D_{s}\vert e\vert Z^{*}V_{0}}$$ (4.60)

Figure 4.14: Grid adaptation in the case of void collision with the refractory layer.

$l$ is the characteristical length of via geometry. An appropriate choice of the time step ensures that the evolving void-metal interface will stay inside the fine grid belt during the simulation. The dynamics of the evolving void-metal interface simulated with a the presented numerical scheme complies with the mass conservation law, the void area (where $\phi=-1$) remains approximately the same during the whole simulation. Notable area deviations during the simulation appear only, if a relatively large factor $\epsilon$ has been chosen. As scaling length we took $l=1 \mu m$ and for the initial void radius $r_0=0.035 l$, $r_1=0.045 l$, and $r_2=0.065 l$. Our simulations have shown that for all considered initial void radii, voids follow the electric current direction (Figure 4.13) and do not transform in slit or wedge like formations which have been found to be a main cause for a complete interconnect failure [54]. Already with $\epsilon = 0.003 l$ good approximations are achieved. The number of elements on the cross section of the void-metal interface was chosen between 6 and 10 with the interface width of $0.0015 l\pi$.

4.6.2 Estimating the Void Growth Time and Resistance Change

An initial void with some small radius $r_0$ is placed at some characteristic position inside the three-dimensional interconnect structure (Figure 4.15). Since most of the fatal voids are nucleated in the vicinity or in the area of interconnect vias we consider in particular these cases. The configurable initial void volume is $V_0$ which is smaller than $4 \pi r_0^3/3$ because the void area is confined by sphere and boundary of the interconnect (Figure 4.15). Starting from the initial void radius $r_0$, the void radius is gradually incremented $r_0, r_1=r_0+\Delta r_1, r_2=r_1+\Delta r_2,...$, with $\Delta r_1\geq \Delta r_2\geq ...\geq\Delta r_n$. For each void radius the electrical field in the interconnect structure is calculated by means of the finite element method using a diffuse interface approach. To obtain the distribution of the electrical potential inside the interconnects the Poisson equation (4.50) has to be solved.

In order to obtain sufficient accuracy the scalar field $\phi(x,y,z,t)$ must be resolved on a locally refined mesh (Figure 4.16). For an electrical field calculated in such a way, the resistance of the interconnect via is also calculated [76,91]. With growing void size the resistance increases. The whole process is stopped when a void radius is reached for which $100\times (R_{actual}/R_{initial}-1)> 20 \%$.

Figure 4.15: Position of the growing void. Initial position and volume are chosen on the basis of experimental results.

The primary driving force of material transport at the void surface is electromigration proportional to the tangent component of the vector current density. Since the diffuse interface approach for the calculation of the current density ensures physical behavior of the electrical field in the vicinity of the isolating void, the normal component of the current density on the void surface is always zero and we can apply the formula

$$J_{m,i}=2\frac{\int_{V}\Vert\mathbf{J}\Vert[1-\phi^2_i(x,y,z)] dV}{\int_{V}[1-\phi^2_i(x,y,z)] dV},$$ (4.61)

for the average current density over a void with radius $r_i$. (4.61) expresses the averaging of the current density weighted with finite element volume inside the interconnect. Since $\phi_i(x,y,z)=1$ in metal and $\phi_i(x,y,z)=-1$ in the void area, the term $1-\phi^2_i(x,y,z)$ is non-zero only in void-metal interface area.

Figure 4.16: Typical current density distribution picture in the vicinity of the spherical void. The red area marks peak values of the current density. The mesh is locally refined in the void-metal interface.

The evolution of the void is caused by material transport on the void surface and in the vicinity of the void surface. The mass conservation law gives the mean propagation velocity $v_i$ of the evolving void-metal interface

$$v_i=\frac{D_V}{kT\rho(\phi)}eZ^{*} J_{m,i},$$ (4.62)

here $D_v$ is the vacancy diffusivity and $Z^{*}$ effective charge number of vacancies. (4.62) is valid for all void shapes.

Figure 4.17: Change of the average current density and via resistance depending on the void radius.

As we can see from Figure 4.17, the average current density on the void surface increases with the void size. Both, current density and resistance, exhibit a very similar dynamic behavior. The dynamical resistance increase is in accordance with the measurement results presented in [52]. Compared with the earlier result [92], which assumes cubical void shapes, our approach enables more realistic simulations. An open question is how to use the obtained average current density (Figure 4.17) for the estimation of the void growing time ($t_E$) up to the critical void size. In [92] a simple formula is applied

$$t_{E}=\frac{V_c-V_0}{v_m A_s}.$$ (4.63)

In this equation $V_c$ is the critical void size, $V_0$ is the initial void size, $A_s$ is the cross section of the interconnect in the vicinity of the growing void, and $v_m$ is the mean velocity of the evolving void-metal interface. However, this formula is only valid in the case of cubical void which is a very rough approximation of the real situation. According to newer experimental results [8] the real void shape is significantly better approximated by a spherical sector. In this case $t_E$ can be estimated as

$$t_{E}=\sum_{i}\frac{\Delta r_{i+1}}{v_{i}},$$ (4.64)

assuming that for sufficiently small $\Delta r_{i+1}$, the void radius grows from $r_i$ to $r_{i+1}$ with a constant velocity $v_i$.

As we can see from ( $\ref{est0}$), the velocity $v_i$ depends on the vacancy diffusivity $D_v$ which itself has significantly varying values depending on the diffusion path. The electromigration assisted self-diffusion of copper is a complex process which includes simultaneous diffusion through the crystal bulk, along grain boundaries, along the copper/barrier interfaces, and along the copper/cap-layer interface. Therefore, the diffusivity applied in ( $\ref{est0}$) must be a cumulative value (see equation 4.12). For a feasible estimation of $t_E$, reliable, experimentally determined values for all relevant diffusivities are needed and these are until now unfortunately not available [55,8].

4.6.3 A Study of Electromigration Promoting Factors and Vacancy Dynamics

In this section we present the simulation of the influence of electromigration promoting factors on the electromigration behavior for the case of a complex interconnect structure produced with advanced process simulation tools.

The results of the transient electro-thermal simulation (Section 4.3.2) are used to setup a thermo-mechanical (Section 4.3.2) and material transport simulation. The model applied is similiar to that presented in Section 4.4.5 but without taking into account the connection between the local vacancy dynamic and overall mechanical stress distribution.

An effective vacancy diffusivity is chosen assuming the copper/TiN barrier interface as the dominant diffusion path [55].

An interconnect structure with 6 metal (copper) lines crossing 4 other lines at a lower level and connected with vias/contacts was generated using a typical damascene process flow. This was performed by using a combination of DEP3D [93] (only for deposition of TiN barrier layer and copper deposition for the second metal layer and for dielectric deposition of the topmost layer) with DEVISE [94] (used for emulating the other deposition, etching and polishing steps). The barrier layer is 0.02 $\mu$m thick; each metal line in the first layer is 3 $\mu$m long with a width of 0.25 $\mu$m and aspect ratio of 1.2. For the second layer, each metal line is 2.2 $\mu$m long with a width 0.3 $\mu$m and aspect ratio 1.5. The contacts/vias have a width of 0.2 $\mu$m with an aspect ratio of 2.5. The dielectric medium surrounding the metal lines is silicon-dioxide. The contacts for electrical characterization were defined using DEVISE.

The simulation of the electro-thermal problem was carried out with the SAP [91] tool and the simulation of mechanical stress and electromigration problem with FEDOS tool. The simulation results are presented in the Figure 4.18 - Figure 4.20. We can see in Figure 4.18 electrical potential distribution. The gradient of the electrical potential determines direction and intensity of the electromigration driving force.

The corresponding temperature distribution is presented in Figure 4.19. The temperature has twofold influence on the electromigraton. First, temperature changes diffusion coefficeints on all diffusion paths according to Arrhenius' law. Secondly, for the case that due to specific interconnect layout and operating conditions a substantial gradients in the temperature distribution are available, these gradients represents additional driving force which can promote or opose electromigration (4.32). Additionally, thermal mismatch between copper, barrier and silicon-dioxide passivation layer causes a build-up of thermo-mechanical stress.

The peak values of the vacancy concentration indicate the spots where serious material transport discontinuity takes place. These hot spots are the locations where the void nucleation can be expected (Figure 4.20). Analysis of the experimental results [55,95,77] shows that the material transport discontunity can be intialised by the geometrical properties of layout itself, like in the cases where the segments of the copper interconnects are divided by barrier layer, and promoted due to electromigration unbalanced by other driving fources. A significant part plays also the grain boundary structure of the copper acting as a network of fast diffusivity paths.

During the operating of an interconnect generally the redistribution of vacancies under the influence of the promoting factors takes place. Starting from the uniformly distributed vacancies, during the operating of the interconnect local peaks of the vacancy concentration are built (Figure 4.20).

If the balance of promoting factors and temperature gradients which characterize operating conditions of the interconnect is reached, the vacancy concentration will remain at some value independent of the simulation duration: the interconnect structure is virtually immortal. However, if the electrical field is so high that electromigration overrides the material reaction presented through concentration and mechanical stress gradients, the vacancy concentration will increase steadily and at some point (e.g. vacancy concentration reaches unrealistic high values) the applied model is no more valid. If the vacancy concentration is higher than the lattice sites concentration throughout some volume of the macroscopic dimensions (for example $1/10$ of the interconnect width), we can assume that void is already nucleated and for further simulation void evolution models should be applied (Sections 4.6.1 and 4.6.2). The interconnect layout and operating conditions choosen in our simulation example produce a situation where the vacancy concentration remains stable at $16 \cdot C_V^{eq}\approx 10^{17}$cm$^{-3}$ and under the lattice sites concentration ($10^{23}$ cm$^{-3}$).

Figure 4.18: Electrical potential distribution [mV].

Figure 4.19: Temperature distribution in silicon-dioxide passivation layer [K].

Figure 4.20: Electromigration induced vacancy distribution relativ to the equilibrum distribution. Peak values indicate spots of the possible void nucleation.