Numerical Analysis and Innovative Simulation
Techniques for Designing Advanced MRAM

3.2 Effective Magnetic Field

The total energy of a ferromagnetic system with respect to its magnetization consists of several contributions that depend on the material’s intrinsic properties. Some of these, such as the demagnetization (or magnetostatic) energy, can be described within the framework of classical magnetostatics. Others, such as exchange energy and magnetocrystalline anisotropy energy, originate from quantum-mechanical interactions.

To derive expressions for the individual contributions to the effective magnetic field used in the LLG equation, one must first define the magnetic free Gibbs energy \(\mathcal {E}\) of the ferromagnetic domain. The energy functional comprises several contributions, most notably the exchange energy, the anisotropy energy, the Zeeman energy from an externally applied magnetic field, and the magnetostatic energy. These components will be introduced in more detail in the following sections.

Once the total free Gibbs energy is defined, the equilibrium state of the magnetization corresponds to a configuration that minimizes this energy, leading to the variational problem:

\begin{equation} \min _{|\mathbf {m}| = 1} \mathcal {E}(\mathbf {m}) \label {eq::energy_minimization} \end{equation}

A stable magnetization configuration thus represents a balance among competing energy terms. The effective magnetic field \(\mathbf {H}_{\mathrm {eff}}\), which governs the magnetization dynamics in the LLG equation, is defined as the functional derivative of the energy with respect to the magnetization:

\begin{equation} \mathbf {H}_{\mathrm {eff}} = -\frac {1}{\mu _0 M_s} \frac {\delta \mathcal {E}}{\delta \mathbf {m}} \label {eq::heff_definition} \end{equation}

From the Euler-Lagrange equations associated with the minimization in Equation (3.16), one obtains the boundary value problem:

\begin{align} \mathbf {H}_{\mathrm {eff}} & = 0, \label {eq::heff_zero} \\ \frac {\delta \mathcal {E}(\mathbf {m})}{\delta \nabla m_i} \cdot n & = 0, \quad \text {for } i = 1, 2, 3 \label {eq::variational_derivative} \end{align} where \(\mathbf {n}\) denotes the unit normal vector at the boundary of the ferromagnetic domain. These conditions imply that, at equilibrium, the magnetization vector is locally aligned with the effective field. Consequently, the LLG equation describes a relaxation process toward a state of minimal magnetic Gibbs energy.

In this thesis, the effective field is composed of the following contributions:

\begin{equation} \mathbf {H}_{\mathrm {eff}} = \mathbf {H}_{\mathrm {exc}} + \mathbf {H}_{\mathrm {ani}} + \mathbf {H}_{\mathrm {ext}} + \mathbf {H}_{\mathrm {demag}} + \mathbf {H}_{\mathrm {curr}} + \mathbf {H}_{\mathrm {th}} + \mathbf {H}_{\mathrm {iec}} \label {eq::heff_components} \end{equation}

Here, \(\mathbf {H}_{\mathrm {exc}}\) is the exchange field, \(\mathbf {H}_{\mathrm {ani}}\) the anisotropy field, \(\mathbf {H}_{\mathrm {ext}}\) the externally applied field, \(\mathbf {H}_{\mathrm {demag}}\) the demagnetizing (magnetostatic) field, \(\mathbf {H}_{\mathrm {curr}}\) the ampère field, \(\mathbf {H}_{\mathrm {th}}\) the thermal field, \(\mathbf {H}_{\mathrm {iec}}\) the interlayer exchange field. Each of these contributions will be discussed in the subsequent sections.

3.2.1 Exchange Field

A defining property of ferromagnetic materials is their ability to retain a nonzero macroscopic magnetization, known as remanent magnetization, even without an external field. In contrast, systems in which magnetic moments interact only via dipole-dipole coupling do not exhibit such behavior. Their net magnetization vanishes when the external field is removed, as predicted by classical electrodynamics  [95].

This distinctly ferromagnetic behavior arises from the quantum mechanical exchange interaction. First proposed by W. Heisenberg in 1926  [96], the exchange interaction favors parallel alignment of neighboring electron spins. Its origin lies in the antisymmetrization requirement imposed on the two-electron wavefunction by the Pauli exclusion principle. Depending on the spin alignment, the symmetry of the total wavefunction affects the spatial part, thereby altering the expected Coulomb energy. A symmetric spatial distribution (corresponding to parallel spins) leads to a larger average separation and consequently a lower Coulomb repulsion. This energetic preference underpins ferromagnetism.

The classical description of this effect is given by the Heisenberg model, in which the exchange energy between two localized unit spins \(\mathbf {s}_i\) and \(\mathbf {s}_j\) is expressed as:

\begin{equation} \mathcal {E}^\text {ex}_{ij} = -J\, \mathbf {s}_i \cdot \mathbf {s}_j, \label {eq::heisenberg_exchange} \end{equation}

where \(J\) is the exchange integral, a material- and distance-dependent parameter. For a detailed derivation of the Heisenberg model, the reader is referred to standard quantum mechanics textbooks  [97].

In a discrete spin system, the total exchange energy at a spin site \(\mathbf {x}\) is given by:

\begin{equation} \mathcal {E}^\text {ex}_\mathbf {x} = -\sum _i \frac {J_i}{2} \, \mathbf {m}(\mathbf {x}) \cdot \mathbf {m}(\mathbf {x} - \boldsymbol {\xi }_i), \end{equation}

where \(\boldsymbol {\xi }_i\) denotes the displacement vector to neighboring spins and \(\mathbf {m}\) is the normalized magnetization.

Using a Taylor expansion and assuming identical coupling \(J_i = J\) and uniform spacing \(\boldsymbol {\xi }_i\), this discrete expression can be approximated as:

\begin{equation} \mathcal {E}^\text {ex}_\mathbf {x} \approx -\sum _i \frac {J}{2} \left [1 - \frac {1}{2} (\nabla \mathbf {m} \cdot \boldsymbol {\xi }_i)^2\right ] + \mathcal {O}(\boldsymbol {\xi }_i^3). \end{equation}

Neglecting constant terms and higher-order corrections, and transitioning to the continuum limit under the assumption of a regular lattice, the exchange energy becomes:

\begin{equation} \mathcal {E}^\text {ex} = \int _{\Omega } \sum _{i,j,k} A_{jk} \, \frac {\partial m_i}{\partial x_j} \frac {\partial m_i}{\partial x_k} \, \mathrm {d}\mathbf {x}, \end{equation}

where \(A_{jk}\) are the components of a symmetric positive-definite matrix that reflects the anisotropy of the exchange interaction and \(\Omega \) is the volume of the ferromagnetic domain. The constant offset from integrating the zeroth-order term is typically omitted, as it does not influence the magnetization dynamics.

The matrix \(A_{jk}\) can always be diagonalized by a suitable choice of coordinate system  [98], yielding a simplified form:

\begin{equation} \mathcal {E}^\text {ex} = \int _{\Omega } \sum _{i,j} A \left (\frac {\partial m_i}{\partial x_j}\right )^2 \, \mathrm {d}\mathbf {x}. \end{equation}

For cubic and isotropic lattice structures, the exchange coefficients further reduce to a scalar exchange stiffness \(A\), giving the widely used micromagnetic expression:

\begin{equation} \mathcal {E}^\text {ex} = A \int _{\Omega } \|\nabla \mathbf {m}(\mathbf {x})\|_F^2 \, \mathrm {d}\mathbf {x}, \label {eq::exchange_energy} \end{equation}

where \(\|\nabla \mathbf {m}\|_F^2\) is the Frobenius norm squared of the gradient tensor. Although this form was originally derived for localized spins and isotropic lattices, it has been shown to effectively describe a broad range of materials, including band magnets and anisotropic systems  [99].

Equation (3.25) represents the lowest-order phenomenological expression that penalizes spatial inhomogeneities in the magnetization vector field. It plays a central role in micromagnetic simulations. The general relation for the effective field is:

\begin{equation} \mathbf {H}_\text {eff} = -\frac {1}{\mu _0 M_s} \frac {\delta E}{\delta \mathbf {m}}, \label {eq::heff_general} \end{equation}

the exchange field contribution becomes:

\begin{equation} \label {eq::exchange_field} \mathbf {H}_\text {exc} = \frac {2A}{\mu _0 M_s} \nabla ^2 \mathbf {m}. \end{equation}

3.2.2 Zeeman Field

An external magnetic field \(\mathbf {H}_\mathrm {ext}\) couples directly to the magnetization and tends to align it along the field direction. The associated energy contribution, known as the Zeeman energy, penalizes deviations of the magnetization from the external field and is expressed as:

\begin{equation} \mathcal {E}_\mathrm {ext} = -\mu _0 M_s \int _{\Omega } \mathbf {H}_\mathrm {ext} \cdot \mathbf {m} \, \mathrm {d}\mathbf {x}. \label {eq::zeeman_energy} \end{equation}

Here, \(\mathbf {m}\) denotes the normalized magnetization vector and \(M_s\) the saturation magnetization. The resulting field contribution to the effective field follows directly from functional differentiation:

\begin{equation} \mathbf {H}_\mathrm {ext} = -\frac {1}{\mu _0 M_s}\frac {\delta \mathcal {E}_\mathrm {ext}}{\delta \mathbf {m}}. \label {eq::zeeman_field} \end{equation}

This term represents the most straightforward contribution to the effective field, as it simply reflects the action of an externally applied magnetic field on the magnetic system.

3.2.3 Anisotropy Field

Another vital contribution to the total magnetic free energy is the anisotropy energy, which favors magnetization alignment along specific directions known as easy axes. The energy required to reorient the magnetization in anisotropic materials depends on its direction relative to these preferred axes. This direction-dependent behavior arises from several sources, including:

• (i) anisotropic crystal structures (magnetocrystalline anisotropy),

• (ii) lattice deformation or spin-orbit coupling at interfaces (interface anisotropy),

• (iii) and the geometry of the magnetic body itself (shape anisotropy)  [99].

These anisotropies are often described phenomenologically by constructing energy expressions based on symmetry arguments. A key property is that the energy must be invariant under magnetization reversal:

\begin{equation} \mathcal {E}_\mathrm {ani}(\mathbf {m}) = \mathcal {E}_\mathrm {ani}(-\mathbf {m}). \end{equation}

3.2.3.1 Uniaxial Anisotropy

A single easy axis \(\mathbf {e}\) exists in the simplest case of uniaxial anisotropy. The anisotropy energy density depends only on the angle between the magnetization \(\mathbf {m}\) and the easy axis \(\mathbf {e}\). A general expansion in even powers yields:

\begin{equation} \mathcal {E}_\mathrm {ani}^{\mathrm {uni}} = \int _{\Omega } \left [ K_{u1} (\mathbf {m} \cdot \mathbf {e})^2 + K_{u2} (\mathbf {m} \cdot \mathbf {e})^4 + \mathcal {O}((\mathbf {m} \cdot \mathbf {e})^6) \right ] \, \mathrm {d}\mathbf {x}, \label {eq::uniaxial_energy} \end{equation}

where \(K_{u1}\) and \(K_{u2}\) are the first- and second-order uniaxial anisotropy constants. For many applications, the first term suffices:

\begin{equation} \mathcal {E}_\mathrm {ani}^{\mathrm {uni}} = -\int _{\Omega } K_u (\mathbf {m} \cdot \mathbf {e})^2 \, \mathrm {d}\mathbf {x}. \label {eq::uniaxial_simplified} \end{equation}

Here, \(K_u > 0\) favors alignment of \(\mathbf {m}\) along \(\pm \mathbf {e}\), while \(K_u < 0\) leads to easy-plane anisotropy, where the magnetization tends to lie perpendicular to the axis.

The corresponding anisotropy field is obtained via functional differentiation:

\begin{equation} \mathbf {H}_\mathrm {ani}^{\mathrm {uni}} = \frac {2K_1}{\mu _0 M_s} (\mathbf {m} \cdot \mathbf {e}) \mathbf {e}. \label {eq::uniaxial_field} \end{equation}

Uniaxial anisotropy is typically observed in hexagonal or tetragonal crystal structures, such as cobalt. An example energy landscape for uniaxial anisotropy is shown in Figure 3.2 (a).

3.2.3.2 Planar Anisotropy

Planar anisotropy is a special case of uniaxial anisotropy in which the magnetization prefers to lie in a plane perpendicular to a given axis \(\mathbf {e}\). This corresponds to a negative uniaxial anisotropy constant \(K_u < 0\). The anisotropy energy density takes the same mathematical form as in Equation (3.32), but with the opposite sign:

\begin{equation} \mathcal {E}_\mathrm {ani}^{\mathrm {planar}} = \int _{\Omega } K_p (\mathbf {m} \cdot \mathbf {e})^2 \, \mathrm {d}\mathbf {x}, \label {eq::planar_energy} \end{equation}

where \(K_p = -K_u > 0\) encourages magnetization to remain orthogonal to \(\mathbf {e}\).

The corresponding anisotropy field is given by:

\begin{equation} \mathbf {H}_\mathrm {ani}^{\mathrm {planar}} = -\frac {2K_p}{\mu _0 M_s} (\mathbf {m} \cdot \mathbf {e}) \mathbf {e}. \label {eq::planar_field} \end{equation}

Planar anisotropy is common in ultrathin magnetic films or multilayers, where the shape or surface anisotropy causes the magnetization to align in-plane. It is effectively modeled by taking the uniaxial formulation with \(K_u < 0\). Figure 3.2 (b) illustrates the planar anisotropy surface, where the minimum energy lies in the plane orthogonal to the easy axis.

3.2.3.3 Cubic Anisotropy

In materials with cubic lattice symmetry, such as body-centered cubic iron or face-centered cubic nickel, the magnetocrystalline anisotropy energy depends on the projections of \(\mathbf {m}\) onto three orthogonal easy axes \(\mathbf {e}_1\), \(\mathbf {e}_2\), and \(\mathbf {e}_3\), satisfying \(\mathbf {e}_i \cdot \mathbf {e}_j = \delta _{ij}\). The anisotropy energy is then expressed as:

\begin{equation} \begin{aligned} \mathcal {E}_\mathrm {ani}^{\mathrm {cub}} = \int _{\Omega } & \left [ K_{c1} \left ( m_1^2 m_2^2 + m_2^2 m_3^2 + m_3^2 m_1^2 \right ) \right . \left . + K_{c2} m_1^2 m_2^2 m_3^2 \right ] \, \mathrm {d}\mathbf {x}, \label {eq::cubic_energy} \end {aligned} \end{equation}

where \(m_i = \mathbf {m} \cdot \mathbf {e}_i\) are the projections of the magnetization along the easy axes. The coefficients \(K_{c1}\) and \(K_{c2}\) are material-dependent anisotropy constants.

As with uniaxial anisotropy, the sign of \(K_{c1}\) and \(K_{c2}\) determines whether the magnetization prefers to align with or avoid the easy axes. Figure 3.2 (c) visualizes the corresponding energy surface. For instance, negative constants in nickel (fcc structure) lead to multiple easy axes (four in total).

The effective anisotropy field for cubic symmetry can be written in matrix form:

\begin{equation} \mathbf {H}_\mathrm {ani}^{\mathrm {cub}} = -\frac {2}{\mu _0 M_s} \mathbf {D} \, \mathbf {m}, \label {eq::cubic_field} \end{equation}

with the matrix \(\mathbf {D}\) defined as:

\[ \scalebox {0.9}{$ \mathbf {D} = \begin {bmatrix} K_{c1} (m_2^2 + m_3^2) + K_{c2} m_2^2 m_3^2 & 0 & 0 \\ 0 & K_{c1} (m_1^2 + m_3^2) + K_{c2} m_1^2 m_3^2 & 0 \\ 0 & 0 & K_{c1} (m_1^2 + m_2^2) + K_{c2} m_1^2 m_2^2 \end {bmatrix} $} \]

Anisotropy may originate at the interface in thin films or multilayers due to broken crystal symmetry, lattice strain, or modified electronic structure. This surface anisotropy is localized at the boundary \(\partial \Omega \) and described similarly to the bulk anisotropy energy, but with surface integrals:

\begin{equation} \mathcal {E}_\mathrm {ani}^{\mathrm {surf}} = \int _{\partial \Omega } f(\mathbf {m}, \mathbf {e}) \, \mathrm {d}S, \end{equation}

where \(f\) represents an anisotropy density similar in form to the volume expressions in Equations (3.32), (3.34), and (3.36).

Although the above formulations are phenomenological, they are well established and accurately capture experimental observations  [99]. In many practical cases, the lowest-order terms (e.g., \(K_{u2} = 0\), \(K_{c2} = 0\)) provide a sufficient approximation for modeling anisotropy effects.

3.2.4 Demagnetizing Field

In ferromagnetic materials, the magnetization generates a field contribution known as the demagnetizing field, which arises from long-range dipole-dipole interactions. This field opposes the magnetization, tending to reduce the systems total magnetic moment. The associated energy term is called the demagnetization energy (also referred to as magnetostatic energy or stray-field energy), and it plays a key role in the formation of magnetic domains.

Unlike the exchange or anisotropy fields, the demagnetizing field is non-local. It depends on both the shape of the magnetic body and the spatial distribution and orientation of the magnetization, thereby complicating its evaluation. Nonetheless, it can be derived from Maxwell’s macroscopic equations under the assumption of no free electric currents (\(\mathbf {j}_e = \mathbf {0}\)), which simplifies them to:

\begin{align} \nabla \cdot \mathbf {B} & = 0, \label {eq::maxwell_divB} \\ \nabla \times \mathbf {H}_\mathrm {d} & = \mathbf {0}, \label {eq::maxwell_curlH} \end{align} where \(\mathbf {B}\) is the magnetic flux density and \(\mathbf {H}_\mathrm {d}\) is the demagnetizing field. Using the constitutive relation:

\begin{equation} \mathbf {B} = \mu _0 (\mathbf {H}_\mathrm {d} + \mathbf {M}), \end{equation}

and noting from Equation (3.39b) that \(\mathbf {H}_\mathrm {d}\) is conservative, one can define a scalar magnetostatic potential \(u\) such that:

\begin{equation} \mathbf {H}_\mathrm {d} = -\nabla u. \label {eq::hd_potential} \end{equation}

Substituting into Equation (3.39a) yields:

\begin{equation} \nabla \cdot (-\nabla u + \mathbf {M}) = 0, \end{equation}

which simplifies to a Poisson equation for \(u\):

\begin{equation} \Delta u = \nabla \cdot \mathbf {M}. \label {eq::poisson_u} \end{equation}

Equation (3.43) is solved over the entire space \(\mathbb {R}^3\), with the open boundary condition:

\begin{equation} u(\mathbf {x}) = \mathcal {O}\left (\frac {1}{|\mathbf {x}|}\right ), \quad \text {as } |\mathbf {x}| \to \infty , \label {eq::open_bc} \end{equation}

ensuring that the magnetostatic potential vanishes at infinity [100, 101]. In analogy to electrostatics, one can define effective magnetic charges:

\begin{align} \rho _m & = -\nabla \cdot \mathbf {M} \quad \text {(volume magnetic charge)}, \\ \sigma _m & = \mathbf {M} \cdot \mathbf {n} \quad \text {(surface magnetic charge)}, \end{align} where \(\mathbf {n}\) is the unit normal vector to the boundary \(\partial \Omega \) of the magnetic body. Using these definitions, the magnetostatic potential \(u\) at position \(\mathbf {x}\) can be expressed as:

\begin{equation} u(\mathbf {x}) = \frac {1}{4\pi } \left ( \int _{\Omega } \frac {\rho _m(\mathbf {x}')}{|\mathbf {x} - \mathbf {x}'|} \, \mathrm {d}x' + \int _{\partial \Omega } \frac {\sigma _m(\mathbf {x}')}{|\mathbf {x} - \mathbf {x}'|} \, \mathrm {d}s' \right ). \label {eq::u_green} \end{equation}

From Equation (3.41), the demagnetizing field then follows as:

\begin{equation} \mathbf {H}_\mathrm {d}(\mathbf {x}) = -\nabla u(\mathbf {x}). \label {eq::hd_from_u} \end{equation}

The demagnetization (or magnetostatic) energy represents the interaction energy of the magnetization with its stray field and is given by:

\begin{equation} \mathcal {E}_\mathrm {d} = -\frac {\mu _0 M_s}{2} \int _{\Omega } \mathbf {m} \cdot \mathbf {H}_\mathrm {d} \, \mathrm {d}\mathbf {x}, \label {eq::ed} \end{equation}

where the factor \(1/2\) avoids double-counting the self-interaction.

Alternatively, in numerical micromagnetics (e.g., finite difference methods), it is often convenient to express the demagnetizing field as a convolution:

\begin{equation} \mathbf {H}_\mathrm {d}(\mathbf {x}) = - \int _{\Omega } \mathbf {N}(\mathbf {x} - \mathbf {x}') \cdot \mathbf {m}(\mathbf {x}') \, \mathrm {d}\mathbf {x}', \label {eq::hd_convolution} \end{equation}

where \(\mathbf {N}\) is the demagnetizing tensor kernel that encodes the dipole-dipole interaction  [102]. The components \(N_{ij}\) of this \(3 \times 3\) tensor are given by  [103]:

\begin{equation} N_{ij}(\mathbf {x} - \mathbf {x}') = \frac {M_S}{4\pi } \frac {3\,r_i\, r_j - \delta _{ij}}{|\mathbf {x} - \mathbf {x}'|^3}, \label {eq::demag_tensor} \end{equation}

where \(r_i = (x_i - x'_i)/|\mathbf {x} - \mathbf {x}'|\) and \(r_j = (x_j - x'_j)/|\mathbf {x} - \mathbf {x}'|\) are the directional cosines of the displacement vector, and \(i, j = x, y, z\) denote Cartesian components. This long-range, shape- and geometry-sensitive contribution is crucial in domain formation and stability analysis in micromagnetic systems.

3.2.5 Ampère Field

When a current flows through the magnetic domain, it generates an additional contribution to the effective magnetic field, usually referred to as the Ampère field. This field arises from the interaction between the electric current and the surrounding magnetic material and is governed by the classical Biot-Savart law [97, 104]. For a given current density \(\mathbf {J}_\mathrm {C}\), the Ampère field \(\mathbf {H}_\mathrm {curr}\) at position \(\mathbf {x}\) is computed as:

\begin{equation} \mathbf {H}_\mathrm {curr}(\mathbf {x}) = \frac {1}{4\pi } \int _{\Omega } \frac {\mathbf {J}_\mathrm {C}(\mathbf {x}') \times (\mathbf {x} - \mathbf {x}')}{|\mathbf {x} - \mathbf {x}'|^3} \, \mathrm {d}\mathbf {x}', \label {eq::ampere_field} \end{equation}

where the integral is taken over the current-carrying domain \(\Omega \). This field contributes to the total effective field used in the LLG equation and becomes particularly relevant in spintronic devices such as MRAM.

3.2.6 Thermal Field

At a nonzero temperature, the magnetization of a ferromagnetic material experiences thermal agitation. These thermal fluctuations can be incorporated into the LLG formalism by introducing an auxiliary, randomly varying term in the effective magnetic field, commonly referred to as the thermal field. This additional field accounts for the stochastic nature of thermal excitations and is modeled as a Gaussian random vector field with specified statistical properties.

The thermal field \(\mathbf {H}_{\mathrm {th}}\) must satisfy the following conditions  [105, 106]:

\begin{align} \langle H^i_{\mathrm {th}}(t) \rangle & = 0, \label {eq::thermal_mean} \\ \langle H^i_{\mathrm {th}}(t), H^j_{\mathrm {th}}(t') \rangle & = D\, \delta _{ij} \delta (t - t'), \label {eq::thermal_correlation} \end{align} where \(i\) and \(j\) denote Cartesian components, \(\langle \cdot \rangle \) indicates the ensemble average over different realizations of the stochastic process, \(\delta _{ij}\) is the Kronecker delta, and \(\delta (t - t')\) is the Dirac delta function. These relations describe a field that is uncorrelated in both time and space.

The prefactor \(D\) quantifies the strength of the thermal fluctuations and is derived from the fluctuation-dissipation theorem  [107]:

\begin{equation} D = \frac {2\, \alpha \, k_B T}{\mu _0\, M_s\, \gamma \, V}, \label {eq::thermal_D} \end{equation}

where \(\alpha \) is the Gilbert damping parameter, \(k_B\) is the Boltzmann constant, \(T\) is the absolute temperature, \(M_s\) is the saturation magnetization, \(\gamma \) is the gyromagnetic ratio, and \(V\) is the volume of the computational cell. The thermal field is generated using Gaussian white noise with a zero mean and a variance proportional to \(D\). Its presence leads to a stochastic extension of the LLG equation, often referred to as the stochastic LLG or Langevin LLG model.

3.2.7 Interlayer Exchange Coupling

Magnetic layers within a multilayer stack can experience exchange coupling when separated by a nonmagnetic spacer (NMS), and even across thin TBs. This interaction, initially proposed by Ruderman and Kittel  [108], is mediated by conduction electrons of the nonmagnetic spacer. The theory was later extended by Kasuya  [109] and Yosida  [110], leading to the well-known RKKY interaction, which predicts an oscillatory exchange coupling as a function of spacer thickness, alternating between FM and antiferromagnet (AFM) behavior.

The nature of IEC depends on the electronic properties of the spacer:

• (i) For metallic spacers, the coupling is well-described by the RKKY theory, which explains the indirect interaction between localized spins in the \(d\)- or \(f\)-orbitals via conduction electrons  [111].

• (ii) For semiconducting spacers, exchange interactions are governed by mechanisms such as variable-range hopping or resonant tunneling, mediated by localized states within the band gap  [17, 112].

• (iii) In the case of insulating spacers, strong and non-oscillatory IEC is observed, often attributed to spin-dependent tunneling  [113, 114].

Two theoretical formulations are commonly used to account for these behaviors across different spacer types: discrete spin models such as the Heisenberg exchange interaction, and effective energy expressions such as the one proposed by Bruno  [111]. While the Heisenberg model describes pairwise interactions between localized spins, Bruno’s formulation captures the angular dependence of the coupling energy in a continuum framework. It effectively coarse-grains the microscopic exchange by summarizing spin-dependent electron reflections at magnetic interfaces  [111].

Bruno proposed the following expression for the interlayer coupling energy, based on the angle \(\Delta \phi \) between the magnetization of adjacent magnetic layers (see Figure 3.3):

\begin{equation} \mathcal {E} = -J_1 \cos (\Delta \phi ) - J_2 \cos ^2(\Delta \phi ), \end{equation}

where \(J_1\) and \(J_2\) are coupling constants. The term \(J_1\) represents the bilinear coupling (dominant in most systems), responsible for the oscillatory behavior seen in metallic spacers:

• (i) \(J_1 > 0\): FM coupling, favoring parallel alignment (\(\Delta \phi = 0\))

• (ii) \(J_1 < 0\): AFM coupling, favoring anti-parallel alignment (\(\Delta \phi = \pi \))

The biquadratic term \(J_2\) is often smaller and associated with structural inhomogeneities such as surface roughness or defects, and can induce non-collinear alignment of the magnetizations  [115]. When \(J_2\) is negligible, the expression simplifies to:

\begin{equation} \mathcal {E} = -J_1 \cos (\Delta \phi ), \end{equation}

highlighting the pure bilinear coupling as in classical RKKY-type interactions.

In the micromagnetic continuum approximation, IEC is modeled as an interaction between the two magnetic interfaces \(\partial \Omega _1\) and \(\partial \Omega _2\) on either side of the spacer, which are assumed to be of equal area. Integrating the discrete Heisenberg interaction (3.20) over these surfaces yields the following interlayer exchange energy:

\begin{equation} \mathcal {E}_{\mathrm {iex}} = - \int _{\partial \Omega _1} A\, \mathbf {m}(\mathbf {x}) \cdot \mathbf {m}(P(\mathbf {x}))\, \mathrm {d}s, \label {eq::continuum_iec} \end{equation}

\(A\) is the effective interlayer exchange constant, and \(P:~\partial \Omega _1 \rightarrow \partial \Omega _2\) is a mapping between the corresponding surface points on both interfaces. This continuum representation is consistent with the bilinear term in Bruno’s angular model and arises naturally from the coarse-grained limit of the Heisenberg interaction.

This energy formulation underpins the design of SAFs, where the spacer thickness is tuned to promote AFM coupling (i.e., negative \(A\)). SAFs are widely implemented in spintronic devices such as MTJs and STT-MRAMs, offering enhanced thermal stability and reduced stray fields.

While Equation (3.56) defines the physical energy contribution, its numerical evaluation in a finite element framework poses a challenge: the non-magnetic spacer is often too thin (\(<1\) nm) to be explicitly discretized with volume elements without inducing high computational costs. Therefore, in this work, this interaction is implemented as a specialized boundary condition in the weak form of the effective field. The details of this implementation, including the interface mapping algorithm \(P(\mathbf {x})\), are presented in Section 6.4.

In summary, the microscopic (Heisenberg) and mesoscopic (Bruno) representations of IEC are deeply connected. The discrete Heisenberg exchange interaction captures the fundamental spin-spin coupling mechanism at the atomic level. Bruno’s angular model, derived using quantum perturbation theory applied to spin-dependent electron reflection at magnetic interfaces  [111], effectively approximates this interaction and expresses the total IEC energy as a function of angular misalignment between magnetizations.

This interaction is subsequently formulated within a micromagnetic framework, bridging theory and computation via a surface integral across magnetic interfaces. This continuum form provides a spatially averaged representation of the bilinear IEC and serves as the theoretical foundation for our numerical approach.