Computation of Torques
in Magnetic Tunnel Junctions

Chapter 6 Spin and Charge Drift-Diffusion for the Computation of Spin Torques and Finite Element Implementation

The driving factor of the writing process in modern MRAM cells is the torque generated by the polarization process of electrons transiting through the ferromagnetic layers in an MTJ stack. The simulation of the switching process can be achieved by solving the LLG equation for magnetization dynamics, with the inclusion of a term describing the torque acting on the magnetization. As seen in the previous chapters, such term can be defined by employing a Slonczewski-like torque approach [16]. This, however, allows to approximately simulate the magnetization dynamics of a thin free layer only. A more complete description of the process can be obtained by computing the non-equilibrium spin accumulation $\mathbf{S}$ across the whole structure. The following sections demonstrate that a more general expression for the torque term ${\mathbf{T}}_{\text{S}}$ entering (3.44) can be given as [25], [119], [120]

$$\begin{array}{}\text{(6.1)}& {\mathbf{T}}_{\text{S}}=-D_{\text{e}}\frac{\mathbf{m}\times \mathbf{S}}{{\lambda }_J^2}-D_{\text{e}}\frac{\mathbf{m}\times (\mathbf{m}\times \mathbf{S})}{{\lambda }_{\phi }^2},\end{array}$$

where $D_{\text{e}}$ is the electron diffusion coefficient. The first term describes precession around the exchange field and is characterized by the exchange length ${\lambda }_J$, and the second term describes the dephasing process of the spins of the transiting electrons, and is characterized by the dephasing length ${\lambda }_{\phi }$.

The spin accumulation $\mathbf{S}$ describes the deviation of the polarization of the conducting electrons from the equilibrium configuration created by a charge current density ${\mathbf{J}}_{\text{C}}$, in units of the transported magnetic moment (A/m). Thus, by definition $\mathbf{S}$ is non-zero only when an electric current is flowing through the system [83]. A solution for $\mathbf{S}$ in all nonmagnetic and ferromagnetic layers of an MRAM cell can be obtained by means of the spin and charge drift-diffusion formalism.

6.1 Spin and Charge Drift-Diffusion

Spin and charge drift-diffusion equations in the presence of a spin-valve structure where the ferromagnetic layers have co-linear (parallel or anti-parallel) magnetization were first derived by T. Valet and A. Fert [121]. A formulation of the equations for the computation of $\mathbf{S}$ in the presence of arbitrary magnetization orientation was later reported by S. Zhang, P. Levy and A. Fert [98], and will be hereby summarized. In a magnetic multilayer with the current perpendicular to the plane of the layer, along the x-direction, the linear response of the current to the electric field can be written in spinor form as

$$\begin{array}{}\text{(6.2)}& \hat{j}=\hat{C}E(x)-\hat{D}\frac{\partial \hat{n}}{\partial x},\end{array}$$

where $E$ is the electric field, which can be expressed in terms of the electric potential $V$ as $E=-\partial V/\partial x$, and $\hat{j}$, $\hat{C}$, $\hat{D}$, and $\hat{n}$ are the $2\times 2$ spinor matrices representing the current density, the conductivity, the diffusion coefficient, and the accumulation at a given position. For a metal, the diffusion constant and the conductivity are connected via the Einstein relation $\hat{C}=e^2{\hat{N}}_{\text{F}}\hat{D}$, with ${\hat{N}}_{\text{F}}$ the spin-dependent density of states at the Fermi level. These matrices can in general be expressed in terms of the Pauli matrix vector $\mathit{\sigma }$:

$$\begin{array}{}\text{(6.3)}& \hat{C}& =C_0\hat{I}+\mathbf{C}\mathbf{\cdot }\hat{\mathit{\sigma }}\text{(6.4)}& \hat{D}& =D_0\hat{I}+\mathbf{D}\mathbf{\cdot }\hat{\mathit{\sigma }}\text{(6.5)}& \hat{n}& =n_0\hat{I}+\mathbf{S}\mathbf{\cdot }\hat{\mathit{\sigma }}\end{array}$$ $2n_0$ is the charge accumulation, $C_0$ is related to the conductivity $\sigma$ by $\sigma =2C_0$, and $D_0$ is related to the electron diffusion coefficient by $D_{\text{e}}=2D_0$. The two vectors $\mathbf{C}$ and $\mathbf{D}$ are related to the different conductivity and diffusion constant for the majority and minority electrons, labeled ${\sigma }^{\uparrow }$, $D_{\text{e}}^{\uparrow }$ and ${\sigma }^{\downarrow }$, $D_{\text{e}}^{\downarrow }$, respectively. This difference is characterized by the conductivity polarization parameter ${\beta }_{\sigma }=({\sigma }^{\uparrow }-{\sigma }^{\downarrow })/({\sigma }^{\uparrow }+{\sigma }^{\downarrow })$ and the diffusion polarization parameter ${\beta }_{\text{D}}=(D_{\text{e}}^{\uparrow }-D_{\text{e}}^{\downarrow })/(D_{\text{e}}^{\uparrow }+D_{\text{e}}^{\downarrow })$. For a transition ferromagnet, $\mathbf{C}$ is then defined as $\mathbf{C}={\beta }_{\sigma }{C}_0\mathbf{m}$ and $\mathbf{D}={\beta }_{\text{D}}{D}_0\mathbf{m}$. The two polarization parameters are different when the density of states for majority and minority electrons is not the same. By inserting (6.3)-(6.5) into (6.2), one obtains the following expressions for the charge current density $J_{\text{C}}$ and spin polarization current density ${\mathbf{J}}_{\text{S}}$ along the x-direction:

$$\begin{array}{}\text{(6.6a)}& J_{\text{C}}\equiv \text{Re}\left(\text{Tr}\left[\hat{j}\right]\right)& =\sigma E+\frac{e}{{\mu }_{\text{B}}}{\beta }_{\text{D}}{D}_{\text{e}}\frac{\partial \mathbf{S}}{\partial x}\mathbf{\cdot }\mathbf{m}\text{(6.6b)}& {\mathbf{J}}_{\text{S}}\equiv \text{Re}\left(\text{Tr}\left[\hat{\mathit{\sigma }}\hat{j}\right]\right)& =-\frac{{\mu }_{\text{B}}}{e}{\beta }_{\sigma }\mathbf{m}\sigma E-D_{\text{e}}\frac{\partial \mathbf{S}}{\partial x}\end{array}$$ Terms involving the spatial derivative of the charge accumulation have been ignored, as in metals the distribution of charges can be considered uniform. The factor ${\mu }_{\text{B}}/e$ converts from the units of electric charge, $A/m^2$, to the ones of the spin polarization current density, $A/s$. The latter can be converted to the typical units of spin current density by multiplication with the factor $\hslash /(2{\mu }_{\text{B}})$. The spin polarization current density will be referred to as spin current for brevity from now on.

Equation (6.6) only describes the spin and charge currents on the axis perpendicular to a multilayer structure. The generalization to a three-dimensional formulation leads to the following expressions:

$$\begin{array}{}\text{(6.7a)}& {\mathbf{J}}_{\text{C}}& =\sigma \mathbf{E}+\frac{e}{{\mu }_{\text{B}}}{\beta }_{\text{D}}{D}_{\text{e}}{(\mathbf{\nabla }\mathbf{S})}^{\text{T}}\mathbf{m}\text{(6.7b)}& {\tilde{\mathbf{J}}}_{\text{S}}& =-\frac{{\mu }_{\text{B}}}{e}{\beta }_{\sigma }\mathbf{m}\otimes \left(\sigma \mathbf{E}\right)-D_{\text{e}}\mathbf{\nabla }\mathbf{S}\end{array}$$ $\otimes$ is the outer product, ${\mathbf{J}}_{\mathbf{C}}$ is the charge current density vector, ${\tilde{\mathbf{J}}}_{\text{S}}$ is the spin current tensor, where the components $J_{\text{S,ij}}$ indicate the flow of the i-th component of spin polarization in the j-th direction, and $\mathbf{\nabla }\mathbf{S}$ is the vector gradient of $\mathbf{S}$, with components ${(\mathbf{\nabla }\mathbf{S})}_{\text{ij}}=\partial S_{\text{i}}/\partial x_{\text{j}}$. The term ${(\mathbf{\nabla }\mathbf{S})}^{\text{T}}\mathbf{m}$ is a vector with components ${\left({(\mathbf{\nabla }\mathbf{S})}^{\text{T}}\mathbf{m}\right)}_{\text{i}}=\sum _{\text{j}}\left(\partial S_{\text{j}}/\partial x_{\text{i}}\right)m_j$. In case the current density distribution is known, the spin current can be expressed in terms of the charge current by inserting (6.7a) into (6.7b):

$$\begin{array}{}\text{(6.8)}& {\tilde{\mathbf{J}}}_{\text{S}}=-\frac{{\mu }_B}{e}{\beta }_{\sigma }\mathbf{m}\otimes ({\mathbf{J}}_{\text{C}}-\frac{e}{{\mu }_B}{\beta }_D{D}_e{(\mathbf{\nabla }\mathbf{S})}^{\text{T}}\mathbf{m})-D_e\mathbf{\nabla }\mathbf{S}\end{array}$$

The equation of motion for the spin accumulation in a ferromagnet is determined by the interaction of the local magnetization moment and the spin of the itinerant electrons. This interaction can be described by the Hamiltonian term

$$\begin{array}{}\text{(6.9)}& H_{\text{int}}=-J\mathbf{S}\mathbf{\cdot }\mathbf{m},\end{array}$$

where $J$ is the coupling strength between local moments and itinerant electrons. The following equation of motion for the spin accumulation can thus be derived:

$$\begin{array}{}\text{(6.10)}& \frac{\mathrm{d}\mathbf{S}}{\mathrm{d}t}-\frac{J}{\hslash }\mathbf{m}\times \mathbf{S}=-\frac{\mathbf{S}}{{\tau }_{\text{sf}}}\end{array}$$

${\tau }_{\text{sf}}$ is the spin-flip relaxation time of the conduction electrons. The second term on the left hand side describes the processional motion of the accumulation due to the exchange interaction. It plays a role when the non-equilibrium spin accumulation and the local magnetization moment are not parallel. The rate of change of $\mathbf{S}$ can be expressed in terms of its partial derivative in time and the divergence of the spin current as $\mathrm{d}\mathbf{S}/\mathrm{d}t=\partial \mathbf{S}/\partial t+\mathbf{\nabla }\mathbf{\cdot }{\tilde{\mathbf{J}}}_{\text{S}}$. By replacing this expression in (6.10), the equation for the spin accumulation becomes

$$\begin{array}{}\text{(6.11)}& \frac{\partial \mathbf{S}}{\partial t}=-\mathbf{\nabla }\mathbf{\cdot }{\tilde{\mathbf{J}}}_{\text{S}}-D_{\text{e}}\frac{\mathbf{S}}{{\lambda }_{sf}^2}+D_{\text{e}}\frac{\mathbf{m}\times \mathbf{S}}{{\lambda }_J^2},\end{array}$$

where the two length parameters ${\lambda }_{\text{sf}}$ and ${\lambda }_{\text{J}}$ are defined as ${\lambda }_{\text{sf}}=\sqrt{D_{\text{e}}{\tau }_{\text{sf}}}$ and ${\lambda }_{\text{J}}=\sqrt{D_{\text{e}}\hslash /J}$. $\mathbf{\nabla }\mathbf{\cdot }{\tilde{\mathbf{J}}}_{\text{S}}$ is the divergence of ${\tilde{\mathbf{J}}}_{\text{S}}$ with components ${(\mathbf{\nabla }\mathbf{\cdot }{\tilde{\mathbf{J}}}_{\text{S}})}_{\text{i}}=\sum _{\text{j}}\partial J_{\text{S,ij}}/\partial x_{\text{j}}$.

Equation (6.11) can be further simplified by taking into consideration the fact that the typical time-scales for the spin accumulation and the magnetization are very different. While the former is of the order of picoseconds, the latter is of the order of nanoseconds [98]. For the computation of the spin torque, to be added to the LLG equation (3.44), it is thus sufficient to consider a steady-state expression for the spin accumulation. This assumption was numerically verified in [122]. With $\partial \mathbf{S}/\partial t=0$, the equation describing the spin accumulation becomes

$$\begin{array}{}\text{(6.12)}& -\mathbf{\nabla }\mathbf{\cdot }{\tilde{\mathbf{J}}}_{\text{S}}-D_{\text{e}}\frac{\mathbf{S}}{{\lambda }_{sf}^2}+D_{\text{e}}\frac{\mathbf{m}\times \mathbf{S}}{{\lambda }_J^2}=0.\end{array}$$

An expression for the torque can be extracted from (6.12) by considering a conservation argument [123]. In the original work of Slonczewski [14], it was argued that the spin current lost by the conducting electrons, described by the divergence term, must have been gained by the magnetization, due to conservation of the total angular momentum. This is always valid in the absence of additional relaxation mechanisms for the spin accumulation [99]. However, in the presence of spin-orbit coupling, part of the angular momentum is transferred to the lattice. This process is responsible for the spin-flip relaxation, described by the ${\lambda }_{\text{sf}}$ term. Spin current lost due to this mechanism does not contribute to the torque acting on the magnetization. The torque can then be obtained as

$$\begin{array}{}\text{(6.13)}& {\mathbf{T}}_{\text{S}}=-\mathbf{\nabla }\mathbf{\cdot }{\tilde{\mathbf{J}}}_{\text{S}}-D_{\text{e}}\frac{\mathbf{S}}{{\lambda }_{\text{sf}}^2}=-D_{\text{e}}\frac{\mathbf{m}\times \mathbf{S}}{{\lambda }_J^2}.\end{array}$$

6.1.1 Spin Dephasing Term

The previous section showed the derivation of a spin torque expression where both the absorption and decay of the transverse spin accumulation components are governed by the exchange length ${\lambda }_{\text{J}}$. Another possible mechanism for the absorption of the transverse components is the dephasing process. Dephasing occurs when, after traveling for a certain distance, different spins have precessed different amounts, so that their transverse components tend to cancel out [120]. In the context of spin transfer torques, this can occur from the variation of electron velocities over the Fermi surface [124] or from spins precessing at the same rate but arriving at different times due to scattering [125]. The behavior of the spin accumulation with both precessional and dephasing terms was described in terms of the Continuous Random Matrix Theory (CRMT) in [125], and the equivalence of CMRT and the spin and charge drift-diffusion formalism was shown. This section gives a brief overview of the derivation thereby presented.

The Landauer-Büttiker formalism, or scattering formalism, is the base for CRMT, and is expressed through a scattering matrix $S$ relating the initial state and the final state of a physical system undergoing a scattering process [126]. The form of the $S$ matrix is derived by considering a system where two leads are connected to a central scattering region. The leads are quasi-one-dimensional waveguides. By assuming the solutions to the Schrödinger equation ${\mathrm{\Psi }}_{\text{n}}^{\text{±}}=\psi (y)\exp (\text{±}i{k}_{\text{n}}x)$ are known outside of the scattering region, the wave functions ${\mathrm{\Psi }}^{\text{±}}(x,y{)}_{\text{L/R}}$ on the left (L) and right (R) side can be expressed as a linear combination of them:

$$\begin{array}{}\text{(6.14)}& {\mathrm{\Psi }}^{\text{±}}(x,y{)}_{\text{L/R}}=\sum _{\text{n}}{a}_{\text{n,L/R}}^{\text{±}}\psi (y)\exp (\text{±}i{k}_{\text{n}}x)\end{array}$$

${\mathrm{\Psi }}^{+}(x,y{)}_{\text{L}}$ and ${\mathrm{\Psi }}^{-}(x,y{)}_{\text{R}}$ are the incoming wave functions, while ${\mathrm{\Psi }}^{-}(x,y{)}_{\text{L}}$ and ${\mathrm{\Psi }}^{+}(x,y{)}_{\text{R}}$ are the outgoing wave functions. By expressing the vectors containing the expansion coefficients $a_{\text{n,L/R}}^{\text{±}}$ as ${\mathrm{\Psi }}_{\text{L/R}}^{\text{±}}$, the incoming and outgoing vectors can be related by a matrix $S$ composed of reflection ($r,r^{\prime }$) and transmission ($t,t^{\prime }$) submatrices:

$$\begin{array}{}\text{(6.15)}& \left(\begin{array}{c}{\mathrm{\Psi }}_{\text{L}}^{-}\\ {\mathrm{\Psi }}_{\text{R}}^{+}\end{array}\right)=\left(\begin{array}{cc}{r}^{\prime }& t\\ t^{\prime }& r\end{array}\right)\left(\begin{array}{c}{\mathrm{\Psi }}_{\text{L}}^{+}\\ {\mathrm{\Psi }}_{\text{R}}^{-}\end{array}\right),\left(\begin{array}{cc}{r}^{\prime }& t\\ t^{\prime }& r\end{array}\right)=S\end{array}$$ If two conductors A and B described by this scattering formalism are in series, the transmission and reflection of the whole system can be written as

$$\begin{array}{}\text{(6.16a)}& t_{\text{AB}}=t_{\text{B}}\frac{1}{1-r_{\text{A}}{r}_{\text{B}}^{\prime }}{t}_{\text{A}},\text{(6.16b)}& r_{\text{AB}}=r_{\text{B}}+t_{\text{B}}\frac{1}{1-r_{\text{A}}{r}_{\text{B}}^{\prime }}{r}_{\text{A}}{t}_{\text{B}}^{\prime }.\end{array}$$ The conductance of a system can be expressed in terms of its transmission matrix $t$ as

$$\begin{array}{}\text{(6.17)}& G=\frac{e^2}{h}\text{T},\text{T}=\text{Tr}\left[t{t}^{†}\right].\end{array}$$

The value $\text{T}$ corresponds to the sum over the probabilities of being transmitted from one side of the system to the other, and $h$ is the Planck constant. The same expression holds true for the reflection probabilities, represented by the quantity $\text{R}=\text{Tr}\left[r{r}^{†}\right]$. By taking a semi-classical limit and assuming that successive reflections mix together different channels, it can be assumed that each electron traveling through the system picks up a certain phase. The transmission through the whole system takes then the form

$$\begin{array}{}\text{(6.18)}& t_{\text{AB}}=t_{\text{B}}\frac{1}{1-\left|r_{\text{A}}{r}_{\text{B}}^{\prime }\right|z}{t}_{\text{A}},\end{array}$$

where $z=\exp (i\phi )$ describes the phase acquired by the traveling electrons. An expression for the conductance is obtained by taking the average over all values of $\phi$:

$$\begin{array}{}\text{(6.19)}& ⟨G_{\text{AB}}⟩=\frac{e^2}{h}⟨{\text{T}}_{\text{AB}}⟩\end{array}$$

The average value can be computed through the residual theorem by integrating over the phase in the complex space. This leads to the expression

$$\begin{array}{}\text{(6.20)}& ⟨{\text{T}}_{\text{AB}}⟩={\text{T}}_{\text{B}}\frac{1}{1-{\text{R}}_{\text{A}}{\text{R}}_{\text{B}}^{\prime }}{\text{T}}_{\text{A}},\end{array}$$

which has the same form of (6.16a), but it relates probabilities instead of amplitudes. The same is also true for the quantity $⟨{\text{R}}_{\text{AB}}⟩$. This means that, within the semi-classical limit, transport can be described by a scattering-like formalism concerning only probabilities.

The same theoretical description, employed until now only in the presence of multiple channels, can also be applied in the presence of the additional spin degree of freedom. In this case, the transmission and reflection matrices have dimensions of $4N_{\text{ch}}\times 4N_{\text{ch}}$, with $N_{\text{ch}}$ the number of channels. The random matrix theory (RMT) employed in [126] can be applied to account for both mixing of the contributions in different channels and loss of coherence of the waves traveling through the system, in a similar fashion to the average over the phase previously shown. The resulting theory presents a structure similar to the scattering formalism, with $4\times 4$ spin-dependent transmission and reflection matrices, indicated by a hat. The block matrix $\hat{S}$ describing transport through a conductor takes the form

$$\begin{array}{}\text{(6.21)}& \hat{S}=\left(\begin{array}{cc}{\hat{r}}^{\prime }& \hat{t}\\ {\hat{t}}^{\prime }& \hat{r}\end{array}\right).\end{array}$$

Each of the four submatrices is defined from the $4N_{\text{ch}}\times 4N_{\text{ch}}$ matrices of the scattering theory, relating different spins and channels, by the following expression:

$$\begin{array}{}\text{(6.22)}& {\hat{t}}_{\sigma \eta ,{\sigma }^{\prime }{\eta }^{\prime }}=\frac{1}{N_{\text{ch}}}{\text{Tr}}_{N_{\text{ch}}}\left[t_{\sigma \eta }{t}_{{\sigma }^{\prime }{\eta }^{\prime }}^{†}\right]\end{array}$$

$\sigma$ and $\eta$ indicate the up ($\uparrow$) or down ($\downarrow$) spin orientation. Similar expressions hold for the other transmission and reflection matrices. The most important contributions to the obtained matrix are the probabilities ${\text{T}}_{\sigma \eta }={\hat{t}}_{\sigma \eta ,\sigma \eta }$ and the mixing coefficients ${\text{T}}_{\text{mx}}={\hat{t}}_{\sigma \sigma ,-\sigma -\sigma }$. The other entries in the matrix have a negligible contribution to the overall transport description and can be disregarded [126]. The final form of the $\hat{t}$ matrix is then given by

$$\begin{array}{}\text{(6.23)}& \hat{t}=\left(\begin{array}{cccc}{\text{T}}_{\uparrow \uparrow }& 0& 0& {\text{T}}_{\uparrow \downarrow }\\ 0& {\text{T}}_{\text{mx}}& 0& 0\\ 0& 0& {\text{T}}_{\text{mx}}^{\ast }& 0\\ {\text{T}}_{\downarrow \uparrow }& 0& 0& {\text{T}}_{\downarrow \downarrow }\end{array}\right).\end{array}$$

Similar expressions hold true for the other matrices entering (6.21). The coefficients ${\text{T}}_{\sigma \eta }={\hat{t}}_{\sigma \eta ,\sigma \eta }$ describe the probability of an electron with spin $\sigma$ to be transmitted with spin $\eta$. The mixing coefficients are related to the precession and loss of coherence of electrons with a spin non-colinear with the magnetization:

$$\begin{array}{}\text{(6.24)}& {\text{T}}_{\text{mx}}\approx \exp (-d/l_{\perp }+id/l_{\text{L}})\end{array}$$

$d$ is the thickness of the ferromagnetic conductor, and $l_{\text{L}}$ is the Larmor precession length. The term containing $l_{\perp }$ is related to the fact that the amount of precession depends on the direction of incidence of the electron into the ferromagnet. Two electrons with different incidence angles will undergo a dephasing process. This results is an exponentially decaying spin density orthogonal to the local magnetization which absorbs the transverse components of the spin current. $l_{\perp }$ indicates the typical length over which these absorption process occurs.

Most of the quantities described until now can be recast in terms of semi-classical properties. The conductance can now be expressed as

$$\begin{array}{}\text{(6.25)}& G=({\text{T}}_{\uparrow \uparrow }+{\text{T}}_{\uparrow \downarrow }+{\text{T}}_{\downarrow \uparrow }+{\text{T}}_{\downarrow \downarrow })/R_{\text{sh}},\end{array}$$

where $R_{\text{sh}}=h/\left(e^2{N}_{\text{ch}}\right)$ is the Sharvin resistance. It corresponds to the resistance of a perfectly transparent waveguide, with $N_{\text{ch}}$ propagating transversal modes. As shown in [126], the addition laws for transmission and reflection in equation (6.20) are still valid. Moreover, it is shown that the number of degrees of freedom can be reduced from the incoming and outgoing wave function coefficients to vectors of length 4, storing only the charge and spin degrees of freedom. These vectors are noted $P_{\text{4,L/R}}^{\text{±}}$, and the incoming and outgoing vectors are related through the $\hat{S}$ matrix (6.21) as

$$\begin{array}{}\text{(6.26)}& \left(\begin{array}{c}{P}_{\text{4,L}}^{-}\\ P_{\text{4,R}}^{+}\end{array}\right)=\hat{S}\left(\begin{array}{c}{P}_{\text{4,L}}^{+}\\ P_{\text{4,R}}^{-}\end{array}\right).\end{array}$$

As shown in [125], one can demonstrate the equivalence of the random matrix description to the spin and charge drift-diffusion formalism. First, transport in a discretized one-dimensional conductor, divided into $N$ parts, is considered. The $P_{\text{4}}$ vectors of each part can be related through the node equations

$$\begin{array}{}\text{(6.27a)}& P_{\text{4,L}}^{+}[i+1]& =P_{\text{4,R}}^{+}\left[i\right],\text{(6.27b)}& P_{\text{4,L}}^{-}[i+1]& =P_{\text{4,R}}^{-}\left[i\right],\end{array}$$ where $i$ is an index going from $1$ to $N$. Through a variable change, the $P_{\text{4}}$ vectors can be employed to obtain the four dimensional potential ${\mu }_{\text{4}}$ and current $j_{\text{4}}$ in the system [127]:

$$\begin{array}{}\text{(6.28a)}& {\mu }_{\text{4}}=\frac{1}{2}(P_{\text{4,L}}^{+}+P_{\text{4,L}}^{-})\text{(6.28b)}& j_{\text{4}}=\frac{1}{2}(P_{\text{4,L}}^{+}-P_{\text{4,L}}^{-})\end{array}$$ After the change of variables, and by assuming that transmission and reflection properties do not depend on the direction, equations (6.26) and (6.27) can be combined to obtain the following relations:

$$\begin{array}{}\text{(6.29a)}& 2({\mu }_{\text{4}}[i+1]-{\mu }_{\text{4}}\left[i\right])& =(\hat{r}-\hat{t}+1)({\mu }_{\text{4}}[i+1]-{\mu }_{\text{4}}\left[i\right]-j_{\text{4}}[i+1]-j_{\text{4}}\left[i\right])\text{(6.29b)}& 2(j_{\text{4}}[i+1]-j_{\text{4}}\left[i\right])& =(\hat{r}+\hat{t}-1)({\mu }_{\text{4}}[i+1]+{\mu }_{\text{4}}\left[i\right]-j_{\text{4}}[i+1]+j_{\text{4}}\left[i\right])\end{array}$$

In order to go from the discrete description presented until now to a continuous one which recovers the drift-diffusion formalism, one must go to the limit where the length of each material slice goes to the infinitesimal quantity $\mathrm{d}x$. This gives the possibility of describing the thin slices’ $\hat{S}$ matrix components with the following parametrization [128]:

$$\begin{array}{}\text{(6.30a)}& \hat{t}& =1-{\mathrm{\Lambda }}_{\text{t}}\mathrm{d}x\text{(6.30b)}& \hat{r}& ={\mathrm{\Lambda }}_{\text{r}}\mathrm{d}x\end{array}$$ ${\mathrm{\Lambda }}_{\text{t}}$ and ${\mathrm{\Lambda }}_{\text{r}}$ are $4\times 4$ matrices. These expressions can be inserted back into equation (6.29), to obtain the following formulation:

$$\begin{array}{}\text{(6.31a)}& 2({\mu }_{\text{4}}(x+\mathrm{d}x)-{\mu }_{\text{4}}\left(x\right))& =\mathrm{d}x({\mathrm{\Lambda }}_{\text{r}}+{\mathrm{\Lambda }}_{\text{t}})({\mu }_{\text{4}}(x+\mathrm{d}x)-{\mu }_{\text{4}}\left(x\right)-j_{\text{4}}(x+\mathrm{d}x)-j_{\text{4}}\left(x\right))\text{(6.31b)}& 2(j_{\text{4}}(x+\mathrm{d}x)-j_{\text{4}}\left(x\right))& =\mathrm{d}x({\mathrm{\Lambda }}_{\text{r}}-{\mathrm{\Lambda }}_{\text{t}})({\mu }_{\text{4}}(x+\mathrm{d}x)+{\mu }_{\text{4}}\left(x\right)-j_{\text{4}}(x+\mathrm{d}x)+j_{\text{4}}\left(x\right))\end{array}$$ The equations can then be developed to the first order in $\mathrm{d}x$ to obtain the continuous differential equations

$$\begin{array}{}\text{(6.32a)}& -\frac{\partial {\mu }_{\text{4}}}{\partial x}\left(x\right)& =2{\mathrm{\Lambda }}_{\mathrm{\Omega }}{j}_{\text{4}}(x),\text{(6.32b)}& -\frac{\partial j_{\text{4}}}{\partial x}\left(x\right)& =2{\mathrm{\Lambda }}_{\mathrm{\Xi }}{\mu }_{\text{4}}(x),\end{array}$$ where $2{\mathrm{\Lambda }}_{\mathrm{\Omega }}={\mathrm{\Lambda }}_{\text{t}}+{\mathrm{\Lambda }}_{\text{r}}$ and $2{\mathrm{\Lambda }}_{\mathrm{\Xi }}={\mathrm{\Lambda }}_{\text{t}}-{\mathrm{\Lambda }}_{\text{r}}$. These matrices, defined so that they fit the Valet-Fert theory, take the form

$$\begin{array}{}\text{(6.33a)}& {\mathrm{\Lambda }}_{\mathrm{\Omega }}& =\left(\begin{array}{cccc}{\mathrm{\Gamma }}_{\uparrow }& 0& 0& 0\\ 0& {\mathrm{\Gamma }}_{\text{mx,eff}}& 0& 0\\ 0& 0& {\mathrm{\Gamma }}_{\text{mx,eff}}^{\ast }& 0\\ 0& 0& 0& {\mathrm{\Gamma }}_{\downarrow }\end{array}\right),\text{(6.33b)}& {\mathrm{\Lambda }}_{\mathrm{\Xi }}& =\left(\begin{array}{cccc}{\mathrm{\Gamma }}_{\text{sf}}& 0& 0& -{\mathrm{\Gamma }}_{\text{sf}}\\ 0& {\mathrm{\Gamma }}_{\tau }& 0& 0\\ 0& 0& {\mathrm{\Gamma }}_{\tau }^{\ast }& 0\\ -{\mathrm{\Gamma }}_{\text{sf}}& 0& 0& {\mathrm{\Gamma }}_{\text{sf}}\end{array}\right).\end{array}$$ The coefficients of these two matrices are linked to the characteristic length scales of the material. ${\mathrm{\Gamma }}_{\uparrow (\downarrow )}=1/l_{\uparrow (\downarrow )}$ is the mean free path for majority (minority) electrons. It can be rewritten in terms of the mean free path $l_{\ast }={(1/l_{\uparrow }+1/l_{\downarrow })}^{-1}$ and spin asymmetry of the material $\beta$, introduced by Valet-Fert in [121] as ${\mathrm{\Gamma }}_{\uparrow (\downarrow )}=(1\mp \beta )/(2l^{\ast })$. These length scales are connected to the spin-dependent resistivities of the Valet-Fert theory [121] as $l_{\uparrow (\downarrow )}=R_{\text{sh}}/{\rho }_{\uparrow (\downarrow )}$ and $l_{\ast }=R_{\text{sh}}/(4{\rho }_{\ast })$. ${\mathrm{\Gamma }}_{\text{mx,eff}}={\mathrm{\Gamma }}_{\text{mx}}/2+({\mathrm{\Gamma }}_{\uparrow }+{\mathrm{\Gamma }}_{\downarrow })/2$ is the component describing the behavior of the electrons with spin direction transverse to the local magnetization. ${\mathrm{\Gamma }}_{\text{mx}}=1/l_{\perp }-i/l_{\text{L}}$ is of ballistic origin, and describes the spin-dependent part of the transport, with the Larmor precession length and spin penetration length introduced in (6.24). The second part is related to the spin-independent part of the transport. ${\mathrm{\Gamma }}_{\text{sf}}=l_{\ast }/(4l_{\text{sf}}^2)$ is connected to spin-flip processes, and $l_{\text{sf}}$ is referred to as spin diffusion length. Finally, ${\mathrm{\Gamma }}_{\tau }={\mathrm{\Gamma }}_{\text{mx}}/2+2{\mathrm{\Gamma }}_{\text{sf}}$ describes the behavior of the non-conserved part of the transverse spins, including the fact that they are also affected by spin-flip processes.

It is then possible to recover charge and spin potentials (${\mu }_{\text{c}}$ and $\mathit{\mu }$, respectively) and charge and spin currents ($j_c$ and $\mathbf{j}$, respectively) from the variables ${\mu }_{\text{4}}$ and $j_{\text{4}}$. The process, described in detail in [125], involves a change of basis based on the set of $4\times 4$ matrices ${\text{I}}_{\text{ij}}={\sigma }_{\text{i}}\otimes {\sigma }_{\text{j}}^{\ast }$, where $i,j=0,1,2,3$, ${\sigma }_0$ is the identity matrix and ${\sigma }_{1,2,3}$ are the Pauli matrices. The resulting equations can be expressed in the same units employed for the Zhang-Levi-Fert formalism (6.6) through a change of variables. All the potentials and currents described until now are expressed in units of energy. One can express the charge potentials in units of V as $V=e{\mu }_{\text{c}}$ and convert the spin potential to spin accumulation, expressed in units of A/m, as $\mathbf{S}=-\left({\mu }_{\text{B}}\sigma /\left(e^2{D}_{\text{e}}\right)\right)\mathit{\mu }$ [83]. The charge current expressed in the usual units of A/m${}^2$ is $J_{\text{C}}=\left(e/\left(4R_{\text{sh}}\right)\right)j_{\text{c}}$, and the spin polarization current density, expressed in units of $A/s$, is ${\mathbf{J}}_{\text{S}}=-\left(4{\mu }_{\text{B}}/\left(e^2{R}_{\text{sh}}\right)\right)\mathbf{j}$ [83], [123]. Furthermore, the relations between the Zhang-Levy-Fert and Valet-Fert theories derived in [83] can be extended to include the length scales coming from the CMRT derivation:

$$\begin{array}{}\text{(6.34a)}& \sigma & =\frac{1}{{\rho }_{\ast }(1-{\beta }^2)}\text{(6.34b)}& {\lambda }_{\text{sf}}^2& =\frac{l_{\text{sf}}^2}{1-{\beta }^2}\text{(6.34c)}& {\lambda }_{\text{J}}^2& =\frac{l_{\ast }{l}_{\text{L}}}{1-{\beta }^2}\text{(6.34d)}& {\lambda }_{\phi }^2& =\frac{l_{\perp }}{l_{\text{L}}}{\lambda }_J^2\text{(6.34e)}& {\lambda }^2& =\frac{l_{\ast }^2}{1-{\beta }^2}\end{array}$$ ${\lambda }_{\phi }$ is the spin dephasing length introduced in (6.1), $\lambda =\sqrt{D_{\text{e}}\tau }$ is the momentum relaxation length, and $\tau$ is the time for momentum relaxation [98], [129]. The spin asymmetry parameter $\beta$ of the Valert-Fert model is equivalent to the conductivity polarization ${\beta }_{\sigma }$ in the Zhang-Levy-Fert approach. Finally, the spin and charge drift-diffusion equations for transport along the x-direction obtained through the formalism of [125] are, after the change of variables:

$$\begin{array}{}\text{(6.35a)}& J_{\text{C}}=\sigma E+\frac{e}{{\mu }_{\text{B}}}{\beta }_{\text{D}}{D}_{\text{e}}\frac{\partial \mathbf{S}}{\partial x}\mathbf{\cdot }\mathbf{m}\text{(6.35b)}& {\mathbf{J}}_{\text{S}}-\frac{{\lambda }^2}{{\lambda }_J^2}\mathbf{m}\times {\mathbf{J}}_{\text{S}}-\frac{{\lambda }^2}{{\lambda }_{\phi }^2}\mathbf{m}\times (\mathbf{m}\times {\mathbf{J}}_{\text{S}})=-\frac{{\mu }_{\text{B}}}{e}{\beta }_{\sigma }\mathbf{m}\sigma E-D_{\text{e}}\frac{\partial \mathbf{S}}{\partial x}\text{(6.35c)}& \frac{\partial J_{\text{C}}}{\partial x}=0\text{(6.35d)}& -\frac{\partial {\mathbf{J}}_{\text{S}}}{\partial x}-D_{\text{e}}\frac{\mathbf{S}}{{\lambda }_{sf}^2}+D_{\text{e}}\frac{\mathbf{m}\times \mathbf{S}}{{\lambda }_J^2}+D_{\text{e}}\frac{\mathbf{m}\times (\mathbf{m}\times \mathbf{S})}{{\lambda }_{\phi }^2}=0\end{array}$$ A generalization of these equations to three-dimensions leads to the following expressions [119], [123], [125]:

$$\begin{array}{}\text{(6.36a)}& {\mathbf{J}}_{\text{C}}=\sigma \mathbf{E}+\frac{e}{{\mu }_{\text{B}}}{\beta }_{\text{D}}{D}_{\text{e}}{(\mathbf{\nabla }\mathbf{S})}^{\text{T}}\mathbf{m}\text{(6.36b)}& {\tilde{\mathbf{J}}}_{\text{S}}-\frac{{\lambda }^2}{{\lambda }_J^2}{\left[\mathbf{m}\right]}_{\times }{\tilde{\mathbf{J}}}_{\text{S}}-\frac{{\lambda }^2}{{\lambda }_{\phi }^2}{\left[\mathbf{m}\right]}_{\times \times }{\tilde{\mathbf{J}}}_{\text{S}}=-\frac{{\mu }_{\text{B}}}{e}{\beta }_{\sigma }\mathbf{m}\otimes \left(\sigma \mathbf{E}\right)-D_{\text{e}}\mathbf{\nabla }\mathbf{S}\text{(6.36c)}& \mathbf{\nabla }\mathbf{\cdot }{\mathbf{J}}_{\text{C}}=0\text{(6.36d)}& -\mathbf{\nabla }\mathbf{\cdot }{\tilde{\mathbf{J}}}_{\text{S}}-D_{\text{e}}\frac{\mathbf{S}}{{\lambda }_{sf}^2}+D_{\text{e}}\frac{\mathbf{m}\times \mathbf{S}}{{\lambda }_J^2}+D_{\text{e}}\frac{\mathbf{m}\times (\mathbf{m}\times \mathbf{S})}{{\lambda }_{\phi }^2}=0\end{array}$$ ${\left[\mathbf{m}\right]}_{\times }$ and ${\left[\mathbf{m}\right]}_{\times \times }$ indicate the matrices associated with the cross-product and double cross-product with $\mathbf{m}$:

$$\begin{array}{}\text{(6.37)}& {\left[\mathbf{m}\right]}_{\times }=\left(\begin{array}{ccc}0& -m_{\text{1}}& m_{\text{2}}\\ m_{\text{3}}& 0& -m_{\text{1}}\\ -m_{\text{2}}& m_{\text{1}}& 0\end{array}\right),{\left[\mathbf{m}\right]}_{\times \times }=\left(\begin{array}{ccc}-m_{\text{2}}^2-m_{\text{3}}^2& m_{\text{1}}{m}_{\text{2}}& m_{\text{1}}{m}_{\text{3}}\\ m_{\text{1}}{m}_{\text{2}}& -m_{\text{1}}^2-m_{\text{3}}^2& m_{\text{2}}{m}_{\text{3}}\\ m_{\text{1}}{m}_{\text{3}}& m_{\text{2}}{m}_{\text{3}}& -m_{\text{1}}^2-m_{\text{2}}^2\end{array}\right)\end{array}$$

The result for the divergence of the spin current in (6.36c) stems from the absence of electric current sources inside the metallic layers, due to the fast redistribution of any charge imbalance [24], [25]. Equation (6.36d) differs from (6.12) by the presence of the spin dephasing term, dependent on ${\lambda }_{\phi }$, which leads to the torque expression (6.1) presented at the beginning of the chapter. Another important difference of this formulation is the presence of the additional terms on the left side of equation (6.36b), which mix up the orthogonal spin current components depending on the local magnetization orientation. These terms originate from the underlying ballistic origin of the transverse spin precession or dephasing [123], and depend mainly on the ratio between the momentum relaxation length $\lambda$ and the transverse absorption lengths ${\lambda }_{\text{J}}$ and ${\lambda }_{\phi }$. For transition metal ferromagnets, $\lambda$ and ${\lambda }_J$ have the same order of magnitude [98]. The present dissertation, following [25], focuses on a finite element implementation which does not include these additional terms, showing how an appropriate treatment of the tunneling layer and tuning of the system parameters are able to reproduce the most important properties of the torque expected in MTJs, while still retaining the capability of obtaining all torque contributions in several ferromagnetic layers from a unified expression. Some considerations about the possible influence of the additional spin current terms on the computed spin torque will be presented in Section 8.5 of Chapter 8.

By employing the spin current expression (6.8), the final set of equations defining the spin accumulation and torque is:

$$\begin{array}{}\text{(6.38a)}& -\mathbf{\nabla }\mathbf{\cdot }{\tilde{\mathbf{J}}}_{\text{S}}-D_{\text{e}}\frac{\mathbf{S}}{{\lambda }_{sf}^2}-{\mathbf{T}}_{\text{S}}=0\text{(6.38b)}& {\tilde{\mathbf{J}}}_{\text{S}}=-\frac{{\mu }_B}{e}{\beta }_{\sigma }\mathbf{m}\otimes ({\mathbf{J}}_{\text{C}}-\frac{e}{{\mu }_B}{\beta }_D{D}_e{(\mathbf{\nabla }\mathbf{S})}^{\text{T}}\mathbf{m})-D_e\mathbf{\nabla }\mathbf{S}\text{(6.38c)}& {\mathbf{T}}_{\text{S}}=-D_{\text{e}}\frac{\mathbf{m}\times \mathbf{S}}{{\lambda }_J^2}-D_{\text{e}}\frac{\mathbf{m}\times (\mathbf{m}\times \mathbf{S})}{{\lambda }_{\phi }^2}\end{array}$$

6.2 Finite Element Implementation

The LLG equation, coupled with the spin and charge drift-diffusion formalism for the computation of the torques, allows to describe the magnetization dynamics of structures containing an arbitrary number of different layers. Analytical solutions to these set of partial differential equations (PDEs) are only possible in simplified scenarios, and numerical methods are necessary to resolve the dynamics in a more general sense. In Chapter 4, a FD implementation of the LLG equation was shown. The FD solver is, however, only able to compute the magnetization dynamics of a thin free layer, with an approximate elliptical shape and thickness of one cell. The FE method is naturally able to handle meshes with complex geometries and several material domains [130], [131], and was therefore employed for the implementation of a solver capable of handling charge, spin accumulation and magnetization dynamics. The implementation was carried out by employing the open-source C++ FE library MFEM [132], [133]. The software responsible for the solution of the spin and charge drift-diffusion equations was developed by the author of the present dissertation. The following sections provide a brief overview of the basic principles of the FE method, followed by the weak formulation of both the LLG and spin and charge drift-diffusion equations, necessary in order to obtain a FE solution.

6.2.1 The Finite Element Method

The FE method is a numerical approach for the computation of approximate solutions to partial differential equations. In order to solve a problem, the simulation domain is first divided into an equivalent system of smaller and simpler bodies, referred to as finite elements, interconnected at nodes (points common to two or more elements) and boundary lines or surfaces. This is achieved by the construction of a mesh of the object, containing a finite number of points. The FE formulation of a given problem approximates the solution of differential equations with the solution of a system of algebraic equations. Instead of solving the problem for the entire body in one operation, the equations are first formulated for each element, and are then combined to obtain a solution valid for the whole domain.

Weak Formulation

The first step in the derivation of a FE representation of a given set of PDEs is to write the equations in the so-called weak formulation. To give an example of the derivation of the weak form of a given problem, the following Poisson equation is considered:

$$\begin{array}{}\text{(6.39)}& -{\mathrm{\nabla }}^{2}u=q\text{in }\mathrm{\Omega }\end{array}$$

$\mathrm{\Omega }$ indicates the integration domain, and $u$ is the unknown quantity being solved for. The boundary conditions are

$$\begin{array}{}\text{(6.40a)}& u& =u_{\text{D}}& \text{on }\partial {\mathrm{\Omega }}_{\text{D}},\text{(6.40b)}& \mathbf{\nabla }u\mathbf{\cdot }\mathbf{n}& =g& \text{on }\partial {\mathrm{\Omega }}_{\text{N}},\end{array}$$ where $\mathbf{n}$ is the external boundary normal, $\partial {\mathrm{\Omega }}_{\text{D}}$ is the portion of the external boundary $\partial \mathrm{\Omega }$ where Dirichlet boundary conditions (6.40a) are applied, and $\partial {\mathrm{\Omega }}_{\text{N}}$ is the portion of the external boundary where Neumann conditions (6.40b) are applied. A Dirichlet condition prescribes the value of the solution, while a Neumann condition prescribes the value of the normal derivative of the solution. A boundary condition needs to be specified for each external boundary of the domain. To obtain the weak formulation, both sides of (6.39) are multiplied by smooth functions, referred to as test functions, and integrated over the simulation domain:

$$\begin{array}{}\text{(6.41)}& -{\int }_{\mathrm{\Omega }}{\mathrm{\nabla }}^{2}uv\mathrm{d}\mathbf{x}={\int }_{\mathrm{\Omega }}qv\mathrm{d}\mathbf{x}\end{array}$$

$v$ indicates a test function. By applying Gauss theorem, stating that ${\int }_{\mathrm{\Omega }}\mathbf{\nabla }\mathbf{\cdot }\mathbf{p}\mathrm{d}\mathbf{x}$ $=$ ${\int }_{\partial \mathrm{\Omega }}\mathbf{p}\mathbf{\cdot }\mathbf{n}\mathrm{d}\mathbf{x}$, to the function $\mathbf{p}=\mathbf{\nabla }uv$, and by considering the product rule of derivation, one obtains the following equivalence:

$$\begin{array}{}\text{(6.42)}& -{\int }_{\mathrm{\Omega }}{\mathrm{\nabla }}^{2}uv={\int }_{\mathrm{\Omega }}\mathbf{\nabla }u\mathbf{\cdot }\mathbf{\nabla }v\mathrm{d}\mathbf{x}-{\int }_{\partial \mathrm{\Omega }}\mathbf{\nabla }u\mathbf{\cdot }\mathbf{n}v\mathrm{d}\mathbf{x}\end{array}$$

By putting this expression back into (6.41), and by considering the boundary conditions (6.40), one obtains the following equation:

$$\begin{array}{}\text{(6.43)}& {\int }_{\mathrm{\Omega }}\mathbf{\nabla }u\mathbf{\cdot }\mathbf{\nabla }v\mathrm{d}\mathbf{x}-{\int }_{\partial {\mathrm{\Omega }}_{\text{D}}}\mathbf{\nabla }u\mathbf{\cdot }\mathbf{n}v\mathrm{d}\mathbf{x}={\int }_{\mathrm{\Omega }}qv\mathrm{d}\mathbf{x}+{\int }_{\partial {\mathrm{\Omega }}_{\text{N}}}gv\mathrm{d}\mathbf{x}\end{array}$$

All the terms containing the unknown quantity $u$ where put on the left-hand side, while the known terms are on the right-hand side. The Dirichlet conditions are then assigned in a strong sense, imposing $u=u_{\text{D}}$ on $\partial {\mathrm{\Omega }}_{\text{D}}$, and the test function $v$ is chosen as to satisfy $v=0$ on $\partial {\mathrm{\Omega }}_{\text{D}}$. This way, the weak formulation of the Poisson equation takes the form

$$\begin{array}{}\text{(6.44)}& {\int }_{\mathrm{\Omega }}\mathbf{\nabla }u\mathbf{\cdot }\mathbf{\nabla }v\mathrm{d}\mathbf{x}={\int }_{\mathrm{\Omega }}qv\mathrm{d}\mathbf{x}+{\int }_{\partial {\mathrm{\Omega }}_{\text{N}}}gv\mathrm{d}\mathbf{x}.\end{array}$$

The test function and the solution are both assumed to belong to Hilbert spaces, and in the so-called Galerkin method, employed for the results presented in this work, they belong to the same Hilbert space. In particular, a proper choice for such space, referred to as $V$, is the following:

$$\begin{array}{}\text{(6.45)}& V:=\{v\in L_2(\mathrm{\Omega })|\mathbf{\nabla }v\in [L^2(\mathrm{\Omega }){]}^d\mathrm{\&}v{|}_{\begin{array}{c}\partial \mathrm{\Omega }\end{array}}\in L^2(\partial \mathrm{\Omega })\}\end{array}$$

$d$ is the number of spacial dimensions under consideration, and $L^2(\mathrm{\Omega })$ is the set of $L^2$-integrable functions in $\mathrm{\Omega }$. The weak form (6.44) is required to hold for all test functions $v$ in $V$. This formulation is referred to as "weak" because it only requires equality in an integral sense, while in the original form (6.39) all the terms must be well defined in all points.

From the Weak Formulation to a System of Algebraic Equations

The next step is the derivation of a numerical approximation of the weak form (6.44). For this purpose, the original domain $\mathrm{\Omega }$ is divided into smaller regular subdomains $T$. The set of all subdomains contained in $\mathrm{\Omega }$ is labeled $\mathcal{T}$ and referred to as triangulation. Typical choices for the shape of the subdomains $T$ are triangles ot squares in two dimensions and tetrahedrons or hexahedrons in three dimensions, but more complex shapes can also be employed. The set of points defining the discretization, referred to as nodes, is indicated as $\mathcal{N}=\left\{{\mathbf{x}}_{\text{i}}\right\}$, where ${\mathbf{x}}_{\text{i}}$ is the coordinate vector of node i. By means of the triangulation, a subspace of $V$ where to look for the approximate numerical solution $u_{\text{h}}\approx u$ can be defined as

$$\begin{array}{}\text{(6.46)}& V_{\text{h}}:=\{v\in C(\mathrm{\Omega })|v{|}_{\begin{array}{c}T\end{array}}\text{ is affine }\mathrm{\forall }T\in \mathcal{T}\}.\end{array}$$

$V_{\text{h}}$ is referred to as finite element space. Each function $v_{\text{h}}$ belonging to $V_{\text{h}}$ is uniquely defined by its values at the nodes. The set of nodes can be divided into two subsets, one containing the nodes at the Dirichlet boundary, ${\mathcal{N}}_{\text{D}}$, and the other containing the remaining (free) nodes, ${\mathcal{N}}_{\text{f}}$. In order to get a discretized version of (6.44), the functions $v_{\text{h}}$ can be represented in terms of the nodal basis $\left\{{\phi }_{\text{i}}\right\}$ of $V_{\text{h}}$, characterized as

$$\begin{array}{}\text{(6.47)}& {\phi }_{\text{i}}\left({\mathbf{x}}_{\text{j}}\right)={\delta }_{ij}.\end{array}$$

The test functions and the finite element solution $u_{\text{h}}$ can be represented with respect to this basis:

$$\begin{array}{}\text{(6.48)}& u_{\text{h}}\left(\mathbf{x}\right)=\sum _{i=1}^N{u}_{\text{i}}{\phi }_{\text{i}}\left(\mathbf{x}\right)\end{array}$$

$N$ is the total number of nodes, and $u_{\text{i}}=u_{\text{h}}({\mathbf{x}}_{\text{i}})$ are the values assumed by the approximate solution at the nodes. Fig. 6.1 illustrates a possible approximation of the solution $u$ through linear basis functions in a one-dimensional scenario. $u$ could be taken to represent the electrical potential in a rod with non-uniform charge distribution. The rod is split into eight equal parts (which represent its triangulation), and the basis functions associated with each node have a value of one at their respective node, and 0 at all other nodes. It can be seen how $u_{\text{h}}$, as defined by (6.48), can be used to approximate the continuous solution $u$. One of the advantages of employing the finite element method is that it offers great freedom in the selection of discretization, both in the choice of elements and basis functions [130]. Moreover, it is possible to locally refine the mesh to obtain a better approximations in the regions where the solution varies quickly.

Determining the coefficients $u_{\text{i}}$ gives a unique solution to the problem at hand. The values on ${\mathcal{N}}_{\text{D}}$ are determined by the Dirichlet boundary conditions:

$$\begin{array}{}\text{(6.49)}& u_{\text{i}}=u_{\text{D}}({\mathbf{x}}_{\text{i}})\mathrm{\forall }{\mathbf{x}}_{\text{i}}\in \partial {\mathrm{\Omega }}_{\text{D}}\end{array}$$

All the remaining values can be determined by inserting $u_{\text{h}}$ and $v_{\text{h}}$ instead of $u$ and $v$ in (6.44), and then expressing both of them in terms of their nodal values through (6.48). Fulfilling (6.44) for the whole space $V_{\text{h}}$ is equivalent to only considering the basis associated with the free nodes $\{{\phi }_{\text{i}}|{\mathbf{x}}_{\text{i}}\in {\mathcal{N}}_{\text{f}}\}$, leading to the following expression:

$$\begin{array}{}\text{(6.50)}& \sum _i\{{\int }_{\mathrm{\Omega }}\mathbf{\nabla }{\phi }_{\text{i}}\mathbf{\cdot }\mathbf{\nabla }{\phi }_{\text{j}}\mathrm{d}\mathbf{x}\}{u}_{\text{i}}={\int }_{\mathrm{\Omega }}q{\phi }_{\text{j}}\mathrm{d}\mathbf{x}+{\int }_{\partial {\mathrm{\Omega }}_{\text{N}}}g{\phi }_{\text{j}}\mathrm{d}\mathbf{x}\end{array}$$

This equivalence needs to be valid for the whole set of basis functions associated with the free nodes. A matrix $A\in {\mathbb{R}}^N\times {\mathbb{R}}^N$ for the right-hand side of (6.50) and a vector $f\in {\mathbb{R}}^N$ for the left-hand side can be defined as

$$\begin{array}{}\text{(6.51a)}& A_{\text{ji}}& ={\int }_{\mathrm{\Omega }}\mathbf{\nabla }{\phi }_{\text{i}}\mathbf{\cdot }\mathbf{\nabla }{\phi }_{\text{j}}\mathrm{d}\mathbf{x},\text{(6.51b)}& f_{\text{j}}& ={\int }_{\mathrm{\Omega }}q{\phi }_{\text{j}}\mathrm{d}\mathbf{x}+{\int }_{\partial {\mathrm{\Omega }}_{\text{N}}}g{\phi }_{\text{j}}\mathrm{d}\mathbf{x}.\end{array}$$ Both $A$ and $f$ can be represented in terms of their components acting on the Dirichlet and free nodes:

$$\begin{array}{}\text{(6.52)}& A=\left(\begin{array}{cc}{A}_{\text{DD}}& A_{\text{Df}}\\ A_{\text{fD}}& A_{\text{ff}}\end{array}\right),f=\left(\begin{array}{c}{f}_{\text{D}}\\ f_{\text{f}}\end{array}\right)\end{array}$$

As only neighboring nodes have overlapping basis functions, $A$ is a sparse matrix, with non-zero terms only around the diagonal. By considering the constraint (6.49), and the solution vector $u=(u_{\text{i}})\in {\mathbb{R}}^N$, $u={(u_{\text{D}},u_{\text{f}})}^{\text{T}}$, (6.50) can be written as the following system of linear equations:

$$\begin{array}{}\text{(6.53)}& \left(\begin{array}{cc}I& 0\\ A_{\text{fD}}& A_{\text{ff}}\end{array}\right)\left(\begin{array}{c}{u}_{\text{D}}\\ u_{\text{f}}\end{array}\right)=\left(\begin{array}{c}{u}_{\text{D}}\\ f_{\text{f}}\end{array}\right)\end{array}$$

$I$ stands for the identity matrix. As the $u_{\text{D}}$ values are known, the final symmetric system of $N_{\text{f}}$ equations for the determination of the unknown $u_{\text{f}}$ values is

$$\begin{array}{}\text{(6.54)}& A_{\text{ff}}{u}_{\text{f}}=f_{\text{f}}-A_{\text{fD}}{u}_{\text{D}}.\end{array}$$

6.2.2 Weak Formulation of the Micromagnetic Equations

The implementation of a FE algorithm for the solution of the LLG equation has been of great interest for the solution of micromagnetic problems with meshes of varying complexity. First schemes for a numerical approximation of the weak LLG solution were proposed in [134], [135], by considering only the exchange field contribution to the effective field. In reference [136] the so-called tangent plane integrator scheme, solving for the discrete time derivative, was introduced, and was later generalized to include the full effective field [137], [138]. The unconditional convergence of the tangent plane integrator scheme and of the FE implementation of the spin and charge drift-diffusion equations was proven in [139], and later applied to metallic spin-valves in [23]. Here, the weak formulation of the tangent plane integrator and of the spin and charge transport equations, employed to obtain the results presented in these thesis, is reported. The flow-chart of the simulation process is presented in Fig. 6.2.

LLG Equation

The tangent plane scheme solves the LLG equation for the magnetization derivative $\partial \mathbf{m}/\partial t=\mathbf{v}$. The weak formulation of the LLG equation comes from the expression

$$\begin{array}{}\text{(6.55)}& \alpha \frac{\partial \mathbf{m}}{\partial t}+\mathbf{m}\times \frac{\partial \mathbf{m}}{\partial t}={\gamma }_0{\mathbf{H}}_{\text{eff}}-{\gamma }_0(\mathbf{m}\mathbf{\cdot }{\mathbf{H}}_{\text{eff}})\mathbf{m},\end{array}$$

which can be obtained by cross-multiplying (3.9) with $\mathbf{m}$, and using the product rule $\mathbf{a}\times (\mathbf{b}\times \mathbf{c})=(\mathbf{c}\mathbf{\cdot }\mathbf{a})\mathbf{b}-(\mathbf{a}\mathbf{\cdot }\mathbf{b})\mathbf{c}$ together with the constraint $\left|\mathbf{m}\right|=1$. For the FE implementation, the magnetization is taken to be a piecewise affine, globally continuous function [83]. Each component belongs to the Sobolev space $H^1$. It is the space of functions in $L^2$ that additionally admit a weak gradient which also belongs to $L^2$ [138]. The weak gradient is a generalization of the gradient, defined with the help of the test functions and valid in an integral sense, and its existence is needed to be able to evaluate terms containing spatial derivatives in the weak formulation. The notation for the vector space of the magnetization is ${\mathbf{H}}^1$. Instead of looking for the solution $\mathbf{v}$ in the same space of the magnetization, the solution space $V_{\text{T}}$ is restricted to vectors tangent to the magnetization, so that $V_{\text{T}}=\{\mathit{w}\in {\mathbf{H}}^1|\mathbf{m}\mathbf{\cdot }\mathit{w}=0\}$. The test functions are also restricted to the same space, so that the weak formulation of (6.55) results in

$$\begin{array}{}\text{(6.56)}& {\int }_{\omega }(\alpha \mathbf{v}+\mathbf{m}\times \mathbf{v})\mathbf{\cdot }\mathit{w}\mathrm{d}\mathbf{x}={\gamma }_0{\int }_{\omega }{\mathbf{H}}_{\text{eff}}(\mathbf{m})\mathbf{\cdot }\mathit{w}\mathrm{d}\mathbf{x},\end{array}$$ where $\omega$ is the subdomain containing only the magnetic parts of the structure under analysis. The last term on the right-hand side of (6.55) is not present, as the test functions are restricted to the tangent space $V_{\text{T}}$. The time derivative $\mathbf{v}$ at a certain time $t^{\text{k}}$ is obtained by setting [83]

$$\begin{array}{}\text{(6.57)}& {\mathbf{m}}^{\text{k+1}}={\mathbf{m}}^{\text{k}}+\theta \delta t\mathbf{v},\end{array}$$

where $\delta t$ is the time-step, and $\theta$ is a parameter between 0 and 1, with the value 0 leading to a fully explicit scheme and 1 to a fully implicit one. Each effective field contribution can be treated with a different value of $\theta$. In the implementation employed for this work, $\theta$ is different from 0 only for the exchange field contribution (3.17), where the value 1 is employed for stability reasons [140]. The weak formulation employed for the computation of the magnetization dynamics by the FE solver, with the inclusion of the torque terms coming from (6.1), is then

$$\begin{array}{c}{\int }_{\omega }(\alpha \mathbf{v}+{\mathbf{m}}^{\text{k}}\times \mathbf{v})\mathbf{\cdot }\mathit{w}\mathrm{d}\mathbf{x}+\frac{2A\gamma }{M_S}\delta t{\int }_{\omega }\mathbf{\nabla }\mathbf{v}:\mathbf{\nabla }\mathit{w}\mathrm{d}x=\\ -\frac{2A\gamma }{M_S}{\int }_{\omega }\mathbf{\nabla }{\mathbf{m}}^{\text{k}}:\mathbf{\nabla }\mathit{w}\mathrm{d}\mathbf{x}+{\gamma }_0{\int }_{\omega }({\mathbf{H}}_{\text{ext}}+{\mathbf{H}}_{\text{ani}}({\mathbf{m}}^{\text{k}})+{\mathbf{H}}_{\text{demag}}({\mathbf{m}}^{\text{k}}))\mathbf{\cdot }\mathit{w}\mathrm{d}x+\\ \text{(6.58a)}& +\frac{D_{\text{e}}}{M_{\text{S}}}{\int }_{\omega }(\frac{{\mathbf{S}}^{\text{k}}}{{\lambda }_J^2}+\frac{{\mathbf{m}}^{\text{k}}\times {\mathbf{S}}^{\text{k}}}{{\lambda }_{\phi }^2})\mathbf{\cdot }\mathit{w}\mathrm{d}\mathbf{x},\text{(6.58b)}& {\mathbf{m}}^{\text{k+1}}=\frac{{\mathbf{m}}^{\text{k}}+\delta t\mathbf{v}}{|{\mathbf{m}}^{\text{k}}+\delta t\mathbf{v}|}.\end{array}$$ The right-hand side of (6.58b) is evaluated nodewise, $\omega$ indicates the magnetic subdomains, and $\mathbf{\nabla }\mathbf{a}:\mathbf{\nabla }\mathbf{b}=\sum _{\text{ij}}\left(\partial a_{\text{i}}/\partial x_{\text{j}}\right)\left(\partial b_{\text{i}}/\partial x_{\text{j}}\right)$ is the Frobenius inner product of two matrices. The exchange field contribution (3.17), together with (6.57), gives rise to the second term on the left-hand side and to the first term on the right-hand side. The boundary integrals arising from the weak formulation of this contribution are put to zero by the Neumann boundary condition (3.41), applied on the external boundary of the magnetic region $\partial \omega$. The equations are subject to the initial condition $\mathbf{m}(0)={\mathbf{m}}^0$. The system of equations resulting from the FE implementation of this weak formulation includes the tangent plane constraint $\mathbf{m}\mathbf{\cdot }\mathit{w}=0$, and the solution at each time-step is computed through a solver based on the generalized minimal residual (GMRES) method, provided by the library MFEM. The GMRES method is an iterative algorithm for the numerical solution of an indefinite nonsymmetric system of linear equations, as is the case of FE implementation of (6.58a) due to the presence of the cross-product terms. Material parameters that can differ from subdomain to subdomain are treated as piecewise constant functions, unless differently stated. The contribution of the demagnetizing field ${\mathbf{H}}_{\text{demag}}$ is evaluated only on the disconnected magnetic domain by using a hybrid approach combining the boundary element method and the FE method [141], [142].

Spin and Charge Drift-Diffusion Equations

The weak formulation for the computation of the spin accumulation is derived from (6.38). The equations are solved for the magnetization ${\mathbf{m}}^{\text{k}}$ at each time-step, and the resulting spin accumulation ${\mathbf{S}}^{\text{k}}$ is employed in (6.58a).

For the computation of the charge potential and current, a Laplace equation is solved for the whole domain $\mathrm{\Omega }$, with the weak formulation

$$\begin{array}{}\text{(6.59a)}& {\int }_{\mathrm{\Omega }}\sigma \mathbf{\nabla }V\mathbf{\cdot }\mathbf{\nabla }v\mathrm{d}\mathbf{x}=0,\text{(6.59b)}& {\int }_{\mathrm{\Omega }}{\mathbf{J}}_{\text{C}}\mathbf{\cdot }\mathit{v}\mathrm{d}\mathbf{x}=-{\int }_{\mathrm{\Omega }}\sigma \mathbf{\nabla }V\mathbf{\cdot }\mathit{v}\mathrm{d}\mathbf{x},\end{array}$$ where $v$ represents a test function belonging to $H^1$, while $\mathit{v}$ is a test function in ${\mathbf{H}}^1$. $V$ belongs to $H^1$, and equation (6.59b) is employed to obtain a projection of ${\mathbf{J}}_{\text{C}}$ in the ${\mathbf{H}}^1$ function space [83]. Dirichlet conditions are applied to prescribe the voltage at the contacts. The Neumann condition $\sigma \mathbf{\nabla }V\mathbf{\cdot }\mathbf{n}=0$ is assumed on external boundaries not containing an electrode. Special treatment of the conductivity $\sigma$ results in an implementation capable of reproducing the charge current dependence on the TMR and relative direction of the magnetization vectors in the FL and RL of an MTJ, as will be shown in Chapter 7. The solution of the system of equations resulting from the FE implementation of (6.59) is computed through a solver based on the conjugate gradient (CG) method, provided by the library MFEM. The CG method is an algorithm for the numerical solution of systems of linear equations whose matrix is positive-definite, as is the case of the FE implementation of (6.59).

The weak formulation of the spin drift-diffusion equations presented in [23] is generalized to apply to (6.38), which includes the additional spin dephasing term, resulting in the following expression:

$$\begin{array}{c}{D}_{\text{e}}{\int }_{\mathrm{\Omega }}(\mathbf{\nabla }\mathbf{S}-{\beta }_{\sigma }{\beta }_{\text{D}}\mathbf{m}\otimes \left({(\mathbf{\nabla }\mathbf{S})}^{\text{T}}\mathbf{m}\right)):\mathbf{\nabla }\mathit{v}\mathrm{d}\mathbf{x}+\\ +D_{\text{e}}{\int }_{\mathrm{\Omega }}(\frac{\mathbf{S}}{{\lambda }_{sf}^2}+\frac{\mathbf{S}\times \mathbf{m}}{{\lambda }_J^2}+\frac{\mathbf{m}\times (\mathbf{S}\times \mathbf{m})}{{\lambda }_{\phi }^2})\mathbf{\cdot }\mathit{v}\mathrm{d}\mathbf{x}=\\ \text{(6.60)}& \frac{{\mu }_B}{e}{\beta }_{\sigma }{\int }_{\omega }(\mathbf{m}\otimes {\mathbf{J}}_{\text{C}}):\mathbf{\nabla }\mathit{v}\mathrm{d}\mathbf{x}-\frac{{\mu }_B}{e}{\beta }_{\sigma }{\int }_{\partial \mathrm{\Omega }\cap \partial \omega }\left((\mathbf{m}\otimes {\mathbf{J}}_{\text{C}})\mathbf{n}\right)\mathbf{\cdot }\mathit{v}\mathrm{d}\mathbf{x}\end{array}$$ $\mathit{v}$ represents again a test function belonging to ${\mathbf{H}}^1$, $\mathbf{S}$ also belongs to ${\mathbf{H}}^1$, and $\partial \mathrm{\Omega }\cap \partial \omega$ indicates the shared external boundary of the whole domain $\mathrm{\Omega }$ and magnetic subdomain $\omega$. The boundary integrals arising from the weak formulation are put to zero by the Neumann condition $(\mathbf{\nabla }\mathbf{S})\mathbf{n}=\mathbf{0}$, assumed on external boundaries. For contacting regions longer than the spin-flip length, this condition is equivalent to an exponential decay of $\mathbf{S}$ towards the electrodes [23], [25]. The charge current ${\mathbf{J}}_{\text{C}}$ is the one computed from (6.59). The solution of the system of equations resulting from the FE implementation of (6.60) is computed through a solver based on the GMRES method, provided by the library MFEM, as due to the cross product terms the FE implementation of (6.60) results in a nonsymmetric system of linear equations.