Numerical Analysis and Innovative Simulation
Techniques for Designing Advanced MRAM

4.3 Finite Element Micromagnetic Solver

Figure 4.4b illustrates the general workflow of the micromagnetic simulation software ViennaMag  [143], a finite element (FE)-based tool for the simulation of MRAM devices. Written in C++, the software is built on the open-source FE library MFEM  [144, 145], which provides essential components for finite element simulations, including mesh handling, solvers, data structures, and result visualization. On top of this foundation, ViennaMag adds features specific to micromagnetics, such as the computation of the demagnetizing field, spin-accumulation, IEC, and the solution of the LLG equation.

The simulation begins by loading input parameters, including the computational mesh that describes the device geometry. Data structures are initialized to represent the magnetization in the ferromagnetic layers, and a method for computing the demagnetizing field is selected. Once setup is complete, the main simulation loop is executed.

In each iteration, the demagnetizing field is computed from the current magnetization.

The spin-accumulation is calculated in two steps: first, the electric potential is solved for to determine the charge-current distribution, this result is then used to solve the spin-accumulation equations. These quantities serve as inputs to the LLG equation, which is assembled and solved to update the magnetization. The system evolves over discrete time steps until the final simulation time.

The LLG equation, when coupled with the spin and charge drift-diffusion formalism, enables the modeling of magnetization dynamics in multilayer structures of arbitrary complexity. Since analytical solutions to this system of PDEs are limited to idealized cases, numerical approaches are essential. Section 4.1 presents an FDM-based implementation of the LLG equation, suitable only for thin FL layers with simple elliptical geometries. In contrast, the FEM naturally supports complex geometries and multiple material domains  [146, 147], making it the preferred method for implementing a unified solver that models charge transport, spin-accumulation, and magnetization dynamics. This functionality is realized using the C++ library MFEM. The weak formulations and discretization strategies for these coupled systems are presented in Section 4.2.

4.3.1 Efficient Algorithms for
Demagnetizing Field Computation

As discussed in Section 3.2.4, the demagnetizing field potential satisfies Poisson’s equation with open boundary conditions, requiring the problem to be solved over the entire space \(\mathbb {R}^3\). This presents a significant computational challenge, as standard FEM methods are restricted to finite domains and cannot directly handle open-boundary problems. In practice, calculating the demagnetizing field is the most computationally demanding task in micromagnetic simulations because it must account for long-range dipole interactions among magnetic moments.

Various numerical approaches have been developed to address this challenge. The truncation approach is a common technique in which a sufficiently large finite region replaces the infinite exterior domain. Alternatively, the hybrid FE-BEM allows accurate treatment of open boundaries by combining FEM in the magnetized region with boundary element method (BEM) in the exterior domain, significantly reducing computational effort while maintaining accuracy.

In the following, we present and compare these methods for computing the demagnetizing field, focusing on their accuracy, efficiency, and their ability to account for the long-range nature of magnetostatic interactions.

4.3.1.1 Truncation Approach

The most straightforward method for approximating the open boundary condition is the truncation approach. Here, the magnetized region \(\Omega _m\) is embedded in a larger, finite, non-magnetic external region \(\Omega _e\), such that the total computational domain is \(\Omega = \Omega _m \cup \Omega _e\). Since the magnetic scalar potential decays with distance outside the magnetized region, the open boundary conditions can be approximated by applying homogeneous Dirichlet boundary conditions (i.e., \(u = 0\)) on the outer surface \(\partial \Omega _e\), thereby mimicking the decay to zero at infinity.

In this setting, the demagnetizing field is computed as:

\begin{equation} \mathbf {H}_{\mathrm {dem}} = -\nabla u, \label {eq::truncation_Hdem} \end{equation}

where \(u\) is the scalar magnetic potential, obtained by solving the Poisson equation derived from Equation (3.43):

\begin{equation} \nabla \cdot (\nabla u - \mathbf {M}_s) = 0 \quad \text {in } \Omega , \end{equation}

subject to the boundary condition:

\begin{equation} u = 0 \quad \text {on } \partial \Omega _e. \end{equation}

Multiplying the equation by a test function \(v \in V_0\) and integrating over the computational domain \(\Omega = \Omega _m \cup \Omega _e\), we obtain the weak form:

\begin{equation} \int _{\Omega } \left ( \nabla u - \mathbf {M}_s \right ) \cdot \nabla v \, \mathrm {d}x = 0 \quad \forall v \in V_0. \label {eq::weak_preibp} \end{equation}

Applying integration by parts to the first term (and using the fact that \(v = 0\) on \(\partial \Omega _e\)) yields:

\begin{equation} \int _{\Omega } \nabla u \cdot \nabla v \, \mathrm {d}x = \int _{\Omega } \mathbf {M}_s \cdot \nabla v \, \mathrm {d}x \quad \forall v \in V_0. \label {eq::weak_truncation} \end{equation}

To accurately account for the discontinuity of the magnetization \(\mathbf {m}\) at the boundary \(\partial \Omega _m\), the right-hand side of Equation (4.63) is typically restricted to integration over the magnetic region \(\Omega _m\):

\begin{equation} \int _{\Omega } \nabla u \cdot \nabla v \, \mathrm {d}x = \int _{\Omega _m} \mathbf {M}_s \cdot \nabla v \, \mathrm {d}x \quad \forall v \in V_0. \end{equation}

This restriction is significant when the magnetization \(\mathbf {m}\) is discretized using piecewise-affine, globally continuous finite-element functions, as it enables a consistent treatment of the discontinuity at \(\partial \Omega _m\). The function space \(V_0\) is defined as:

\begin{equation} V_0 = \{ v \in V : v(x) = 0 \ \forall x \in \partial \Omega _e \}. \end{equation}

Imposing a zero potential at a finite boundary is inherently approximate, since the true potential vanishes only in the limit \(r \to \infty \). Thus, the size of the external domain \(\Omega _e\) significantly influences the accuracy of the solution. In practice, choosing \(\Omega _e\) to be approximately five times larger than \(\Omega _m\) in each spatial direction offers a good compromise between computational efficiency and accuracy  [148] (see Figure 4.5). To further reduce computational effort without sacrificing precision, the mesh in \(\Omega _e\) can be progressively coarsened toward the boundary \(\partial \Omega _e\).

4.3.1.2 Hybrid Finite Element–Boundary Element Method

To reduce the number of DOFs in the computation of the demagnetizing field, it is advantageous to restrict spatial discretization to the magnetic domain \(\Omega _m\).
This is achieved by a hybrid FE-BEM, initially proposed by Fredkin and Koehler  [149]. Within this hybrid approach, the Poisson equation is solved in the magnetic domain using FEM, while the behavior of the potential at infinity is guaranteed by a double-layer potential placed on the boundary of the domain. Similar to the truncation method, the FE-BEM approach uses the scalar magnetic potential \(u\) as an intermediate quantity, from which the demagnetizing field is obtained via:

\begin{equation} \mathbf {H}_{\mathrm {dem}} = -\nabla u. \end{equation}

The potential is split as:

\begin{equation} u = u_1 + u_2, \end{equation}

where \(u_1\) is defined as the solution of Poisson’s equation inside the magnetic domain:

\begin{align} \nabla ^2 u_1(\mathbf {r}) & = \nabla \cdot \mathbf {M} & & \text {for } \mathbf {r} \in \Omega _m, \\ u_1 & = 0 & & \text {for } \mathbf {r} \notin \Omega _m, \end{align} with the Neumann boundary condition on \(\partial \Omega _m\):

\begin{equation} \frac {\partial u_1}{\partial n} = \mathbf {M} \cdot \mathbf {n}. \end{equation}

While \(u_1\) satisfies the correct normal jump condition at the boundary, it does not ensure continuity of the potential. Therefore, a correction term \(u_2\) is introduced such that:

\begin{align} \nabla ^2 u_2 & = 0 & & \text {in } \mathbb {R}^3, \\ \left [ \nabla u_2 \cdot \mathbf {n} \right ] & = 0 & & \text {on } \partial \Omega _m, \\ \left [ u_2 \right ] & = u_1 & & \text {on } \partial \Omega _m, \end{align} where \([f] = f_{\mathrm {in}} - f_{\mathrm {out}}\) denotes the jump of \(f\) across the boundary.

A double-layer potential of the form fulfils these conditions:

\begin{equation} u_2(\mathbf {r}) = \frac {1}{4\pi } \oint _{\partial \Omega _m} \frac {\partial u_1(\mathbf {r}')}{\partial n'} \cdot \frac {1}{|\mathbf {r} - \mathbf {r}'|} \, \mathrm {d}S'. \end{equation}

Although \(u_2\) could theoretically be evaluated in the entire domain via this integral, doing so is computationally expensive due to its \(\mathcal {O}(n^2)\) complexity. Instead, the integral is only evaluated on the boundary \(\partial \Omega _m\) using the BEM, yielding the boundary values of \(u_2\) through a matrix-vector product:

\begin{equation} \mathbf {u}_{2}\big |_{\partial \Omega _m} = \mathbf {B} \, \mathbf {u}_{1}\big |_{\partial \Omega _m}. \end{equation}

These boundary values then serve as Dirichlet boundary conditions for solving the Laplace equation \(\nabla ^2 u_2 = 0\) inside \(\Omega _m\) via the FEM, which provides \(u_2\) throughout the magnetic domain. The total potential \(u = u_1 + u_2\) is then used to compute the demagnetizing field.

The resulting matrix \(\mathbf {B}\) is dense and scales quadratically with the number of boundary nodes \(n_B\). Hierarchical matrix compression techniques  [150, 151] are often employed to mitigate the associated memory and computational costs, reducing complexity to \(\mathcal {O}(n_B \log n_B)\).

This hybrid approach significantly reduces the computational domain by discretizing only the magnetic region \(\Omega _m\). It remains effective even for disconnected or complex magnetic geometries. The trade-off, however, is the use of dense BEM operators, which differ from the sparse matrices typically encountered in FEM.

4.3.1.3 Verification and Benchmarking

The Micromagnetic Modeling Activity Group² at the NIST Center for Theoretical and Computational Materials Science has defined a set of standard problems to enable consistent comparison between micromagnetic simulation codes and published results. Following these recommendations, we assess the correctness and accuracy of our demagnetizing-field solvers using magnetization states from the \(\mu \)MAG Standard Problem 3  [152] together with geometries that admit analytical solutions.

The benchmark consists of three test cases:

• (i) Single-domain tests: A unit cube with three magnetization states (uniform, flower, vortex) to isolate solver accuracy for single-domain ferromagnets.

• (ii) External-field test: A two-disc configuration to evaluate the field outside the magnetized body.

• (iii) Device application: A realistic STT-MRAM cell with disconnected magnetic domains.

For a uniformly magnetized cube, the demagnetizing field and thus the demagnetizing energy can be obtained in closed form. The energy is given by:

\begin{equation} \mathcal {E}_d = -\frac {\mu _0 M_S}{2} \int _\Omega \mathbf {m}(\mathbf {r}) \cdot \mathbf {H}_d(\mathbf {r}) \,\mathrm {d}\mathbf {r}. \end{equation}

The scalar magnetic potential reaches its maximum near the cube surfaces, where the divergence of \(\mathbf {M}\) is largest.

The flower and vortex initializations follow  [153]. The flower state is defined as:

\begin{align} m_x(\mathbf {r}) & = \tfrac {1}{a}\,xz, \\ m_y(\mathbf {r}) & = \tfrac {1}{c}\,yz + \tfrac {1}{b^3}\,y^3 z^3, \\ m_z(\mathbf {r}) & = 1, \end{align} with parameters \(a = c = 1\) and \(b = 2\), and the cube centered at \((0,0,0)\).

The vortex state is defined as:

\begin{align} m_x(\mathbf {r}) & = -\frac {y}{r}\sqrt {1 - \exp \left (-\frac {4r^2}{r_C^2}\right )}, \\ m_y(\mathbf {r}) & = \frac {x}{r}\sqrt {1 - \exp \left (-\frac {4r^2}{r_C^2}\right )}, \\ m_z(\mathbf {r}) & = \exp \left (-\frac {2r^2}{r_C^2}\right ), \end{align} where \(r = \sqrt {x^2 + y^2}\) and \(r_C = 0.14\) as defined in  [153].

Figure 4.6 illustrates the three magnetization configurations, uniform, flower, and vortex, together with their corresponding magnetic scalar potentials. Figure 4.7 presents the demagnetizing energies computed using both the truncation and hybrid FE-BEM approaches, compared against the analytical reference values from  [153]. For quantitative comparison, we define the relative error as:

\begin{equation} \mathrm {err}_{\mathcal {E}}\,[\%] = 100 \,\frac {|\mathcal {E}_d^{\mathrm {ref}} - \mathcal {E}_d^{\mathrm {num}}|}{\max (\mathcal {E}_d^{\mathrm {ref}}) + |\min (\mathcal {E}_d^{\mathrm {ref}})|}. \label {eq::relerr_energy} \end{equation}

The hybrid FE-BEM agrees within \(<1\%\) for the uniform and flower states and approximately \(5\%\) for the vortex case at identical mesh resolution, whereas the truncation approach consistently exhibits errors greater than \(10\%\). These results confirm that, for a similar mesh in the magnetic domain, the hybrid FE-BEM approach is superior to the truncation approach in both accuracy and computational cost.

To evaluate the demagnetizing field outside the magnetized body, we consider a two-disc configuration with known analytical solution  [155]. The on-axis field of a single uniformly magnetized disc of radius \(R\) and thickness \(L\), magnetized along its symmetry axis, is given by:

\begin{equation} H_{d,x}^{\mathrm {disc}}(x) = \frac {M_S}{2}\left ( \frac {x}{\sqrt {x^2 + R^2}} - \frac {x - L}{\sqrt {(x - L)^2 + R^2}} \right ), \label {eq::disc_analytic} \end{equation}

which follows from the surface-charge representation \(\sigma _m = \mathbf {M} \cdot \mathbf {n}\).

The analytical reference for two identical discs is assembled pointwise from Equation (4.84) by mirroring about the center plane. For axial sampling points \(\{x_j\}_{j=0}^{N-1}\):

\begin{align} H_{d,x}^{\mathrm {P}}[j] & = H_1[j] + H_1[N-1-j], \\ H_{d,x}^{\mathrm {AP}}[j] & = H_1[j] - H_1[N-1-j], \end{align} where \(H_1[j] \equiv H_{d,x}^{\mathrm {disc}}(x_j)\), and P/AP denote parallel and anti-parallel magnetizations, respectively.

To quantify accuracy along the axis, we report the pointwise relative error:

\begin{equation} \mathrm {err}_H[j]\,[\%] = 100\,\frac {\bigl |H_{d,x}^{\mathrm {num}}(x_j) - H_{d,x}^{\mathrm {ref}}(x_j)\bigr |} {\max _k \bigl |H_{d,x}^{\mathrm {ref}}(x_k)\bigr |}, \label {eq::relerr_field} \end{equation}

which avoids divisions by small numbers near zero crossings (relevant in the anti-parallel case).

For the numerical setup we use discs of thickness \(L = 0.1\,\si {\nano \meter }\) and radius \(R = 1\,\si {\nano \meter }\). To obtain field values between the discs with the hybrid FE-BEM (which solves only on the magnetic domain), we insert a thin probe disc between them with negligible saturation magnetization so as not to perturb the solution. Figure 4.8 and Figure 4.9 show the field profiles and relative errors for both configurations. In both cases, the FE-BEM solution tracks the analytical reference within a few percent across the axis. At the same time, the truncation approach shows larger deviations of up to \(10\)–\(15\%\), especially near the disc faces, where gradients are strongest.

Having established accuracy on canonical test cases, we evaluate the demagnetizing field in a realistic STT-MRAM cell. The device has a diameter of \(40\,\si {\nano \meter }\) and consists of two non-magnetic contacts of length \(50\,\si {\nano \meter }\), a ferromagnetic RL of thickness \(1\,\si {\nano \meter }\), a tunnel barrier of thickness \(1\,\si {\nano \meter }\), and a free ferromagnetic layer of thickness \(1.7\,\si {\nano \meter }\).

The magnetic domain is disconnected, consisting of two ferromagnetic layers separated by the tunnel barrier. This geometry provides a thorough test for the solver. Figure 4.10 and Figure 4.11 compare the hybrid FE-BEM and truncation solutions for the P and AP magnetization configurations, respectively. For visualization, the contacts and barrier are hidden, arrows indicate the direction of the demagnetizing field, and the color map encodes the magnetic potential.

Both methods produce very similar fields and potentials, confirming that the FE-BEM implementation remains accurate for disconnected magnetic domains. Given its superior accuracy at comparable mesh resolutions and reduced computational effort relative to large truncated domains, the hybrid FE-BEM is adopted for the complex geometries considered in subsequent chapters, including synthetic antiferromagnets.