3 Discretization of the Three-Stream Mulvaney-Richardson Model

First we consider dopant continuity equation of Mulvaney-Richardson model (3.45). The local change of dopants is due to the divergence of fluxes of dopant-interstitial and dopand vacancy pairs,

$$\frac{\partial C_A}{\partial t}=\nabla\cdot(\mathbf{J}_{AI}+\mathbf{J}_{AV}).$$ (154)

In this model the diffusion coefficient $D_A$ of the point defect-dopant pairs is considered constant and we can write (based on (3.45)),

$$\nabla\cdot\mathbf{J}_{AI}= \frac{f_I D_A}{C_{I}^{eq}}\Delta (C_I C_A) + \frac{f_I D_A}{C_{I}^{eq}}\nabla(C_I C_A \nabla ( n)),$$ (155)

$$\nabla\cdot\mathbf{J}_{AV}=\frac{(1-f_I) D_A}{C_{V}^{eq}}\Delta (C_V C_A) + \frac{(1-f_I) D_A}{C_{V}^{eq}}\nabla(C_V C_A \nabla ( n)).$$ (156)

In order to obtain a weak formulation of (3.92) on the single element $T\in T_h(\Omega)$ we multiply both sides of the equation with the basis function $\varphi_k$, and integrate over $T$,

$$\int_T\nabla\cdot\mathbf{J}_{AI} \varphi_k d\Omega = -\frac{f_... ...a \varphi_k d\Omega +\frac{f_I D_A}{C_{I}^{eq}}\int_T \nabla(C_I C_A \nabla ( n)) \varphi_k d\Omega=$$ (157)

$$=-\frac{f_I D_A}{C_{I}^{eq}}\int_T \nabla(C_I C_A)\nabla \varphi_k d\Omega-\frac{f_I D_A}{C_{I}^{eq}}\int_T C_I C_A \nabla ( n) \nabla \varphi_k d\Omega.$$ (158)

The nonlinear terms of (3.95), $C_I(\mathbf{x}) C_A(\mathbf{x})$, $\nabla(C_I(\mathbf{x}) C_A(\mathbf{x}))$, and ln$n$ are expressed by the basis nodal functions $\varphi_1,\varphi_2,\varphi_3,\varphi_4$ on the following way,

$$C_I(\mathbf{x}) C_A(\mathbf{x}) = \sum_{i,j=1}^4 C_{I,i}C_{A,j}\varphi_i(\mathbf{x}) \varphi_j(\mathbf{x}),$$

$$\nabla(C_I(\mathbf{x}) C_A(\mathbf{x})) = \sum_{i,j=1}^4 C_{I,i}C_{A,j}\nabla(\varphi_i(\mathbf{x}) \varphi_j(\mathbf{x})),$$ (159)

$$\nabla( n)=\frac{1}{\sqrt{4n_{i}^{2}+C_{A,m}^{2}}} \sum_{i=1}^4 C_{A,i}\nabla\varphi_i(\mathbf{x}).$$

In the last expression we used the relationship (3.13) for the electron concentration $n$ in the case of single charged dopant. We substitute the relationships (3.96) into (3.95) and obtain,

$$\int_T \nabla\cdot\mathbf{J}_{AI} \varphi_k d\Omega=-\frac{f_I ... ...^4C_{I,i}C_{A,j}\int_T \nabla (\varphi_i \varphi_j)\nabla \varphi_k d\Omega -$$

$$-\frac{f_I D_A}{C_{I}^{eq}\sqrt{4n_{i}^{2}+C_{A,m}^{2}}}\sum_{i,j... ...}C_{A,l}\int_T \varphi_i \varphi_j \nabla \varphi_l \nabla \varphi_k d\Omega.$$ (160)

To complete the formulation of the weak equation form we have to evaluate following two integrals,

$$\int_T \nabla (\varphi_i \varphi_j)\nabla \varphi_k d\Omega=det... ...}_i \varphi^{e}_j)\mathbf{\Lambda}\nabla_e \varphi^{e}_k d\zeta d\eta d\xi,$$ (161)

$$\int_T \varphi_i \varphi_j \nabla \varphi_l \nabla \varphi_k d\... ...a_e \varphi^{e}_l \mathbf{\Lambda}\nabla_e \varphi^{e}_k d\zeta d\eta d\xi,$$ (162)

here we transfer the calculation to the normed $(\xi,\zeta,\eta)$ coordinate system. Transformation of the gradient operators is done by multiplication with the matrix $\mathbf{\Lambda}=\mathbf{(J^{-1})^T}$.

$$\mathbf{\Lambda}\nabla_e (\varphi^{e}_i \varphi^{e}_j)\mathbf{\La... ...rtial\eta} \\ \frac{\partial(\varphi^{e}_k)}{\partial\zeta} \\ \end{bmatrix}=$$

$$=\Bigl(\Lambda_{0,0}\frac{\partial(\varphi^{e}_i \varphi^{e}_j)}{... ...\partial\eta}+\Lambda_{0,2}\frac{\partial(\varphi^{e}_k)}{\partial\zeta}\Bigr)+$$

$$+ \Bigl(\Lambda_{1,0}\frac{\partial(\varphi^{e}_i \varphi^{e}_j)}... ...tial\eta}\!+\!\Lambda_{1,2}\frac{\partial(\varphi^{e}_k)}{\partial\zeta}\Bigr)+$$

$$+ \Bigl(\Lambda_{2,0}\frac{\partial(\varphi^{e}_i \varphi^{e}_j)}... ...l\eta}+\Lambda_{2,2}\frac{\partial(\varphi^{e}_k)}{\partial\zeta}\Bigr),\quad$$ (163)

Furthermore we rearrange (3.100),

$$\frac{\partial(\varphi^{e}_i \varphi^{e}_j)}{\partial\xi}\frac{\p... ...bda_{0,0}\Lambda_{0,1}+\Lambda_{1,0}\Lambda_{1,1}+\Lambda_{2,0}\Lambda_{2,1}) +$$

$$+\frac{\partial(\varphi^{e}_i \varphi^{e}_j)}{\partial\xi}\frac{\... ...mbda_{0,1}\Lambda_{0,0}+\Lambda_{1,1}\Lambda_{1,0}+\Lambda_{2,1}\Lambda_{2,0})+$$

$$+ \frac{\partial(\varphi^{e}_i \varphi^{e}_j)}{\partial\eta}\frac... ...mbda_{0,1}\Lambda_{0,2}+\Lambda_{1,1}\Lambda_{1,2}+\Lambda_{2,1}\Lambda_{2,2})+$$

$$+\frac{\partial(\varphi^{e}_i \varphi^{e}_j)}{\partial\zeta}\frac... ...mbda_{0,2}\Lambda_{0,1}+\Lambda_{1,2}\Lambda_{1,1}+\Lambda_{2,2}\Lambda_{2,1})+$$

$$+\frac{\partial(\varphi^{e}_i\varphi^{e}_j)}{\partial\zeta}\frac{... ...^{e}_k)}{\partial\zeta}(\Lambda_{0,2}^{2}+\Lambda_{1,2}^{2}+\Lambda_{2,2}^{2}).$$ (164)

and introduce the matrix $\mathbf{K}_{\xi}$, $\mathbf{K}_{\eta}$, $\mathbf{K}_{\zeta}$.

$${K}_{\xi}^{i,j}=\int_{\xi=0}^{1}\int_{\eta=0}^{1-\xi}\int_{\zeta=... ...frac{\partial(\varphi^{e}_i \varphi ^{e}_j)}{\partial\xi} d\zeta d\eta d\xi,$$

$${K}_{\eta}^{i,j}=\int_{\xi=0}^{1}\int_{\eta=0}^{1-\xi}\int_{\zeta... ...rac{\partial(\varphi^{e}_i \varphi ^{e}_j)}{\partial\eta} d\zeta d\eta d\xi,$$ (165)

$${K}_{\zeta}^{i,j}=\int_{\xi=0}^{1}\int_{\eta=0}^{1-\xi}\int_{\zet... ...rac{\partial(\varphi^{e}_i \varphi^{e}_j)}{\partial\zeta} d\zeta d\eta d\xi.$$

The values of the integrals (3.98) and (3.99) can be expressed as functions of local node indexes,

$$\alpha_{T}(i,j,k)=\int_T \nabla (\varphi_i \varphi_j)\nabla \varp... ...i,j}( \varphi_{\xi}^{k}(\Lambda_{0,0}^{2}+\Lambda_{1,0}^{2}+\Lambda_{2,0}^{2})+$$

$$+\varphi_{\eta}^{k}(\Lambda_{0,0}\Lambda_{0,1}+\Lambda_{1,0}\Lamb... ...bda_{0,0}\Lambda_{0,2}+\Lambda_{1,0}\Lambda_{1,2}+\Lambda_{2,0}\Lambda_{2,2}))+$$

$$+ {K}_{\eta}^{i,j}(\varphi_{\xi}^{k}(\Lambda_{0,1}\Lambda_{0,0}+\... ...,0})+\varphi_{\eta}^{k}(\Lambda_{0,1}^{2}+\Lambda_{1,1}^{2}+\Lambda_{2,1}^{2})+$$

$$+ \varphi_{\zeta}^{k}(\Lambda_{0,1}\Lambda_{0,2}+\Lambda_{1,1}\Lambda_{1,2}+\Lambda_{2,1}\Lambda))+$$

$$+{K}_{\zeta}^{i,j}(\varphi_{\xi}^{k}(\Lambda_{0,2}\Lambda_{0,0}+\... ...mbda_{0,2}\Lambda_{0,1}+\Lambda_{1,2}\Lambda_{1,1}+\Lambda_{2,2}\Lambda_{2,1})+$$

$$+ \varphi_{\zeta}^{k}(\Lambda_{0,2}^{2}+\Lambda_{1,2}^{2}+\Lambda_{2,2}^{2}))),$$ (166)

$$\beta_{T}(i,j,l,k)=\int_T \varphi_i \varphi_j \nabla \varphi_l \n... ...^{l}( \varphi_{\xi}^{k}(\Lambda_{0,0}^{2}+\Lambda_{1,0}^{2}+\Lambda_{2,0}^{2})+$$

$$+\varphi_{\eta}^{k}(\Lambda_{0,0}\Lambda_{0,1}+\Lambda_{1,0}\Lamb... ...bda_{0,0}\Lambda_{0,2}+\Lambda_{1,0}\Lambda_{1,2}+\Lambda_{2,0}\Lambda_{2,2}))+$$

$$+ {\varphi}_{\eta}^{i}(\varphi_{\xi}^{k}(\Lambda_{0,1}\Lambda_{0,... ...,0})+\varphi_{\eta}^{k}(\Lambda_{0,1}^{2}+\Lambda_{1,1}^{2}+\Lambda_{2,1}^{2})+$$

$$+ \varphi_{\zeta}^{k}(\Lambda_{0,1}\Lambda_{0,2}+\Lambda_{1,1}\Lambda_{1,2}+\Lambda_{2,1}\Lambda))+$$

$$+{\varphi}_{\zeta}^{i}(\varphi_{\xi}^{k}(\Lambda_{0,2}\Lambda_{0,... ...mbda_{0,2}\Lambda_{0,1}+\Lambda_{1,2}\Lambda_{1,1}+\Lambda_{2,2}\Lambda_{2,1})+$$

$$+ \varphi_{\zeta}^{k}(\Lambda_{0,2}^{2}+\Lambda_{1,2}^{2}+\Lambda_{2,2}^{2}))).$$ (167)

In (3.103) and (3.104), for the sake of simplicity, we used the abbervations, $\varphi_{X}^k=\frac{\partial\varphi^{e}_k}{\partial X}$ for $X=\xi,\eta,\zeta$. It can easily been shown that the functions $\alpha_T(i,j,k)$ and $\beta_T(i,j,l,k)$ satisfy following relations,

$$\sum_{k=1}^4\alpha_T(i,j,k)=0, \sum_{i=1}^4\alpha_T(i,j,k)=K_{jk}(T), \sum_{j=1}^4\alpha_T(i,j,k)=K_{ik}(T),$$

$$\sum_{k=1}^4\beta_T(i,j,l,k)=0, \sum_{l=1}^4\beta_T(i,j,l,k)=0, \sum_{i,j=1}^4\beta_T(i,j,l,k)=K_{lk}(T).$$ (168)

We introduce now the backward Euler time discretization scheme for (3.91) and by applying the functions $\alpha_{T}(i,j,k)$ and $\beta_{T}(i,j,l,k)$ we write the $\ell_{k,1}$ functions (for $k=1,2,3,4$),

$$\ell_{k,1}^e(\mathbf{c}_A^n,\mathbf{c}_I^n,\mathbf{c}_V^n)=\frac{... ...I\frac{C_{I,i}^n}{C_{I}^{eq}}+(1-f_I)\frac{C_{V,i}^n}{C_{V}^{eq}}\Bigr)\Bigr\}+$$ (169)

$$+\frac{D_A }{\sqrt{4n_{i}^{2}+(C_{A,m}^n)^{2}}}\sum_{i,j,l=1}^4\B... ...I\frac{C_{I,i}^n}{C_{I}^{eq}}+(1-f_I)\frac{C_{V,i}^n}{C_{V}^{eq}}\Bigr)\Bigr\},$$ (170)

and residuum vector $R_{k,1}^e$ as

$$R_{k,1}^e(\mathbf{c}_A^n,\mathbf{c}_I^n,\mathbf{c}_V^n)=\ell_{k,1... ...thbf{c}_I^n,\mathbf{c}_V^n)-\frac{1}{\Delta t}\sum_{i=1}^4 M_{ik}C_{A,i}^{n-1}.$$ (171)

Equations (3.74) and (3.75) of the Mulvaney-Richardson model have the same structure, it is sufficient, if we consider in more detail only (3.74),

$$\frac{\partial C_I}{\partial t}=D_I\Delta C_I - k_{f} ( C_I C_V - C_{I}^{eq} C_{V}^{eq}) +\nabla \cdot \mathbf{J}_{AI}.$$ (172)

We need to consider only $D_I\Delta C_I - k_{f} ( C_I C_V - C_{I}^{eq} C_{V}^{eq})$, the rest of the equation has already been treated above. Again, we use linear basis nodal functions $\varphi_k$ and apply Green's theorem,

$$-D_I\sum_{i,=1}^{4}C_{I,i}\int_T \nabla \varphi_i \nabla Nk d\Om... ... \varphi_j \varphi_kd\Omega + k_fC_{I}^{eq} C_{V}^{eq}\int_T \varphi_k d\Omega.$$ (173)

we set $\gamma_T(i,j,k)=\int_T \varphi_i \varphi_j \varphi_k d\Omega$, and proceed,

$$-D_I\sum_{i=1}^{4}C_{I,i} K_{i,k} - k_f\sum_{i,j=1}^{4}C_{I,i}C_{V,j}\gamma_T(i,j,k) + k_f C_{I}^{eq} C_{V}^{eq}V_k.$$ (174)

Finally, we have the discretizied form of equation (3.74),

$$\ell_{k,2}^e(\mathbf{c}_A^n,\mathbf{c}_I^n,\mathbf{c}_V^n)=\sum_{... ...j=1}^{4}C^{n}_{I,i}C_{V,j}^{n}\gamma_T(i,j,k) - C_{I}^{eq} C_{V}^{eq}V_k\Bigr)+$$

$$+\frac{f_I D_A}{C_{I}^{eq}}\sum_{i,j=1}^{4}\Bigl\{C^{n}_{A,i}\Big... ...{2}+(C^n_{A,m})^2}}\sum_{l=1}^{4}C_{I,l}^{n}\beta_T(i,j,l,k)\Bigr)\Bigr\},\quad$$ (175)

and the residuum vector $R_{k,2}^e$ as,

$$R_{k,2}^e(\mathbf{c}_A^n,\mathbf{c}_I^n,\mathbf{c}_V^n)=\ell_{k,2... ...thbf{c}_I^n,\mathbf{c}_V^n)-\frac{1}{\Delta t}\sum_{i=1}^4 M_{ik}C_{I,i}^{n-1}.$$ (176)

and analogously for the vacancy balance equation,

$$\ell_{k,3}^e(\mathbf{c}_A^n,\mathbf{c}_I^n,\mathbf{c}_V^n)=\sum_{... ...j=1}^{4}C^{n}_{I,i}C_{V,j}^{n}\gamma_T(i,j,k) - C_{I}^{eq} C_{V}^{eq}V_k\Bigr)+$$

$$+\frac{(1-f_I) D_A}{C_{V}^{eq}}\sum_{i,j=1}^{4}\Bigl\{C^{n}_{A,i}... ...{2}+(C^n_{A,m})^2}}\sum_{l=1}^{4}C_{V,l}^{n}\beta_T(i,j,l,k)\Bigr)\Bigr\},\quad$$ (177)

and the corresponding residuum vector $R_{k,3}^e$ as,

$$R_{k,3}^e(\mathbf{c}_A^n,\mathbf{c}_I^n,\mathbf{c}_V^n)=\ell_{k,3... ...thbf{c}_I^n,\mathbf{c}_V^n)-\frac{1}{\Delta t}\sum_{i=1}^4 M_{ik}C_{V,i}^{n-1}.$$ (178)

Using the functions $\ell^e_{k,1}$, $\ell^e_{k,2}$, and $\ell^e_{k,3}$ for $k=1,2,3,4$ we construct the nucleus matrix $\mathbf{\Pi}(T)$ of the diffusion problem,

$$\mathbf{\Pi}(T)=\frac{\partial (\ell_{1,1}^e,\ell_{1,2}^e,\ell_{1... ...3,3}^e,\ell_{3,4}^e)}{\partial (\mathbf{c}_A^n,\mathbf{c}_I^n,\mathbf{c}_V^n)}.$$ (179)

Note that the global matrix $\mathbf{G}$ in this case has rank $3\cdot N$.