Finite Element Approach for the Simulation of Modern MRAM Devices

Abstract

Because of their nonvolatile nature and simple structure, the interest in MRAM devices has been steadily growing in recent years. Reliable simulation tools, capable of handling complex geometries composed of multiple materials, provide valuable help in improving the design of MRAM cells. In this work, we describe a solver based on the finite element implementation of the Landau–Lifshitz–Gilbert equation coupled to the spin and charge drift-diffusion formalism. The torque acting in all layers from different contributions is computed from a unified expression. In consequence of the versatility of the finite element implementation, the solver is applied to switching simulations of recently proposed structures based on spin-transfer torque, with a double reference layer or an elongated and composite free layer, and of a structure combining spin-transfer and spin-orbit torques.

1. Introduction

As the scaling of conventional CMOS technology shows signs of saturation, due to the increase in standby power consumption and leakage currents, the employment of nonvolatile memory components which do not require the memory bits to be refreshed becomes increasingly appealing [1]. One of the most promising candidates as a nonvolatile replacement is magnetoresistive random access memory (MRAM). It possesses a simple structure and is directly compatible with CMOS back-end of line processes. It has shown to be promising for several applications, for example in stand-alone memories and in the embedded automotive and Internet of Things fields, and is expected to replace charge-based devices in frame buffer memory and slow SRAM [2,3,4,5]. Moreover, MRAM devices have shown to be interesting for cryogenic applications, especially for employment in quantum computing systems [6,7,8].

The core of an MRAM cell is the magnetic tunnel junction (MTJ), a stack of two ferromagnetic (FM) layers separated by an oxide layer. The properties of the two FM layers are such that the magnetization in one of them, the reference layer (RL), is fixed, while in the other one, the free layer (FL), it can be switched between the two stable parallel (P) and anti-parallel (AP) states. The resistance of the stack can be employed to store the binary information, as it is higher in the AP state. The percentage difference between the resistance of the two stable states is labeled the tunneling magnetoresistance ratio (TMR).

The writing process in modern devices is performed by relying on spin-transfer torque (STT), spin-orbit torque (SOT), or a combination of both of them. In STT switching, the torque is generated by a current flowing through the MTJ. The electrons are polarized by the RL and transfer their polarization to the FL magnetization, providing the torque [9,10,11]. Examples of structures based on STT are reported in Figure 1a–c. In SOT switching, the torque is generated by passing the current through a heavy metal line (HM) below the FL. The spin Hall effect (SHE) produces a spin current orthogonal to the charge one, which is absorbed by the FL magnetization, providing the torque [12,13]. An example of a three-terminal structure combining STT and SOT is shown in Figure 1d.

The design of modern MRAM cells can be supported by the development of reliable simulation tools. The magnetization dynamics are described by the Landau–Lifshitz–Gilbert (LLG) equation, which must be supplied with a term describing the spin torque. The drift-diffusion formalism offers a way of computing different torque contributions in all the ferromagnetic layers from a unified expression [14,15,16]. The finite element (FE) method, being able to handle structures with several domains of different materials and complex geometries, represents an optimal choice for computing a numerical solution to the micromagnetic equations in modern MRAM cells. In this work, we present an FE-based implementation of the LLG equation coupled with the drift-diffusion formalism, extended to include the charge and torque properties expected in MTJs. The solver was developed by employing the open-source C++ FE library MFEM [17,18], and is applied to the simulation of recently proposed structures based both on STT and SOT switching. The source code is available as an open-source repository [19].

Figure 1. Four examples of multi-layer MRAM cell design: (a) standard STT-MRAM with single MTJ; (b) double RL STT-MRAM, where the second RL provides additional torque to reduce the critical voltage required for switching [20]; (c) ultra-scaled STT-MRAM, where the FM layers are elongated and additional oxide layers are added to improve scalability and benefit from the shape anisotropy [21]; (d) SOT-assisted STT-MRAM, where the switching process is kick-started by an initial current pulse in the HM [22].

Figure 1. Four examples of multi-layer MRAM cell design: (a) standard STT-MRAM with single MTJ; (b) double RL STT-MRAM, where the second RL provides additional torque to reduce the critical voltage required for switching [20]; (c) ultra-scaled STT-MRAM, where the FM layers are elongated and additional oxide layers are added to improve scalability and benefit from the shape anisotropy [21]; (d) SOT-assisted STT-MRAM, where the switching process is kick-started by an initial current pulse in the HM [22].

2. Micromagnetic Modeling

The LLG equation for the description of the magnetization dynamics was first derived by Landau and Lifshitz in 1935 [23] and reworked by Gilbert in a more treatable form in 1955 [24]. With the inclusion of the spin torque ${\mathbf{T}}_{\mathrm{S}}$, it takes the form

$$\frac{\partial \mathbf{m}}{\partial t}=-\left|\gamma \right|{\mu }_0\mathbf{m}\times {\mathbf{H}}_{\mathrm{eff}}+\alpha \mathbf{m}\times \frac{\partial \mathbf{m}}{\partial t}+\frac{1}{M_{\mathrm{S}}}{\mathbf{T}}_{\mathrm{S}}$$

(1)

where $\mathbf{m}=\mathbf{M}/M_{\mathrm{S}}$ is the unit vector in the direction of the local magnetization, $M_{\mathrm{S}}$ is the saturation magnetization, $\gamma$ is the gyromagnetic ratio, and ${\mu }_0$ is the vacuum permeability. ${\mathbf{H}}_{\mathrm{eff}}$ is an effective magnetic field including the contributions of an externally applied field, the exchange coupling, the anisotropy field, and the demagnetizing field. The effects of temperature, which can be included by an additional effective field contribution describing thermal fluctuations [25], are not considered for the switching results presented in this work. The main effect of their inclusion would be to reduce the incubation time necessary for the switching process, while the behavior would remain qualitatively similar. While the external field ${\mathbf{H}}_{\mathrm{ext}}$ can be simply added as an input parameter, the other contributions are intrinsic to the ferromagnets and must be computed from material parameters.

The exchange coupling can be modeled through a field which tends to keep the magnetization vectors aligned throughout the magnetic domain, described by the expression

$${\mathbf{H}}_{\mathrm{exc}}=\frac{2A_{\mathrm{exc}}}{{\mu }_0{M}_{\mathrm{S}}}{\nabla }^2\mathbf{m},$$

(2)

where $A_{\mathrm{exc}}$ is the exchange coefficient.

Modern MRAM cells utilize MTJs with magnetization perpendicular to the stack, by virtue of both interface and shape anisotropy contributions [21]. While the latter is taken into account by the demagnetizing field, the former can be included as a uniaxial anisotropy field with the expression

$${\mathbf{H}}_{\mathrm{ani}}=\frac{2K_{\mathrm{ani}}}{{\mu }_0{M}_{\mathrm{S}}}\left(\mathbf{a}·\mathbf{m}\right)\mathbf{a},$$

(3)

where $\mathbf{a}$ is a unit vector in the direction perpendicular to the stack and $K_{\mathrm{ani}}$ is the anisotropy coefficient, which can be computed from the interface anisotropy $K_{\mathrm{int}}$ as $K_{\mathrm{ani}}=K_{\mathrm{int}}/d_{\mathrm{FM}}$, where $d_{\mathrm{FM}}$ is the thickness of the ferromagnetic layer under consideration.

The demagnetizing field can be computed from the scalar magnetic potential $u_{\mathrm{m}}$ as

$${\mathbf{H}}_{\mathrm{demag}}=-\mathbf{\nabla }{u}_{\mathrm{m}}$$

(4)

$u_{\mathrm{m}}$ is obtained through the solution of the Poisson equation

$$-{\nabla }^2{u}_{\mathrm{m}}=-M_{\mathrm{S}}\mathbf{\nabla }·\mathbf{m},$$

(5)

with the constraint of $u_{\mathrm{m}}$ decaying to zero as $\mathcal{O}(1/|\mathbf{x}{|}^2)$ outside the magnetic domain and the boundary condition $[\mathbf{\nabla }{u}_{\mathrm{m}}·\mathbf{n}]=-M_{\mathrm{S}}\mathbf{m}·\mathbf{n}$, where $\mathbf{n}$ is the unit vector normal to the boundary and $[\dots ]$ denotes a discontinuity across the boundary.

In the presence of a single thin FL, the torque acting on the magnetization can be described by simplified expressions [26,27,28,29], derived by Slonczewski [9,11]. A more general form of the torque term, which allows an arbitrary number of ferromagnetic and non-magnetic layers to be dealt with, can be obtained by computing the non-equilibrium spin accumulation in the structure under study through the solution of spin and charge transport equations.

Spin and Charge Transport

The torque term ${\mathbf{T}}_{\mathrm{S}}$ entering (1) can be computed from the spin accumulation through the following expression [16,30,31,32]:

$${\mathbf{T}}_{\mathrm{S}}=-D_{\mathrm{e}}\frac{\mathbf{m}\times \mathbf{S}}{{\lambda }_J^2}-D_{\mathrm{e}}\frac{\mathbf{m}\times \left(\mathbf{m}\times \mathbf{S}\right)}{{\lambda }_{\phi }^2}$$

(6)

The first term describes the 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 }$. $D_{\mathrm{e}}$ is the electron diffusion coefficient. 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}}_{\mathrm{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 [33]. A solution for $\mathbf{S}$ in all non-magnetic and ferromagnetic layers of an MRAM cell can be obtained by means of the spin and charge drift-diffusion formalism.

Spin and charge drift-diffusion equations in multilayer structures with arbitrary magnetization orientation were reported by S. Zhang, P. Levy, and A. Fert [34], with both the precession and decay of the transverse spin accumulation components governed by the exchange length ${\lambda }_{\mathrm{J}}$. Another possible mechanism for the absorption of the transverse components is the dephasing process [32,35]. The behavior of the spin accumulation with both precessional and dephasing terms was described in terms of the Continuous Random Matrix Theory (CRMT) in [35], and the equivalence of the CMRT and the spin and charge drift-diffusion formalism was shown. The resulting equations for spin and charge currents are [14,16]:

$$\begin{array}{cc}\hfill {\mathbf{J}}_{\mathrm{C}}& =\sigma \mathbf{E}+\frac{e}{{\mu }_{\mathrm{B}}}{\beta }_{\mathrm{D}}{D}_{\mathrm{e}}{\left(\mathbf{\nabla }\mathbf{S}\right)}^{\mathrm{T}}\mathbf{m},\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill {\tilde{\mathbf{J}}}_{\mathrm{S}}& =-\frac{{\mu }_{\mathrm{B}}}{e}{\beta }_{\sigma }\mathbf{m}\otimes \left(\sigma \mathbf{E}\right)-D_{\mathrm{e}}\mathbf{\nabla }\mathbf{S},\hfill \end{array}$$

(8)

where ${\tilde{\mathbf{J}}}_{\mathrm{S}}$ is the spin current tensor, ${\mathbf{J}}_{\mathrm{C}}$ is the charge current density, ${\mu }_{\mathrm{B}}$ is the Bohr magneton, e is the elementary charge, $\sigma$ is the conductivity, and ⊗ is the outer product. ${\beta }_{\sigma }=\left({\sigma }^{\uparrow }-{\sigma }^{\downarrow }\right)/\left({\sigma }^{\uparrow }+{\sigma }^{\downarrow }\right)$ and ${\beta }_{\mathrm{D}}=\left(D_{\mathrm{e}}^{\uparrow }-D_{\mathrm{e}}^{\downarrow }\right)/\left(D_{\mathrm{e}}^{\uparrow }+D_{\mathrm{e}}^{\downarrow }\right)$ are the conductivity and diffusion polarization parameters, respectively, with ${\sigma }^{\uparrow }$, $D_{\mathrm{e}}^{\uparrow }$ (${\sigma }^{\downarrow }$, $D_{\mathrm{e}}^{\downarrow }$) the conductivity and diffusion coefficient for the majority (minority) electrons. $\mathbf{\nabla }\mathbf{S}$ is the vector gradient of $\mathbf{S}$, with components ${\left(\mathbf{\nabla }\mathbf{S}\right)}_{\mathrm{ij}}=\partial S_{\mathrm{i}}/\partial x_{\mathrm{j}}$, and the term ${\left(\mathbf{\nabla }\mathbf{S}\right)}^{\mathrm{T}}\mathbf{m}$ is a vector with components ${\left({\left(\mathbf{\nabla }\mathbf{S}\right)}^{\mathrm{T}}\mathbf{m}\right)}_{\mathrm{i}}={\sum }_{\mathrm{j}}\left(\partial S_{\mathrm{i}}/\partial x_{\mathrm{j}}\right)m_j$. The spin current can be expressed in terms of the charge current by inserting (7) into (8):

$${\tilde{\mathbf{J}}}_{\mathrm{S}}=-\frac{{\mu }_B}{e}{\beta }_{\sigma }\mathbf{m}\otimes \left({\mathbf{J}}_{\mathrm{C}}-\frac{e}{{\mu }_B}{\beta }_D{D}_e{\left(\mathbf{\nabla }\mathbf{S}\right)}^{\mathrm{T}}\mathbf{m}\right)-D_e\mathbf{\nabla }\mathbf{S}$$

(9)

The equation of motion for the spin accumulation is given by

$$\frac{\partial \mathbf{S}}{\partial t}=-\mathbf{\nabla }·{\tilde{\mathbf{J}}}_{\mathrm{S}}-D_{\mathrm{e}}\frac{\mathbf{S}}{{\lambda }_{sf}^2}-{\mathbf{T}}_{\mathrm{S}},$$

(10)

where ${\lambda }_{\mathrm{sf}}$ is the spin-flip length and $\nabla ·{\tilde{\mathbf{J}}}_{\mathrm{S}}$ is the divergence of ${\tilde{\mathbf{J}}}_{\mathrm{S}}$, with components ${\left(\nabla ·{\tilde{\mathbf{J}}}_{\mathrm{S}}\right)}_{\mathrm{i}}={\sum }_{\mathrm{j}}\partial J_{\mathrm{S},\mathrm{ij}}/\partial x_{\mathrm{j}}$. As the typical time scale for the magnetization motion is three orders of magnitude larger than the spin accumulation one [34], it is sufficient to consider a steady-state expression for the spin accumulation. This assumption was numerically verified in [36]. With $\partial \mathbf{S}/\partial t=0$, the equation describing the spin accumulation becomes

$$-\mathbf{\nabla }·{\tilde{\mathbf{J}}}_{\mathrm{S}}-D_{\mathrm{e}}\frac{\mathbf{S}}{{\lambda }_{sf}^2}-{\mathbf{T}}_{\mathrm{S}}=0$$

(11)

3. Finite Element Implementation

The presented set of equations allows the magnetization dynamics of structures containing an arbitrary number of layers of different materials to be described. The FE method, a numerical approach for the computation of approximate solutions to partial differential equations, is naturally able to handle meshes with complex geometries and several domains of different materials [37,38], and was therefore employed for the implementation of a solver capable of handling charge, spin accumulation, and the magnetization dynamics. The implementation was carried out by employing the open-source C++ FE library MFEM [17,18].

The first schemes for a FE implementation of the LLG equation, which considered only the contribution of the exchange coupling, were proposed in [39,40]. A new FE algorithm, referred to as the tangent plane integrator scheme, was introduced in [41] and later generalized in [42,43] to include the contributions of the demagnetizing and anisotropy fields. The unconditional convergence of an algorithm coupling the LLG equation with a FE implementation of the spin and charge drift-diffusion formalism was proven in [44], and the scheme was later successfully applied to metallic spin valves in [14]. We report here an extension of the scheme to MTJs, which includes the spin dephasing contribution and allows both the TMR effect and the expected torque properties to be reproduced [45,46].

3.1. Charge Current Solution

For the computation of the charge current entering (9), the Laplace equation is solved:

$$\begin{array}{c}-\mathbf{\nabla }·\left(\sigma \mathbf{\nabla }V\right)=0\end{array}$$

(12)

$$\begin{array}{c}{\mathbf{J}}_{\mathrm{C}}=-\sigma \mathbf{\nabla }V\end{array}$$

(13)

where V is the electric potential. Dirichlet conditions are applied to prescribe the voltage at the contacts, and the Neumann condition $\sigma \mathbf{\nabla }V·\mathbf{n}=0$ is assumed on external boundaries not containing an electrode, with $\mathbf{n}$ the unit vector normal to the boundary. In order to be able to reproduce the TMR effect, the tunnel barrier is treated as a poor conductor whose local conductivity depends on the relative magnetization orientation in the RL and FL [45,46]:

$${\sigma }_{\mathrm{TB}}={\sigma }_0\left(1+P_{\mathrm{FL}}{P}_{\mathrm{RL}}{\mathbf{m}}_{\mathrm{RL}}·{\mathbf{m}}_{\mathrm{FL}}\right)$$

(14)

where $P_{\mathrm{RL}}$ and $P_{\mathrm{FL}}$ are the Slonczewski polarization parameters [11,47], ${\sigma }_0=({\sigma }_P+{\sigma }_{AP})/2$ is the angle independent portion of the conductivity, ${\sigma }_{P\left(AP\right)}$ is the conductivity in the parallel (anti-parallel) state, and ${\mathbf{m}}_{\mathrm{RL}\left(\mathrm{FL}\right)}$ is the magnetization of the RL(FL) close to the interface. $P_{\mathrm{RL}}$ and $P_{\mathrm{FL}}$ are related to the TMR by Julliere’s formula [48]:

$$\mathrm{TMR}=\frac{R_{\mathrm{AP}}-R_{\mathrm{P}}}{R_{\mathrm{P}}}=\frac{2P_{\mathrm{FL}}{P}_{\mathrm{RL}}}{1-P_{\mathrm{FL}}{P}_{\mathrm{RL}}}$$

(15)

where $R_{\mathrm{P}\left(\mathrm{AP}\right)}$ is the resistance in the parallel (anti-parallel) state.

In order to derive an FE representation, the equations must be written in the so-called weak formulation. For the presented Laplace equation, by using Gauss’s theorem and applying the Neumann boundary conditions, the weak form reduces to

$${\int }_{\mathsf{\Omega }}\sigma \mathbf{\nabla }V·\mathbf{\nabla }v\mathrm{d}\mathbf{x}=0$$

(16)

The test function v and the solution are both assumed to belong to the Sobolev space $H^1$, so that both they and their weak gradients are $L^2$-integrable [43].

In the FE approximation, in order to obtain a discretized version of (16), the original domain $\mathsf{\Omega }$ is divided into smaller regular elements. The discrete solution $V_{\mathrm{h}}$ is defined by its values on the elements’ nodes:

$$V_{\mathrm{h}}\left(\mathbf{x}\right)=\sum _{i=1}^N{V}_{\mathrm{i}}{\phi }_{\mathrm{i}}\left(\mathbf{x}\right)$$

(17)

where N is the total number of nodes, $V_{\mathrm{i}}=V_{\mathrm{h}}\left({\mathbf{x}}_{\mathrm{i}}\right)$ are the values assumed by the approximate solution at the nodes, and ${\mathbf{x}}_{\mathrm{i}}$ is the coordinate vector of node i. ${\phi }_{\mathrm{i}}$ is an affine function of the nodal basis of the mesh, characterized as

$${\phi }_{\mathrm{i}}\left({\mathbf{x}}_{\mathrm{j}}\right)={\delta }_{ij}$$

(18)

Figure 2 illustrates an example of the approximation of a function u through linear basis functions in a one-dimensional scenario.

With the given nodal basis decomposition, the original problem can thus be rewritten as the following system of linear equations:

$$\underline{A_{\mathrm{V}}}\underline{V_{\mathrm{h}}}=\underline{0}$$

(19)

where $\underline{A_{\mathrm{V}}}\in {\mathbb{R}}^N\times {\mathbb{R}}^N$ is the matrix representation of the left-hand side (LHS) of (16), and $\underline{V_{\mathrm{h}}}$ is a vector in ${\mathbb{R}}^N$ composed of the nodal values of $V_{\mathrm{h}}$. As only neighboring nodes have overlapping basis functions, $\underline{A_{\mathrm{V}}}$ is a sparse matrix, with non-zero terms only around the diagonal.

Figure 2. Representation of a continuous function u and its finite element approximation $u_{\mathrm{h}}$ in a one-dimensional setting. The basis functions for all the nodes are reported at the bottom of the graph. The basis function and solution value associated with the node $x_2$ are labeled ${\phi }_2$ and $u_2$, respectively.

Figure 2. Representation of a continuous function u and its finite element approximation $u_{\mathrm{h}}$ in a one-dimensional setting. The basis functions for all the nodes are reported at the bottom of the graph. The basis function and solution value associated with the node $x_2$ are labeled ${\phi }_2$ and $u_2$, respectively.

The weak formulation of (13) is

$${\int }_{\mathsf{\Omega }}{\mathbf{J}}_{\mathrm{C}}·v\mathrm{d}\mathbf{x}=-{\int }_{\mathsf{\Omega }}\sigma \mathbf{\nabla }V·v\mathrm{d}\mathbf{x}$$

(20)

By choosing the test function $\mathit{v}$ so that its components belong to $H^1$, noting the space containing $v$ as ${\mathbf{H}}^1$, this equation allows a projection to be obtained of ${\mathbf{J}}_{\mathrm{C}}$ in the ${\mathbf{H}}^1$ function space [33] so that it can be readily employed for the computation of the spin accumulation. The resulting system of linear equations is

$$\underline{A_{\mathrm{J}}}\underline{{\mathbf{J}}_{\mathrm{C},\mathrm{h}}}=\underline{f_{\mathrm{J}}},$$

(21)

where $\underline{A_{\mathrm{J}}}\in {\mathbb{R}}^{3N}\times {\mathbb{R}}^{3N}$ is the matrix coming from the LHS of (20) and $\underline{f_{\mathrm{J}}}\in {\mathbb{R}}^{3N}$ is the vector coming from the right-hand side (RHS).

The solution to both systems of equations is computed through a solver based on the conjugate gradient method [49], designed for the numerical solution of systems of linear equations whose matrices are positive definite, provided by the library MFEM.

In the scope of the MFEM library, only the data associated with a local element can be accessed during the assembly of the system matrices, while the computation of (14) in the TB requires knowledge of the magnetization vectors in the neighboring FM layers. In order to obtain access to the magnetization values, the coefficient describing the TB conductivity is initialized as follows:

• For each point inside the TB where the conductivity needs to be computed, referred to as an integration point, the solver loops through the integration points of the RL and FL elements closer to the interfaces.
• The RL and FL points near to or at the interface with coordinates closest to the TB point are selected.
• The integration point number and element number associated with the nearest RL and FL points are mapped to the coordinates of the TB points.

In a transient simulation, the search is carried out only during the initialization of the solver. At every time step, the data necessary for the computation of (14) can be accessed through the generated maps, without the need to repeat the search procedure. The computed charge current, consistent with the TMR effect, can then be employed to obtain a solution to the spin accumulation equation.

3.2. Spin Accumulation Solution

The weak form of Equation (11), with the spin current expressed as (9) and the spin torque as (6), takes the form

$$\begin{array}{c}-D_{\mathrm{e}}{\int }_{\mathsf{\Omega }}\left(\mathbf{\nabla }·\left(\mathbf{\nabla }\mathbf{S}-{\beta }_{\sigma }{\beta }_{\mathrm{D}}\mathbf{m}\otimes \left({\left(\mathbf{\nabla }\mathbf{S}\right)}^{\mathrm{T}}\mathbf{m}\right)\right)\right)·v\mathrm{d}\mathbf{x}+\\ +D_{\mathrm{e}}{\int }_{\mathsf{\Omega }}\left(\frac{\mathbf{S}}{{\lambda }_{sf}^2}+\frac{\mathbf{S}\times \mathbf{m}}{{\lambda }_J^2}+\frac{\mathbf{m}\times \left(\mathbf{S}\times \mathbf{m}\right)}{{\lambda }_{\phi }^2}\right)·v\mathrm{d}\mathbf{x}=\\ \frac{{\mu }_B}{e}{\beta }_{\sigma }{\int }_{\omega }\left(\mathbf{\nabla }·\left(\mathbf{m}\otimes {\mathbf{J}}_{\mathrm{C}}\right)\right)·v\mathrm{d}\mathbf{x},\end{array}$$

(22)

where $v$ represents again a test function belonging to ${\mathbf{H}}^1$, $\mathbf{S}$ also belongs to ${\mathbf{H}}^1$, and $\omega$ indicates a subdomain composed of only the ferromagnetic layers. By applying Gauss’s theorem, the first term on the LHS of (22) becomes

$$\begin{array}{c}{D}_{\mathrm{e}}{\int }_{\mathsf{\Omega }}\left(\mathbf{\nabla }\mathbf{S}-{\beta }_{\sigma }{\beta }_{\mathrm{D}}\mathbf{m}\otimes \left({\left(\mathbf{\nabla }\mathbf{S}\right)}^{\mathrm{T}}\mathbf{m}\right)\right):\mathbf{\nabla }v\mathrm{d}\mathbf{x}+\\ -{\int }_{\partial \mathsf{\Omega }}\left(\left(\mathbf{\nabla }\mathbf{S}\right)\mathbf{n}-{\beta }_{\sigma }{\beta }_{\mathrm{D}}\left(\mathbf{m}\otimes \mathbf{m}\right)\left(\left(\mathbf{\nabla }\mathbf{S}\right)\mathbf{n}\right)\right)·v\mathrm{d}\mathbf{x},\end{array}$$

(23)

where $\mathbf{\nabla }\mathbf{a}:\mathbf{\nabla }\mathbf{b}={\sum }_{\mathrm{ij}}\left(\partial a_{\mathrm{i}}/\partial x_{\mathrm{j}}\right)\left(\partial b_{\mathrm{i}}/\partial x_{\mathrm{j}}\right)$ is the Frobenius inner product of two matrices. By assuming the natural Neumann condition $\left(\mathbf{\nabla }\mathbf{S}\right)\mathbf{n}=\mathbf{0}$ on all external boundaries of the whole domain $\mathsf{\Omega }$, the integrals on $\partial \mathsf{\Omega }$ are put to zero. If the contacts are longer than the spin-flip length, the condition is equivalent to an exponential decay of $\mathbf{S}$ towards the electrodes [14,16].

Gauss’s theorem can also be applied to the RHS term of (22), obtaining

$$-\frac{{\mu }_B}{e}{\beta }_{\sigma }{\int }_{\omega }\left(\mathbf{m}\otimes {\mathbf{J}}_{\mathrm{C}}\right):\mathbf{\nabla }v\mathrm{d}\mathbf{x}+\frac{{\mu }_B}{e}{\beta }_{\sigma }{\int }_{\partial \mathsf{\Omega }\cap \partial \omega }\left(\left(\mathbf{m}\otimes {\mathbf{J}}_{\mathrm{C}}\right)\mathbf{n}\right)·v\mathrm{d}\mathbf{x}$$

(24)

$\partial \omega$ indicates the external boundaries of the magnetic subdomains, and $\partial \mathsf{\Omega }\cap \partial \omega$ indicates the shared external boundaries of the subdomain $\omega$ and the whole domain $\mathsf{\Omega }$.

3.2.1. Tunneling Spin Current

The inclusion of appropriate boundary conditions at the TB interface with the RL and FL, together with the employment of a low diffusion coefficient inside the TB, allows the expected properties of the torque acting in MTJs to be reproduced [46]. The additional boundary conditions to be added to the RHS of (22) read

$${\mathrm{RHS}}_{\mathrm{TB}}=-{\int }_{RL|TB}{\mathbf{J}}_{\mathrm{S},\mathrm{TB}}·v\mathrm{d}\mathbf{x}+{\int }_{TB|FL}{\mathbf{J}}_{\mathrm{S},\mathrm{TB}}·v\mathrm{d}\mathbf{x},$$

(25)

where $RL\left|TB\right(TB\left|FL\right)$ indicates the interface of the TB with the RL(FL). These internal boundary conditions prescribe the difference in spin current between the FM layers and the TB, according to the spin current polarization generated by the tunneling process. The tunneling spin current ${\mathbf{J}}_{\mathrm{S},\mathrm{TB}}$ can be expressed as [47,50]

$$\begin{array}{c}{\mathbf{J}}_{\mathrm{S},\mathrm{TB}}=-\frac{{\mu }_B}{e}\frac{{\mathbf{J}}_{\mathrm{C},\mathrm{TB}}·\mathbf{n}}{1+P_{\mathrm{RL}}{P}_{\mathrm{FL}}{\mathbf{m}}_{\mathrm{RL}}·{\mathbf{m}}_{\mathrm{FL}}}(a_{\mathrm{mx}}{P}_{\mathrm{RL}}{\mathbf{m}}_{\mathrm{RL}}+a_{\mathrm{mx}}{P}_{\mathrm{FL}}{\mathbf{m}}_{\mathrm{FL}}+\\ \left.+1/2\left(P_{\mathrm{RL}}{P}_{\mathrm{RL}}^{\eta }-P_{\mathrm{FL}}{P}_{\mathrm{FL}}^{\eta }\right){\mathbf{m}}_{\mathrm{RL}}\times {\mathbf{m}}_{\mathrm{FL}}\right),\end{array}$$

(26)

where $P_{RL}^{\eta }$ and $P_{FL}^{\eta }$ are out-of-plane polarization parameters [47], $a_{mx}$ describes the influence of the interface spin-mixing conductance on the transmitted in-plane spin current [50], ${\mathbf{J}}_{\mathrm{C},\mathrm{TB}}$ is the electric current density at the interface, $\mathbf{n}$ is the interface normal, and ${\mathbf{m}}_{\mathrm{RL}\left(\mathrm{FL}\right)}$ is the local value of the RL(FL) magnetization at the interface.

3.2.2. Spin Hall Effect

When a charge current flows through an HM layer with strong spin–orbit coupling, it generates a spin current perpendicular to it, carrying a spin polarization perpendicular to the direction of both spin and charge currents [12]. This process is known as the spin Hall effect (SHE). If the FL is deposited right above the HM, this spin current can be employed to provide the torque necessary for switching.

In order to reproduce the SHE, the following term must be added to the spin current expression (9) [12,16]:

$${\tilde{\mathbf{J}}}_{\mathrm{S},\mathrm{SHE}}=-{\theta }_{\mathrm{SHA}}\frac{{\mu }_{\mathrm{B}}}{e}\epsilon {\mathbf{J}}_{\mathrm{C}}$$

(27)

where ${\theta }_{\mathrm{SHA}}$ is the spin Hall angle, and $\epsilon$ is the rank-3 unit antisymmetric tensor [16]. With the boundary condition

$$(\mathbf{\nabla }\mathbf{S})\mathbf{n}=-{\theta }_{\mathrm{SHA}}\frac{D_{\mathrm{e}}}{\sigma }\frac{{\mu }_{\mathrm{B}}}{e}\left(\epsilon {\mathbf{J}}_{\mathrm{C}}\right)\mathbf{n}$$

(28)

the weak formulation with the updated spin current expression is the same as (22) with the addition of the following RHS term:

$${\mathrm{RHS}}_{\mathrm{SHE}}=-{\int }_{\mathrm{HM}}{\theta }_{\mathrm{SHA}}\frac{{\mu }_{\mathrm{B}}}{e}\left(\epsilon {\mathbf{J}}_{\mathrm{C}}\right):\mathbf{\nabla }v\mathrm{d}\mathbf{x}$$

(29)

The integral is performed only over the HM layer.

3.2.3. Complete Weak Formulation

The complete weak formulation of the spin accumulation equation takes the form

$$\begin{array}{c}{D}_{\mathrm{e}}{\int }_{\mathsf{\Omega }}\left(\mathbf{\nabla }\mathbf{S}-{\beta }_{\sigma }{\beta }_{\mathrm{D}}\mathbf{m}\otimes \left({\left(\mathbf{\nabla }\mathbf{S}\right)}^{\mathrm{T}}\mathbf{m}\right)\right):\mathbf{\nabla }v\mathrm{d}\mathbf{x}+\\ +D_{\mathrm{e}}{\int }_{\mathsf{\Omega }}\left(\frac{\mathbf{S}}{{\lambda }_{sf}^2}+\frac{\mathbf{S}\times \mathbf{m}}{{\lambda }_J^2}+\frac{\mathbf{m}\times \left(\mathbf{S}\times \mathbf{m}\right)}{{\lambda }_{\phi }^2}\right)·v\mathrm{d}\mathbf{x}=\\ -\frac{{\mu }_B}{e}{\beta }_{\sigma }{\int }_{\omega }\left(\mathbf{m}\otimes {\mathbf{J}}_{\mathrm{C}}\right):\mathbf{\nabla }v\mathrm{d}\mathbf{x}+\frac{{\mu }_B}{e}{\beta }_{\sigma }{\int }_{\partial \mathsf{\Omega }\cap \partial \omega }\left(\left(\mathbf{m}\otimes {\mathbf{J}}_{\mathrm{C}}\right)\mathbf{n}\right)·v\mathrm{d}\mathbf{x}+\\ -{\int }_{RL|TB}{\mathbf{J}}_{\mathrm{S},\mathrm{TB}}·v\mathrm{d}\mathbf{x}+{\int }_{TB|FL}{\mathbf{J}}_{\mathrm{S},\mathrm{TB}}·v-{\int }_{\mathrm{HM}}{\theta }_{\mathrm{SHA}}\frac{{\mu }_{\mathrm{B}}}{e}\left(\epsilon {\mathbf{J}}_{\mathrm{C}}\right):\mathbf{\nabla }v\mathrm{d}\mathbf{x}\end{array}$$

(30)

where $\mathbf{S}$ is the spin accumulation, $\mathbf{m}$ is the unit magnetization vector, ${\mathbf{J}}_{\mathrm{C}}$ is the charge current density, ${\mathbf{J}}_{\mathrm{S},\mathrm{TB}}$ is the tunneling spin current (26), $D_{\mathrm{e}}$ is the electron diffusion coefficient, ${\beta }_{\sigma }$ and ${\beta }_{\mathrm{D}}$ are polarization parameters, ${\lambda }_{\mathrm{sf}}$ is the spin-flip length, ${\lambda }_{\mathrm{J}}$ is the exchange length, ${\lambda }_{\phi }$ is the dephasing length, and ${\theta }_{\mathrm{SHA}}$ is the spin Hall angle. The system of linear equations to be solved in the FE implementation of (30) is

$$\underline{A_{\mathrm{S}}}\underline{{\mathbf{S}}_{\mathrm{h}}}=\underline{f_{\mathrm{S}}},$$

(31)

where $\underline{A_{\mathrm{S}}}\in {\mathbb{R}}^{3N}\times {\mathbb{R}}^{3N}$ is the matrix coming from the LHS of (30) and $\underline{f_{\mathrm{S}}}\in {\mathbb{R}}^{3N}$ is the vector coming from the RHS. The solution of this system of equations is computed through a solver based on the generalized minimal residual (GMRES) method [51], provided by the library MFEM. The GMRES method is designed for indefinite non-symmetric systems of linear equations, as is the case for (31) due to the presence of the cross-product terms in (30). Material parameters that can change between the different subdomains are treated as piecewise constant coefficients.

As is the case for (14), the inclusion of the additional boundary conditions (25) in the MFEM implementation demands special care. The computation of the boundary terms requires knowledge of the magnetization vector on the opposite interfaces. In order to obtain access to these values, the coefficient describing the boundary integral is initialized as follows:

• For each integration point on the RL|TB interface requiring the computation of the tunneling spin current, the solver loops through the integration points of the TB|FL interface.
• The TB|FL point with coordinates closest to the RL|TB one is selected.
• The integration point number and the element number associated with the found TB|FL point are mapped to the coordinates of the RL|TB one.
• The mapping procedure is repeated for the TB|FL interface.

In a transient simulation, the search is carried out only during the initialization of the solver. At every time step, the data necessary for the computation of (26) can be accessed through the generated maps, without the need to repeat the search procedure.

3.3. Magnetization Dynamics Solution

In the tangent plane scheme, the quantity being solved for is the magnetization derivative $\partial \mathbf{m}/\partial t=\mathsf{v}$, with the constraint $\mathbf{m}·\mathsf{v}=0$. By cross-multiplying (1) with $\mathbf{m}$, using the product rule $\mathbf{a}\times \left(\mathbf{b}\times \mathbf{c}\right)=\left(\mathbf{c}·\mathbf{a}\right)\mathbf{b}-\left(\mathbf{a}·\mathbf{b}\right)\mathbf{c}$ and the constraint $\left|\mathbf{m}\right|=1$, the LLG equation can be rewritten in a form employed to derive a weak formulation for the tangent plane scheme:

$$\begin{array}{c}\alpha \frac{\partial \mathbf{m}}{\partial t}+\mathbf{m}\times \frac{\partial \mathbf{m}}{\partial t}=\left|\gamma \right|{\mu }_0{\mathbf{H}}_{\mathrm{eff}}+\frac{D_{\mathrm{e}}}{M_{\mathrm{S}}{\lambda }_{\mathrm{J}}^2}\mathbf{S}+\frac{D_{\mathrm{e}}}{M_{\mathrm{S}}{\lambda }_{\phi }^2}\mathbf{m}\times \mathbf{S}+\\ -\left|\gamma \right|{\mu }_0\left(\mathbf{m}·{\mathbf{H}}_{\mathrm{eff}}\right)\mathbf{m}-\frac{D_{\mathrm{e}}}{M_{\mathrm{S}}{\lambda }_{\mathrm{J}}^2}\left(\mathbf{m}·\mathbf{S}\right)\mathbf{m}\end{array}$$

(32)

The magnetization is taken to belong to ${\mathbf{H}}^1$, while the solution $\mathsf{v}$ and the test functions $\mathsf{w}$ are restricted to a space of vectors orthogonal to the magnetization, $U_{\mathrm{T}}=\mathsf{w}\in {\mathbf{H}}^1|\mathbf{m}·\mathsf{w}=0$. The weak formulation of (32) is then

$$\begin{array}{c}{\int }_{\omega }\left(\alpha \mathsf{v}+\mathbf{m}\times \mathsf{v}\right)·\mathsf{w}\mathrm{d}\mathbf{x}=\left|\gamma \right|{\mu }_0{\int }_{\omega }\left({\mathbf{H}}_{\mathrm{ext}}+{\mathbf{H}}_{\mathrm{exc}}+{\mathbf{H}}_{\mathrm{ani}}+{\mathbf{H}}_{\mathrm{demag}}\right)·\mathsf{w}\mathrm{d}\mathbf{x}+\\ +\frac{D_{\mathrm{e}}}{M_{\mathrm{S}}}{\int }_{\omega }\left(\frac{\mathbf{S}}{{\lambda }_J^2}+\frac{\mathbf{m}\times \mathbf{S}}{{\lambda }_{\phi }^2}\right)·\mathsf{w}\mathrm{d}\mathbf{x},\end{array}$$

(33)

where the last two terms on the RHS of (32) are not present, as their scalar product with the test functions, belonging to the tangent space $U_{\mathrm{T}}$, is zero. By using Gauss’s theorem, the weak form of expression (2) for the exchange contribution can be written as

$$\frac{2A_{\mathrm{exc}}}{{\mu }_0{M}_{\mathrm{S}}}{\int }_{\omega }{\nabla }^2\mathbf{m}·\mathsf{w}\mathrm{d}\mathbf{x}=-\frac{2A_{\mathrm{exc}}}{{\mu }_0{M}_{\mathrm{S}}}{\int }_{\omega }\mathbf{\nabla }\mathbf{m}:\mathbf{\nabla }\mathsf{w}\mathrm{d}\mathbf{x}+\frac{2A_{\mathrm{exc}}}{{\mu }_0{M}_{\mathrm{S}}}{\int }_{\partial \omega }\left(\left(\mathbf{\nabla }\mathbf{m}\right)\mathbf{n}\right)·\mathsf{w}\mathrm{d}\mathbf{x}$$

(34)

The natural Neumann condition $(\mathbf{\nabla }\mathbf{m})\mathbf{n}=\mathbf{0}$ is assumed on $\partial \omega$, so that the boundary integral on the RHS is put to zero.

With the given weak formulation, the time derivative $\mathsf{v}$ at a certain time $t^{\mathrm{k}}$ is obtained by setting [33]

$${\mathbf{m}}^{\mathrm{k}+1}={\mathbf{m}}^{\mathrm{k}}+\theta \delta t\mathsf{v},$$

(35)

where $\delta t$ indicates the time step and $\theta$ is a parameter ranging from 0 to 1. A value of 0 leads to a fully explicit scheme, while a value of 1 gives a fully implicit one. The value of $\theta$ can differ between each effective field contribution. In the implementation reported here, only the exchange field contribution is treated implicitly with $\theta =1$, as this leads to a better stability of the scheme [52]. The weak formulation employed by the FE solver to compute the magnetization dynamics is then expressed as

$$\begin{array}{c}{\int }_{\omega }\left(\alpha \mathsf{v}+{\mathbf{m}}^{\mathrm{k}}\times \mathsf{v}\right)·\mathsf{w}\mathrm{d}\mathbf{x}+\frac{2A_{\mathrm{exc}}\left|\gamma \right|}{M_S}\delta t{\int }_{\omega }\mathbf{\nabla }\mathsf{v}:\mathbf{\nabla }\mathsf{w}\mathrm{d}\mathbf{x}=\\ -\frac{2A_{\mathrm{exc}}\left|\gamma \right|}{M_S}{\int }_{\omega }\mathbf{\nabla }{\mathbf{m}}^{\mathrm{k}}:\mathbf{\nabla }\mathsf{w}\mathrm{d}\mathbf{x}+{\gamma }_0{\int }_{\omega }\left({\mathbf{H}}_{\mathrm{ext}}+{\mathbf{H}}_{\mathrm{ani}}+{\mathbf{H}}_{\mathrm{demag}}\right)·\mathsf{w}\mathrm{d}\mathbf{x}+\\ +\frac{D_{\mathrm{e}}}{M_{\mathrm{S}}}{\int }_{\omega }\left(\frac{{\mathbf{S}}^{\mathrm{k}}}{{\lambda }_J^2}+\frac{{\mathbf{m}}^{\mathrm{k}}\times {\mathbf{S}}^{\mathrm{k}}}{{\lambda }_{\phi }^2}\right)·\mathsf{w}\mathrm{d}\mathbf{x},\end{array}$$

(36)

$$\begin{array}{c}{\mathbf{m}}^{\mathrm{k}+1}=\frac{{\mathbf{m}}^{\mathrm{k}}+\delta t\mathsf{v}}{|{\mathbf{m}}^{\mathrm{k}}+\delta t\mathsf{v}|},\end{array}$$

(37)

with the initial condition $\mathbf{m}\left(0\right)={\mathbf{m}}^0$. Equation (37) is evaluated nodewise. The additional tangent plane constraint $\mathbf{m}·\mathsf{w}=0$ leads to the following saddle point problem [53]:

$$\left(\begin{array}{cc}\underline{A_{\mathrm{M}}}& \underline{C_{\mathrm{M}}^{\mathrm{T}}}\\ \underline{C_{\mathrm{M}}}& \underline{0}\end{array}\right)\left(\begin{array}{c}\underline{{\mathsf{v}}_{\mathrm{h}}}\\ \underline{\lambda }\end{array}\right)=\left(\begin{array}{c}\underline{f_{\mathrm{M}}}\\ \underline{0}\end{array}\right)$$

(38)

where $\underline{A_{\mathrm{M}}}\in {\mathbb{R}}^{3N}\times {\mathbb{R}}^{3N}$ is the matrix coming from the LHS of (36), $\underline{f_{\mathrm{M}}}\in {\mathbb{R}}^{3N}$ is the vector coming from its RHS, $\underline{\lambda }$ is a scalar field, and $\underline{C_{\mathrm{M}}}\in {\mathbb{R}}^N\times {\mathbb{R}}^{3N}$ implements the constraint. A solution of (38) is computed at each time step through a solver based on the GMRES method.

3.4. Demagnetizing Field

The demagnetizing field contribution needs to be computed from the magnetic potential as (4), which in turn is obtained from (5). A direct FE implementation of the latter requires a large computational domain surrounding the magnetic material in order to ensure the proper decay properties of the computed potential. There have been various solutions proposed to solve this open-boundary problem [54,55,56,57], with the truncation of the external domain surrounding the magnetic one at a certain distance being the most straightforward. This approach, however, decreases computational efficiency, as it requires the inclusion of additional degrees of freedom.

High accuracy and reduced computational costs can be achieved by employing a hybrid approach, combining the FE method with the boundary element method (FEM-BEM), allowing $u_{\mathrm{m}}$ to be computed only in the magnetic subdomains [58]. The potential is first split into two parts:

$$u_{\mathrm{m}}=u_{\mathrm{m},1}+u_{\mathrm{m},2}$$

(39)

$u_{\mathrm{m},1}$ satisfies (5) inside the magnetic subdomain, with the boundary condition $\mathbf{\nabla }{u}_{\mathrm{m},1}·\mathbf{n}=M_{\mathrm{S}}\mathbf{m}·\mathbf{n}$, and is zero outside of it. $u_{\mathrm{m},2}$ satisfies the Laplace equation

$${\nabla }^2{u}_{\mathrm{m},2}=0,$$

(40)

with the boundary conditions $[\mathbf{\nabla }{u}_{\mathrm{m},2}·\mathbf{n}]=0$ and $\left[u_{m,2}\right]=u_{m,1}$, where $[\dots ]$ denotes a discontinuity across the boundary. By having $u_{\mathrm{m},2}\to 0$ for $\left|\mathbf{x}\right|\to \infty$, potential theory leads to the following relation between $u_{\mathrm{m},1}$ and $u_{\mathrm{m},2}$ [59]:

$$u_{\mathrm{m},2}={\int }_{\partial \omega }{u}_{\mathrm{m},1}\frac{\partial }{\partial \mathbf{n}}\frac{1}{|{\mathbf{x}}^{\prime }-\mathbf{x}|}\mathrm{d}{\mathbf{x}}^{\prime }$$

(41)

The decomposition allows $u_{\mathrm{m},1}$ to be computed by solving (5) only in the disconnected magnetic layers. The boundary value of $u_{\mathrm{m},2}$ is obtained by solving (41), and is then used as a Dirichlet condition for (40), which is also computed only inside the magnetized portions of the structure.

The weak formulation of (5), after applying Gauss’s theorem and the boundary condition, results in

$${\int }_{\omega }\mathbf{\nabla }{u}_{\mathrm{m},1}·\mathbf{\nabla }v\mathrm{d}\mathbf{x}=M_{\mathrm{S}}{\int }_{\omega }\mathbf{m}·\mathbf{\nabla }v\mathrm{d}\mathbf{x},$$

(42)

while that of (40) is

$${\int }_{\omega }\mathbf{\nabla }{u}_{\mathrm{m},2}·\mathbf{\nabla }v\mathrm{d}\mathbf{x}=0$$

(43)

Both the magnetic potential and the test functions belong to $H^1$. The FE implementation results in a system of equations analogous to the one employed for the charge potential.

A BEM approach is employed to discretize (41) on the magnetic boundary, resulting in the following matrix-vector multiplication:

$$\underline{u_{\mathrm{m},2}^{\mathrm{bdr}}}=\underline{B_{\mathrm{M}}}\underline{u_{\mathrm{m},1}^{\mathrm{bdr}}}$$

(44)

The matrix $\underline{B_{\mathrm{M}}}$ belongs to ${\mathbb{R}}^{N_{\mathrm{bdr}}}\times {\mathbb{R}}^{N_{\mathrm{bdr}}}$ and $\underline{u_{\mathrm{m},1}^{\mathrm{bdr}}}$, $\underline{u_{\mathrm{m},2}^{\mathrm{bdr}}}$ belong to ${\mathbb{R}}^{N_{\mathrm{bdr}}}$, with $N_{\mathrm{bdr}}$ being the number of boundary nodes of the magnetic subdomains. Even though $\underline{B_{\mathrm{M}}}$ is a dense matrix, the employment of matrix compression algorithms [60] can significantly reduce the memory demands [58]. The matrix compression algorithms and BEM functionalities are implemented by employing the H2Lib library [61]. The demagnetizing field is finally computed as the gradient of the magnetic potential $u_{\mathrm{m}}$ by using the same projection approach employed for the charge current in (20).

An example of the magnetic potential and field computed in a structure with three disconnected ferromagnetic layers is shown in Figure 3. Without interaction between the layers, the potential would only vary linearly along the magnetization direction. When applying the described FEM-BEM approach, the interactions are taken into account and the magnetic potential in each layer is shifted due to the stray field contributions of the neighboring ferromagnetic segments. The presented implementation allows the demagnetizing field acting in structures containing multiple ferromagnetic layers, as is typical of modern MRAM cells, to be readily computed.

4. Device Simulation

Recently proposed devices are composed of several layers of ferromagnetic materials, non-magnetic spacers, and tunnel barriers, in order to reduce switching currents and cell size. Due to the capability of computing the torque acting in all layers from a unified expression, the presented FE solver is suitable for the simulation of such structures. The following sections report the results of switching simulations performed in the structures of Figure 1. The parameters employed are presented in

Table 1. Material parameters.

LLG parameters	Value
Saturation magnetization ($M_s$)	$0.81\times {10}^6$ A/m
Exchange constant ($A_{exc}$)	$2.0\times {10}^{-11}$ J/m
Interface anisotropy ($K_{int}$)	$1.29\times {10}^{-3}$ J/m${}^2$
Gilbert damping constant ($\alpha$)	$0.02$
Drift-diffusion parameters	Value
Conductivity polarization, ${\beta }_{\sigma }$	$0.52$
Diffusion polarization, ${\beta }_D$	$0.7$
FM diffusion coefficient, $D_{e,FM}$	${10}^{-3}$ m${}^2$/s
NM diffusion coefficient, $D_{e,NM}$	${10}^{-2}$ m${}^2$/s
HM diffusion coefficient, $D_{e,HM}$	$1.1\times {10}^{-3}$ m${}^2$/s
TB diffusion coefficient, $D_S$	$2.0\times {10}^{-8}$ m${}^2$/s
FM conductivity ${\sigma }_{FM}$	$4.0\times {10}^6$ S/m
NM conductivity ${\sigma }_{NM}$	$5.0\times {10}^6$ S/m
HM conductivity ${\sigma }_{HM}$	$7.0\times {10}^6$ S/m
FM spin-flip length, ${\lambda }_{sf,FM}$	10 nm
NM spin-flip length, ${\lambda }_{sf,NM}$	10 nm
HM spin-flip length, ${\lambda }_{sf,HM}$	$1.4$ nm
Spin exchange length, ${\lambda }_J$	$0.8$ nm
Spin dephasing length, ${\lambda }_{\phi }$	$0.4$ nm
Spin Hall angle, ${\theta }_{\mathrm{SHA}}$	0.19
TB resistance standard and double RL STT-MTJ	Value
Resistance parallel ($R_P$)	$4.3\times {10}^3$ k$\mathsf{\Omega }$
Resistance anti-parallel ($R_{AP}$)	$9.1\times {10}^3$ k$\mathsf{\Omega }$
TB resistance ultra-scaled STT-MTJ	Value
Resistance parallel ($R_P$)	$4.1\times {10}^5$ k$\mathsf{\Omega }$
Resistance anti-parallel ($R_{AP}$)	$7.5\times {10}^5$ k$\mathsf{\Omega }$
TB resistance SOT-assisted STT-MTJ	Value
Resistance parallel ($R_P$)	$1.4\times {10}^4$ k$\mathsf{\Omega }$
Resistance anti-parallel ($R_{AP}$)	$4.2\times {10}^4$ k$\mathsf{\Omega }$

. They are consistent with CoFeB and MgO for the FM layers and TB layers, respectively. The low values of ${\lambda }_J$ and ${\lambda }_{\phi }$ are employed to have complete absorption of the transverse spin accumulation components near the TB interface [46]. The results reported in this paper were obtained by employing tetrahedral elements.

4.1. Double RL STT-MRAM

In order to reduce the critical current required for switching, an additional RL (RL${}_2$) can be deposited on top of the FL [20] (cf. Figure 1b). When RL${}_2$ is anti-parallel to the first RL (RL${}_1$), the torque coming from the two becomes additive, and the switching is made faster [62]. To not compromise the TMR and data read, the second RL is separated from the FL by a non-magnetic metallic spacer (NMS).

We employed the presented solver to perform an AP to P switching simulation of both a regular MTJ with single RL (SMTJ) and the double RL MTJ (DSMTJ). The structure used for the DSMTJ simulation, with a diameter of 40 nm, is reported in Figure 4a. Long NM contacts were employed to allow the spin accumulation to completely decay inside them. The total number of nodes in the mesh was 7338. The magnetization reversal of the FL in both the SMTJ and DSMTJ is reported in Figure 4. A voltage of 1.0 V was applied. We note that the switching of the DSMTJ, presenting more oscillations in the z-component, is less smooth than the one of the SMTJ. This is due to the additional torque and stray field contributions from the second RL, which cause the magnetization to switch less uniformly and to produce the observed non-smooth trajectory of its average components. The results show a substantial reduction in switching time for the same applied voltage in the DSMTJ, in good agreement with the experimental results reported in [20].

4.2. Ultra-Scaled STT-MRAM

The stability of the FL can be increased by adding additional MgO tunneling layers, because of the perpendicular anisotropy provided by their interfaces with CoFeB. Moreover, employing elongated layers with small diameters allows additional stability to be gained from the contribution of the shape anisotropy [21] (cf. Figure 1c). Because of the reduced FL diameter, the scalability of this kind of device is also improved.

We employed the FE solver to investigate the switching behavior of such ultra-scaled MRAM cells. The structure used for the simulation is reported in Figure 5a. The cell had a diameter of 2.3 nm, and the total number of nodes was 9634. The FL was capped by a second TB, and further split into two sections, FL${}_1$ and FL${}_2$, by a third TB. The applied bias voltage was 1.5 V. The magnetization reversal from AP to P computed in this mesh is reported in Figure 5b, where clear steps in the trajectory can be observed. This is caused by the fact that while the static magnetic coupling between the two segments increases the overall stability of the FL and allows them to respond coherently to an applied external field, during STT switching the torque contributions coming from the different TBs make the segments switch one at a time. At the beginning of the process, the torques acting from RL and FL${}_2$ on FL${}_1$ are additive, causing it to switch first and fast. Then, the torque acting from FL${}_1$ on FL${}_2$ makes it switch as well, at a slower pace. The observed behavior can help explain the reduction in critical switching current observed for a quad-interface device in [63].

4.3. SOT Assisted STT-MRAM

Including the SHE in the model allows for the proper treatment of SOT. By interfacing the FL with an HM layer, and running a second current through it, it is possible to assist the STT switching by bringing the magnetization in-plane with SOTs in the starting phase [64] (cf. Figure 1d). With SOT, the incubation time needed for the FL to break its colinearity with the RL is avoided, reducing the overall switching time at the cost of a larger footprint.

We applied the presented solver to the switching simulation of an MRAM cell relying on both STT and SOT. The structure used for the simulation is reported in Figure 6a. The diameter of the MTJ is 40 nm, and the total number of nodes is 15,058. A bias voltage of 2.0 V was used for the STT current and one of 0.3 V for the SOT current. The parameters employed for the HM are consistent with Pt, with a spin Hall angle of 0.19 [65]. The SOT current is only applied for the first 0.2 ns of the simulation. As expected, the magnetization is quickly brought in-plane by the SOT contribution. When the SOT current is turned off, the switching is completed by the STT contribution.

5. Conclusions

We presented the derivation of a finite element solution to the weak formulation of the LLG equation coupled with the spin and charge drift-diffusion formalism. The treatment of the tunneling layers as poor conductors, whose local conductivity depends on the relative magnetization orientation in the ferromagnetic layers, and the addition of appropriate boundary conditions at the tunnel barriers’ interfaces, can account for the properties of both resistance and torque expected in MTJs. The addition of terms accounting for the SHE to the spin equation allows us to also reproduce the contribution of spin-orbit torques. The demagnetizing field is computed by employing a hybrid FEM-BEM approach. The presented solver was successfully used to perform switching simulations of recently proposed structures composed of several ferromagnetic, tunneling and non-magnetic layers, as well as heavy metal lines for the generation of spin-orbit torques, supporting its employment to help investigate and predict the switching performance of newly introduced devices.