Numerical Analysis and Innovative Simulation
Techniques for Designing Advanced MRAM

Chapter 4 Computational Methods

The micromagnetic model introduced in the preceding sections defines a set of nonlinear partial differential equations in space and time that are analytically tractable only for highly simplified cases. In practice, the simulation of static and dynamic micromagnetic phenomena relies on numerical techniques. However, developing efficient numerical methods is nontrivial due to several intrinsic challenges of the underlying equations.

First, the demagnetization field, which represents long-range dipole-dipole interactions, introduces a computational bottleneck. A simple evaluation scales as \(\mathcal {O}(n^2)\) quickly becomes restrictive with the number of simulation cells \(n\). To address this, several approaches have been proposed that reduce the complexity to \(\mathcal {O}(n \log n)\) [126] or even \(\mathcal {O}(n)\) [127].

Second, the exchange interaction imposes a local but highly stiff coupling due to its second-order spatial derivative. This stiffness, combined with the nonlinearity of the governing equations, places strict requirements on the numerical time integration scheme.

Third, the total micromagnetic energy landscape is typically nonconvex and high-dimensional, with multiple local minima. In quasistatic problems such as hysteresis simulations, the system must relax toward a nearby energy minimum. The complexity of the energy landscape necessitates careful selection of minimization algorithms to avoid missing relevant minima.

To solve the micromagnetic equations numerically, spatial and temporal discretizations are typically decoupled. For spatial discretization, the most widely used methods are the finite difference method (FDM) and the finite element method (FEM). Both approaches partition the magnetic domain into cells to compute the local contributions to the effective field \(\mathbf {H}_\mathrm {eff}\) or the total energy \(E\). While FDM is typically based on a regular Cartesian grid, FEM uses tetrahedral elements to accommodate complex geometries (see Figure 4.1).

The mesh resolution in either method must be fine enough to resolve domain walls, whose characteristic width is given by the exchange length:

\begin{equation} \ell = \sqrt {\frac {A}{K_\mathrm {eff}}}, \end{equation}

where \(A\) is the exchange stiffness and \(K_\mathrm {eff}\) denotes the effective anisotropy constant. The latter includes crystalline anisotropy and shape anisotropy contributions induced by the demagnetization field  [99].

Once the system is spatially discretized, the static energy minimum or the dynamic time evolution of the magnetization \(\mathbf {m}\) must be computed. This requires a separate class of algorithms, such as energy-minimization techniques or time integration methods for the LLG equation. The effective field must be evaluated only within the magnetic domain in both cases.

The following sections review two spatial discretization techniques: the FDM (Section 4.1) and the FEM (Section 4.2). Particular attention is given to efficient algorithms for computing the demagnetizing field, including the truncation approach and the hybrid finite element-boundary element method (FE-BEM), which are presented in Section 4.3.1. The numerical time integration of the LLG equation is treated separately in Chapter 5.

4.1 Finite Difference Method

The FDM is a widely used numerical technique for solving differential equations and has become a cornerstone of computational science due to its conceptual simplicity and ease of implementation. At its core, FDM seeks to approximate solutions to differential equations, i.e., to find a function (or a discrete approximation) that satisfies the prescribed relationships between various derivatives on a given spatial and/or temporal domain, subject to boundary conditions. This is a nontrivial problem, and analytical solutions are rarely obtained. The finite difference method addresses this challenge by replacing derivatives with finite difference approximations. This approximation yields a large but finite system of algebraic equations that can be solved numerically on a computer [128, 129, 130].

Before solving such systems, one must first understand how to approximate derivatives of a known function using only discrete pointwise evaluations. This foundational step enables the development of higher-order schemes and introduces essential concepts, including the consistency, stability, and accuracy of finite-difference operators. The discretization typically involves subdividing the spatial domain into an equidistant grid of points. In the 1D case, this corresponds to grid points:

\begin{equation} x_i = ih, \quad \text {for } i = 0, \ldots , N, \end{equation}

where \(h = \frac {1}{N}\) is the uniform grid spacing, and \(N \in \mathbb {N}\) is the number of discretization points. This concept naturally extends to higher dimensions. For example, in a 2D setting, the grid points are given by:

\begin{equation} \mathbf {x}_{i,j} = (ih, jh), \quad \text {for } i,j = 0, \ldots , N, \end{equation}

with the same spacing \(h = \frac {1}{N}\) applied along both axes.

Figure 4.2 illustrates the grid layout in both the 1D and 2D cases. In 1D (Figure 4.2 (a)), the points lie along a line with uniform spacing. In 2D (Figure 4.2 (b)), the points form a Cartesian mesh.

In micromagnetics, FDM has proven to be a robust and efficient tool for simulating magnetization dynamics. Its restriction to regular, structured meshes enables the use of highly optimized algorithms, particularly for computing quantities like the demagnetization field, which are computationally expensive due to their long-range nature. However, this same grid regularity limits the applicability of FDM to geometries that conform well to Cartesian meshes, making FDM less suitable for complex or curved domains, where unstructured methods such as FEM are more appropriate.

FDM is a classical numerical technique for approximating the solution of differential equations by discretizing space into a uniform grid and replacing derivatives with finite difference quotients. If one can efficiently and accurately approximate derivatives at discrete points, then in principle, one can also approximate the solution of differential equations. This process begins by approximating derivatives of a real-valued function \(u \colon [a,b] \to \mathbb {R}\) using difference formulas:

• (i) Forward difference:

  \(\seteqnumber{0}{4.}{3}\)

  \begin{equation} \frac {u(x+h) - u(x)}{h} \end{equation}

• (ii) Backward difference:

  \(\seteqnumber{0}{4.}{4}\)

  \begin{equation} \frac {u(x) - u(x-h)}{h} \end{equation}

• (iii) Central difference:

  \(\seteqnumber{0}{4.}{5}\)

  \begin{equation} \frac {u(x+h) - u(x-h)}{2h} \end{equation}

The central difference formula is generally more accurate, introducing an error term of order \(\mathcal {O}(h^2)\), while forward and backward differences have errors of order \(\mathcal {O}(h)\). These formulas can be derived from Taylor expansions or from interpolating polynomials evaluated at the mesh points.

This idea extends to second-order derivatives. In one spatial dimension, the second derivative of a function \(u(x)\) at a point \(x_i\) is approximated by applying the central difference formula twice, once for the first derivative and again for the derivative of that result. This yields the second-order central difference approximation:

\begin{equation} \frac {\mathrm {d}^2 u}{\mathrm {d}x^2}\Big |_{x_i} \approx \frac {u_{i-1} - 2u_i + u_{i+1}}{h^2} + \mathcal {O}(h^2), \end{equation}

where \(h\) is the grid spacing and \(u_i\) denotes the discrete function value at grid point \(x_i\), i.e. \(u_i \approx u(x_i)\). This approximation introduces a truncation error of order \(\mathcal {O}(h^2)\), making it simple and reasonably accurate for many applications.

To illustrate how FDM works in practice, we consider the Poisson equation, a fundamental equation arising in many scientific and engineering applications:

\begin{align} - \nabla ^2 u(x) & = f(x), \quad x \in \Omega , \\ u(x) & = u_D(x), \quad x \in \partial \Omega . \end{align}

In one dimension, the domain \([0,1]\) is discretized into a uniform mesh of \(N \in \mathbb {N}\) grid points, with spacing \(h = \frac {1}{N}\). Using the finite difference approximation, the continuous Poisson equation becomes:

\begin{equation} - \frac {u_{i-1} - 2u_i + u_{i+1}}{h^2} = f_i, \quad i = 1, \ldots , N-1. \end{equation}

This leads to a linear system of equations:

\begin{equation} A u = r, \end{equation}

where \(A\) is a sparse matrix representing the discrete Laplace operator, \(u\) is the vector of unknowns, and \(r\) is the discrete right-hand side vector.

Dirichlet boundary conditions can be directly incorporated by assigning values to the boundary nodes \(u_0\) and \(u_N\). To evaluate the difference quotient at the domain boundaries, one might use ghost points, with values such as:

\begin{equation} u_{-1} = u_1, \quad u_{N+1} = u_{N-1}. \end{equation}

This ensures the scheme can be applied across the full grid. Once the system is assembled, the solution can be computed as:

\begin{equation} u = A^{-1} r, \end{equation}

in practice, iterative solvers such as Jacobi, Gauss-Seidel, or successive over-relaxation (SOR) are typically used due to the size and sparsity of \(A\) [131, 132].

To extend finite difference methods to higher dimensions, consider the Poisson equation in two dimensions:

\begin{equation} - \nabla ^2 u(x, y) = -\left ( \frac {\partial ^2 u}{\partial x^2} + \frac {\partial ^2 u}{\partial y^2} \right ) = f(x, y), \quad (x, y) \in \Omega = (0,1) \times (0,1), \end{equation}

\begin{equation} u(x, y) = g(x, y), \quad (x, y) \in \partial \Omega . \end{equation}

Here, \(g(x, y)\) defines the Dirichlet boundary condition on the domain boundary. As in the 1D case, we discretize the domain by defining a mesh that includes the boundary:

\begin{equation} \mathcal {M}_{M,N} = \left \{ (x_i, y_j) = \left ( \frac {i}{M}, \frac {j}{N} \right ) \,\middle |\, 0 \leq i \leq M, \ 0 \leq j \leq N \right \}, \end{equation}

where \(M, N \in \mathbb {N}\) are the number of intervals along each axis.

To discretize the differential operator, we approximate the second derivatives in each spatial direction using central finite differences:

\begin{align} \frac {\partial ^2 u}{\partial x^2}(x_i, y_j) & \approx \frac {u_{i+1,j} - 2u_{i,j} + u_{i-1,j}}{h^2}, \\ \frac {\partial ^2 u}{\partial y^2}(x_i, y_j) & \approx \frac {u_{i,j+1} - 2u_{i,j} + u_{i,j-1}}{h^2}. \end{align} By adding these two approximations, the full Laplacian is approximated as:

\begin{equation} \left ( \frac {\partial ^2 u}{\partial x^2} + \frac {\partial ^2 u}{\partial y^2} \right )(x_i, y_j) \approx \frac {u_{i+1,j} + u_{i-1,j} + u_{i,j+1} + u_{i,j-1} - 4u_{i,j}}{h^2}. \end{equation}

Substituting this into the 2D Poisson equation results in a discrete equation for interior points:

\begin{equation} - \frac {u_{i+1,j} + u_{i-1,j} + u_{i,j+1} + u_{i,j-1} - 4u_{i,j}}{h^2} = f_{i,j}. \end{equation}

The boundary values are set using the function \(g(x, y)\), just like in the 1D case. The resulting linear system again has a sparse matrix structure and can be solved using iterative solvers designed for large-scale problems.

The extension of this method to three dimensions is conceptually straightforward. The 3D Laplacian is defined as:

\begin{equation} \nabla ^2 u(x, y, z) = \frac {\partial ^2 u}{\partial x^2} + \frac {\partial ^2 u}{\partial y^2} + \frac {\partial ^2 u}{\partial z^2}. \end{equation}

It can be discretized by independently applying the same finite difference approximations in all three spatial directions. This yields a discrete equation involving the values of \(u\) at neighboring grid points along the \(x\), \(y\), and \(z\) directions. The resulting system of equations has a similar sparse structure and is solved with the same numerical techniques used in lower dimensions.

4.1.1 Discretization of the Landau-Lifshitz-Gilbert Equation

The simulations in this work are conducted using an in-house micromagnetic solver based on the FDM, originally introduced in  [103]. It is applied to the modeling of pMTJ devices. The effective magnetic field is computed directly from the discretized magnetization  [133]. Only the dynamics of the FL are simulated, while the magnetostatic stray field of the RL is treated as a static external contribution.

The computational domain is discretized into parallelepiped cells with volume \(\Delta x \cdot \Delta y \cdot \Delta z\), where \(\Delta x\), \(\Delta y\), and \(\Delta z\) represent the spatial resolution in each Cartesian direction. The number of cells in each direction is denoted by \(N_x\), \(N_y\), and \(N_z\), respectively. A unit magnetization vector \(\mathbf {m}\) is assigned to the center of each cell.

4.1.1.1 External Field

The external magnetic field is specified as a uniform, time-independent input:

\begin{equation} \mathbf {H}_\text {ext}(i, j, k) = \mathbf {H}_\text {ext,input} \end{equation}

4.1.1.2 Exchange Field

The exchange field contribution, being proportional to the Laplacian of the magnetization as shown in Equation (3.27), requires knowledge of the magnetization in the computational cell and its six nearest neighbors, two along each Cartesian direction. In a regular FDM grid, these neighbors are located at the adjacent face, sharing positions, specifically at coordinates \((i\pm 1,j,k)\), \((i,j\pm 1,k)\), and \((i,j,k\pm 1)\).

The first-order discretization of the exchange field at cell \((i,j,k)\) yields the following expression  [107]:

\begin{equation} \mathbf {H}_\text {ex}(i,j,k) = \frac {2A}{\mu _0 M_s} \sum _{(i',j',k')} \frac {\mathbf {m}(i',j',k') - \mathbf {m}(i,j,k)}{\Delta ^2_{i',j',k'}}, \end{equation}

where \(A\) is the exchange stiffness constant, \(\mu _0\) the vacuum permeability, \(M_s\) the saturation magnetization, and \(\Delta _{i',j',k'} \in \{\Delta x, \Delta y, \Delta z\}\) denotes the grid spacing in the direction connecting the central and neighboring cells.

At the boundaries of the computational domain, some neighboring positions may lie outside the simulation grid. To preserve the consistency of the exchange field calculation near these boundaries, so-called ghost cells were introduced. These are virtual cells placed just outside the domain, with their magnetization vectors set equal to the nearest valid interior cell, thereby enforcing Neumann-like boundary conditions as detailed in Section 4.1. This approach ensures that the exchange field can be consistently computed throughout the domain, including at the boundaries  [129].

4.1.1.3 Anisotropy Field

The computation of the discretized anisotropy field contribution only requires the magnetization vector within the computational cell. Equations (3.33) and (3.37) can be directly applied for field evaluation. In the present work, the interface-induced perpendicular anisotropy, typical in pMTJs, is modeled as a uniaxial anisotropy contribution. The corresponding discretized expression for the anisotropy field is given by:

\begin{equation} \mathbf {H}_\text {ani}(i,j,k) = \frac {2K_1}{\mu _0 M_s} (\mathbf {m}(i,j,k) \cdot \mathbf {e}) \mathbf {e}, \end{equation}

where \(\mathbf {e}\) is the unit vector defining the easy axis of magnetization, typically oriented perpendicular to the plane of the structure in systems exhibiting PMA.

For materials characterized by cubic anisotropy, the following expression gives the field contribution:

\begin{equation} \mathbf {H}_\text {ani}(i,j,k) = -\frac {2 D(i,j,k)}{\mu _0 M_s} \mathbf {m}(i,j,k), \end{equation}

where \(D(i,j,k)\) is a second-rank anisotropy matrix that defines the local anisotropy energy landscape, the construction and role of this matrix are described in detail in Section 3.2.3.3.

4.1.1.4 Demagnetizing Field

The demagnetizing field is the most computationally expensive contribution to the effective field. Originating from dipole-dipole interactions, it is inherently non-local: its value at a given location depends on the magnetization vectors throughout the magnetic structure. Consequently, calculating the demagnetizing field at each cell requires summing contributions from all other cells in the domain.

For the computation of the demagnetizing field in the FL, a discretized form of Equation (3.46) was used. The field at each computational node is calculated as:

\begin{equation} \mathbf {H}_\text {d}(i,j,k) = \frac {M_s}{4\pi } \sum _{i'=1}^{N_x} \sum _{j'=1}^{N_y} \sum _{k'=1}^{N_z} \tilde {G}(i,j,k,i',j',k') \cdot \mathbf {m}(i',j',k'), \label {eq::demag_field} \end{equation}

where \(\tilde {G}\) is a space-dependent tensor that encodes the dipole-dipole interaction kernel  [107, 133]. This tensor is symmetric, and given its structure, only six of its nine components need to be evaluated explicitly. Since it is time-independent and solely geometry-dependent, \(\tilde {G}\) is precomputed once at the beginning of the simulation and reused in subsequent time steps.

The magnetostatic coupling field between the RL and the FL shares the exact physical origin as the demagnetizing field in the FL, arising from dipole-dipole interactions. However, because the RL magnetization remains static throughout the simulation, its contribution can be computed once during initialization and treated as a time-invariant external field.

To evaluate this coupling, the RL is discretized into parallelepiped cells, each assigned a unit magnetization vector \(\mathbf {p} = \mathbf {M}_\text {RL} / M_s^\text {RL}\), where \(\mathbf {M}_\text {RL}\) denotes the local magnetization and \(M_s^\text {RL}\) is the saturation magnetization of the RL.

The resulting magnetostatic field at a given FL cell \((i,j,k)\) is then computed using an expression analogous to Equation (4.26):

\begin{equation} \mathbf {H}_\text {ms}(i,j,k) = \frac {M_s}{4\pi } \sum _{i'=1}^{N^{RL}_x} \sum _{j'=1}^{N^{RL}_y} \sum _{k'=1}^{N^{RL}_z} \tilde {G}(i,j,k,i',j',k') \cdot \mathbf {p}(i',j',k'), \end{equation}

where \(N^{RL}_x\), \(N^{RL}_y\), and \(N^{RL}_z\) represent the number of discretized RL cells along the \(x\)-, \(y\)-, and \(z\)-directions, respectively.

By precomputing this static contribution, the magnetostatic interaction between RL and FL can be incorporated efficiently without increasing the computational burden during the dynamic simulation steps.

4.1.1.5 Ampère Field

The Ampère (or Oersted) field, resulting from the magnetic field generated by the flowing current, is computed using the discretized form of Equation (3.46). The field at each computational cell \((i,j,k)\) is obtained by summing the contributions from the current density in all other cells of the simulation domain:

\begin{equation} \mathbf {H}_\text {curr}(i,j,k) = \frac {1}{4\pi } \sum _{i'} \sum _{j'} \sum _{k'} \mathbf {J}_\text {c}(i',j',k') \times \mathbf {G}'(i,j,k,i',j',k'), \end{equation}

where \(\mathbf {J}_\text {c}\) is the local current density vector and \(\mathbf {G}'\) is a space-dependent interaction vector that encodes the geometric relationship between the source and observation cells. The summation excludes the self-term, i.e., cells for which \((i',j',k') = (i,j,k)\) are omitted.

A detailed derivation of the interaction vector \(\mathbf {G}'\), including methods for the efficient numerical evaluation of the required integrals, is provided in  [134]. In practice, the Oersted field is approximated at the start of the simulation under the assumption of a constant current-density distribution. It is then treated as a static external field during time integration.

4.1.1.6 Thermal Field

The thermal field is modeled as spatially and temporally uncorrelated Gaussian noise, consistent with the stochastic formulation of the LLG equation presented in  [135]. Its implementation must also adhere to the fluctuation-dissipation relation described by Equation (3.52b), ensuring proper thermal agitation behavior across the simulation domain.

The discretized thermal field at each computational cell \((i,j,k)\) is computed as:

\begin{equation} \mathbf {H}_\text {th}(i,j,k) = \sigma (i,j,k) \sqrt {\frac {\alpha }{1+\alpha ^2} \cdot \frac {2k_B T}{\gamma _0 \Delta t \Delta V M_s}}, \end{equation}

where \(\alpha \) is the Gilbert damping factor, \(k_B\) the Boltzmann constant, \(T\) the temperature, \(\gamma _0\) the gyromagnetic ratio, \(\Delta t\) the simulation time step, \(\Delta V\) the volume of the computational cell, and \(M_s\) the saturation magnetization. The term \(\sigma (i,j,k)\) represents a random variable drawn from a standard normal distribution, uncorrelated in space and time, with zero mean and unit standard deviation.

The implementation follows the approach in  [136], in which each component of the thermal field vector is independently generated from Gaussian noise. This ensures that the thermal fluctuations correctly reflect the statistical properties of a system in thermal equilibrium and preserve the physical fidelity of thermally activated switching processes.

4.1.1.7 Time Integration Scheme

After the initial setup phase, where simulation parameters are loaded from either a configuration file or the command line, the magnetostatic coupling field \(\mathbf {H}_\text {ms}\) and the Ampère field \(\mathbf {H}_\text {curr}\) are calculated. Since these contributions remain constant throughout the simulation, they are precomputed and stored as static fields.

The core of the time evolution is handled by an explicit fourth-order Runge-Kutta (RK4) integration scheme, used to solve the LLG equation. Within each time step, the magnetization update proceeds via the evaluation of four successive intermediate stages:

\begin{align} \mathbf {k}_1 & = \text {LLG}(\mathbf {m}_n, t_n)\Delta t, \\ \mathbf {k}_2 & = \text {LLG}\left (\mathbf {m}_n + \frac {\mathbf {k}_1}{2}, t_n + \frac {\Delta t}{2}\right )\Delta t, \\ \mathbf {k}_3 & = \text {LLG}\left (\mathbf {m}_n + \frac {\mathbf {k}_2}{2}, t_n + \frac {\Delta t}{2}\right )\Delta t, \\ \mathbf {k}_4 & = \text {LLG}(\mathbf {m}_n + \mathbf {k}_3, t_n + \Delta t)\Delta t. \end{align} Using these coefficients, the magnetization is updated as:

\begin{equation} \mathbf {m}_{n+1} = \mathbf {m}_n + \frac {\Delta t}{6}(\mathbf {k}_1 + 2\mathbf {k}_2 + 2\mathbf {k}_3 + \mathbf {k}_4). \end{equation}

This method balances accuracy and computational efficiency, making it well-suited for simulating magnetization dynamics over long time scales. The time loop continues until the desired total simulation time is reached, with field contributions and magnetization consistently updated at each step.

Further emphasis on the mathematical formulation, stability considerations, and numerical performance of time integration methods will be provided in Chapter 5.

Figure 4.4a illustrates the overall flow of this simulation framework.