3.5 Model Summary

The model equations presented in the chapter are implemented in a TCAD tool (COMSOL Multiphysics^® [38]) in order to simulate the electromigration failure phenomenon in realistic 3D interconnect structures. Since electromigration modeling and simulation constitute a multiphysics problem, they can be conveniently separated in several submodels. In the following a brief summary of the general model, with the main model equations which are numerically solved for the simulations, is presented.

Electro-Thermal Model

The distributions of current density j, electric potential V_e, and temperature T in the interconnect line are determined by the solutions of the non-linear system of equations

\[\begin{equation} \nabla \cdot \vec{j}=0, \end{equation}\] (3.96)

\[\begin{equation} \nabla^2 V_\text{e} =0, \end{equation}\] (3.97)

\[\begin{equation} \nabla^2 T -\cfrac{\rho_\text{m}c_\text{p}}{k_\text{t}} \frac{\partial T}{\partial t} = - \cfrac{\sigma_\text{e}}{k_\text{t}}(\nabla V_\text{e})^2, \end{equation}\] (3.98)

as presented in Section 3.1.

Vacancy Dynamics Model

The vacancy transport J_v, responsible for electromigration failure is induced by different driving forces and is described as follows

\[\begin{equation} \vec J_\text{v} = -D_\text{v}\left(\nabla C_\text{v}+\cfrac{C_\text{v}|Z^*|e}{k_\text{b}T}\nabla V_\text{e}+\cfrac{C_\text{v}Q^*}{k_\text{b}T^2}\nabla T+\cfrac{C_\text{v}f \Omega_\text{a}}{k_\text{b}T}\nabla \sigma\right). \end{equation}\] (3.99)

The change of vacancy concentration C_v caused by the vacancy flux in the structure is expressed by the continuity equation

\[\begin{equation} \frac{\partial{C_\text{v}}}{\partial t}=-\nabla \cdot\vec J_\text{v}+G. \end{equation}\] (3.100)

The effect of the fast diffusivity path n, which affects the vacancy dynamics inside an interconnect line, is taken into account by setting appropriate diffusion coefficients D_vⁿ and source/sink terms as follows

\[\begin{equation} G^\text{n}=\cfrac{1}{\tau^\text{n}_\text{v}}\left(C^\text{n}_{\text{v}\text{,eq}}-C^\text{t}_\text{v}\left(1+\cfrac{\omega_\text{r}}{\omega_\text{t}C_\text{v}}\right)\right). \end{equation}\] (3.101)

By following the model described in Section 3.2, the vacancy concentration distribution in the structure can be determined.

Solid Mechanics Model

It has been shown in Section 3.3 that the mechanical stress σ and displacement u distributions in the line are obtained by solving the following system of equations

\[\begin{equation} \cfrac{\partial \epsilon^\text{v}_\text{ik}}{\partial t} =\cfrac{\Omega_\text{a}}{3}\left((1-f)\nabla \cdot \vec{J_\text{v}}\pm f G\right)\delta_\text{ik}, \end{equation}\] (3.102)

\[\begin{equation} \sigma_\text{ik}=\lambda\delta_\text{ik}(\epsilon_\text{ll}-\epsilon^\text{th}_\text{ll}-\epsilon^\text{v}_\text{ll})+2\mu(\epsilon_\text{ik}-\epsilon^\text{th}_\text{ik}-\epsilon^\text{v}_\text{ik}), \end{equation}\] (3.103)

\[\begin{equation} \nabla \cdot \overline{\overline\sigma}=0, \end{equation}\] (3.104)

\[\begin{equation} \epsilon_\text{ik}=\cfrac{1}{2}\left(\cfrac{\partial u_\text{i}}{\partial x_\text{k}}+\cfrac{\partial u_\text{k}}{\partial x_\text{i}}\right), \end{equation}\] (3.105)

\[\begin{equation} \mu\nabla^2 u_\text{i}+(\lambda+\mu)\cfrac{\partial}{\partial x_\text{i}} (\nabla \cdot \vec u)= B \cfrac{\partial \mathrm{Tr}(\epsilon^\text{th}+\epsilon^\text{v})}{\partial x_\text{i}}. \end{equation}\] (3.106)

From the time evolution of the stress build-up due to electromigration and the accompanying driving forces in the line, the void nucleation condition can be determined. The void nucleates in the structure as soon as the stress reaches the threshold value

\[\begin{equation} \sigma_\text{thr}= \cfrac{2\gamma_\text{$m$}\sin\theta_\text{c}}{R_\text{p}}. \end{equation}\] (3.107)

Void Evolution Model

The diffuse interface model and a semi-empirical model are applied to a structure with an already-nucleated void for the simulation of the electromigration-induced void evolution mechanism, as presented in Section 3.4. The diffuse interface model provides the distribution of the order parameter φ in the interconnect by solving the set of equations

\[\begin{equation} \cfrac{\partial \phi}{\partial t} = -\ \nabla \cdot \vec{J_\text{s}}, \end{equation}\] (3.108)

\[\begin{equation} \vec{J_\text{s}}= -\cfrac{2 D_\text{s}(\phi)}{\epsilon_\text{di} \pi}\left(\nabla\mu_\text{s}-e|Z^*|\nabla V_\text{e}\right). \end{equation}\] (3.109)

The development of a semi-empirical model, based on the void size dependence of the incoming flux of vacancies due to electromigration, provides numerical solutions of the time t needed to grow a void of a given void radius r_v as follows

\[\begin{equation} t = t_\text{0} + \cfrac{\pi}{ \alpha}\int_{r_\text{0}}^{r} \cfrac{r_\text{v}^2}{A_\text{i}(r_\text{v})C_\text{v}(r_\text{v})||{\vec{j}(r_\text{v})}||}\, \text{d}{r_\text{v}}, \end{equation}\] (3.110)

where

\[\begin{equation} \alpha = f\Omega_\text{a} \cfrac{e|Z^*|D_\text{v}\rho_\text{e}}{k_\text{b}T}. \end{equation}\] (3.111)

The model equations are used to determine the change in the interconnect resistance with time, due to the growing void in the structure. Consequently, the prediction of the interconnect lifetime can be obtained.