RECENT improvements in the development of com-
plimentary metal-oxide-semiconductor (CMOS)
devices have been driven by downscaling of their compo-
nents. This scaling, however, shows signs of saturation,
as it is accompanied by an increase in stand-by power
consumption and leakage currents in conventional CMOS
technology [1]. The introduction of nonvolatile memory
components can help mitigate these issues, as they do
not require refreshing of the memory bits. Spin-transfer
torque (STT) magnetoresistive random access memory
(MRAM) is a viable candidate as a nonvolatile compo-
nent, thanks to its simple structure and compatibility
with CMOS back-end of line processes. It is promising
for several stand-alone, embedded automotive and Inter-
net of Things applications, as well as possible employ-
ments in frame buffer memory and slow SRAM [2]–[5].
The driving factor of the writing process in modern
MRAM cells is the torque generated by the polarization
process of electrons transiting through the ferromagnetic
layers in a magnetic tunnel junction (MTJ). Develop-
ment of reliable simulation tools can provide a valuable
support for the design of advanced MRAM cells. The
simulation of the switching process can be achieved by
solving the Landau-Lifshitz-Gilbert (LLG) equation for
magnetization dynamics, with the inclusion of a term
describing the torque acting on the magnetization. Such
term is computed from the non-equilibrium spin accu-
mulation Sin the structure. A solution for Sin all
non-magnetic and ferromagnetic layers of an MRAM cell
can be obtained by means of the spin and charge drift-
diffusion formalism [6]–[8]. Analytical solutions to the
drift-diffusion equations, coupled to the LLG, are only
possible in simplified scenarios, and numerical methods
are necessary to resolve the dynamics in a more general
sense. As the finite element (FE) method is naturally
able to handle meshes with complex geometries and sev-
eral material domains [9],[10], in this work we show how
it can be employed for the implementation of a solver
capable of handling charge, spin accumulation and mag-
netization dynamics. The implementation was carried
out by employing the open-source C++ finite-element li-
brary MFEM [11], and can be applied to several MRAM
structures.
The LLG equation for the description of the magne-
tization dynamics is
∂m
∂t =−|γ|µ0m×Heff +αm×∂m
∂t +1
MS
TS.(1)
m=M/MSis the unit vector in the direction of the lo-
cal magnetization, MSis the saturation magnetization,
γis the gyromagnetic ratio, and µ0is the vacuum per-
meability. Heff is an effective field including several con-
tributions, like the external field, the exchange coupling,
the demagnetizing field and the anisotropy field. TSis
the spin torque, which can be obtained from the spin ac-
cumulation and computed through the spin-charge drift
diffusion equations [6], [8]
−∇ · ˜
JS−De
S
λ2
sf
−TS= 0 ,(2a)
˜
JS=−µB
eβσm⊗JC−e
µB
βDDe(∇S)Tm−De∇S,
(2b)
TS=−De
m×S
λ2
J
−De
m×(m×S)
λ2
φ
.(2c)
˜
JSis the spin current tensor, Deis the electron diffusion
coefficient, λsf is the spin-diffusion length, µBis the Bohr
Finite Element Method for MRAM Switching Simulations
1,2S. FIORENTINI, 1R. L. DE ORIO, 1,2J. ENDER, 1S. SELBERHERR,
1,2M. BENDRA, 1,2N. JØRSTAD, 3WOLFGANG GOES, 1,2V. SVERDLOV
1Christian Doppler Laboratory for Nonvolatile Magnetoresisitive Memory and Logic at the
2Institute for Microelectronics, TU Wien, Gußhausstraße 27–29/E360, 1040 Vienna, AUSTRIA
3Silvaco Europe Ltd., Cambridge, UNITED KINGDOM
Abstract: - The development of reliable simulation tools provides a valuable help in the design of modern
MRAMdevices. Thanks to its versatility in the choice of meshes and discretization, the finite element method
is a useful framework for the numerical solution of the magnetization dynamics. We review a finite element
implementation of both the Landau-Lifshitz-Gilbert equation and the spin and charge drift-diffusion formalism
in a solver employing open source software. The presented approach is successfully applied to emerging multi-
layered MRAM cells.
Keywords: - Finite element method, Landau- Lifshitz-Gilbert equation, spin and charge driftdiffusion, STT-
MRAM
Received: March 22, 2022. Revised: November 15, 2022. Accepted: December 13, 2022. Published: December 31, 2022.
1. Introduction
2. Micromagnetic Equations
WSEAS TRANSACTIONS on SYSTEMS and CONTROL
DOI: 10.37394/23203.2022.17.64
S. Fiorentini, R. L. De Orio, J. Ender,
S. Selberherr, M. Bendra, N. Jørstad,
Wolfgang Goes, V. Sverdlov
E-ISSN: 2224-2856
585
Volume 17, 2022
magneton, eis the elementary charge, βσand βDare
polarization parameters, JCis the charge current density,
λJis the exchange length and λφis the spin dephasing
length.
First schemes for the implementation of a FE algo-
rithm were proposed in [12],[13], by considering only the
exchange field contribution to the effective field. In refer-
ence [14] the so-called tangent plane integrator scheme,
solving for the discrete time derivative, was introduced,
and was later generalized to include the full effective field
[15], [16]. The unconditional convergence of the tangent
plane integrator scheme and of the finite-element imple-
mentation of the spin and charge drift-diffusion equa-
tions was proven in [17], and later applied to metallic
spin-valves in [6]. The weak formulation of the tangent
plane integrator and of the spin and charge transport
equations are hereby reported.
The tangent plane scheme solves the LLG equation
for the magnetization derivative ∂m/∂t =v. The weak
formulation of the standard LLG equation, without the
inclusion of the torque term TS, comes from the expres-
sion
α∂m
∂t +m×∂m
∂t =|γ|µ0Heff + (m·Heff)m,(3)
which can be obtained by cross-multiplying
(1) with m, and using the product rule
a×(b×c)=(c·a)b−(a·b)ctogether with the
constraint |m|= 1. For the finite element implemen-
tation, the magnetization is taken to be a piecewise
affine, globally continuous function [18]. Each com-
ponent belongs to the Sobolev space H1. It is the
space of functions in L2that additionally admit a
weak gradient which also belongs to L2[16]. The
notation for the vector space of the magnetization is
H1. Instead of looking for the solution vin the same
space of the magnetization, the solution space VTis
restricted to vectors tangent to the magnetization, so
that VT=w∈H1|m·w= 0. The test functions
are also restricted to the same space, so that the weak
formulation of (3) results in
Zω
(αv+m×v)·wdx=|γ|µ0Zω
Heff(m)·wdx.
(4)
ωis the subdomain containing only the magnetic parts of
the structure under analysis. The last term on the right-
hand side of (3) is not present, as the test functions are
restricted to the tangent space VT. The time derivative
vat a certain time tkis obtained by setting [18]
mk+1 =mk+θδtv.(5)
δt is the time-step, and θis a parameter between 0 and
1, with the value 0 leading to a fully explicit scheme and
1 to a fully implicit one. Each effective field contribution
can be treated with a different value of θ. In the imple-
mentation employed for this work, θis different from 0
only for the exchange field contribution, were the value 1
is employed for stability reasons [19]. The weak formula-
tion employed for the computation of the magnetization
dynamics by the FE solver is then, with the inclusion of
the torque terms coming from (2):
Zωαv+mk×v·wdx+2A|γ|
MS
δt Zω
∇v:∇wdx=
−2A|γ|
MSZω
∇mk:∇wdx+
γ0Zω
Heff ·wdx+De
MSZωSk
λ2
J
+mk×Sk
λ2
φ·wdx
(6a)
mk+1 =mk+δt v
|mk+δt v|(6b)
The right-hand side of (6b) is evaluated nodewise, and
∇a:∇b=Pij (∂ai/∂xj)(∂bi/∂xj) is the Frobenius in-
ner product of two matrices. The exchange field contri-
bution, together with (5), gives rise to the second term
on the left-hand side and to the first term on the right-
hand side. The boundary integrals arising from the weak
formulation of this contribution are put to zero by the
Neumann boundary condition (∇m)n=0(with nthe
boundary normal), applied on the external boundary of
the magnetic region ∂ω.Heff contains the remaining ef-
fective field contributions. The equations are subject
to the initial condition m(0) = m0. The system of
equations resulting from the FE implementation of this
weak formulation includes the tangent plane constraint
m·w= 0, and the solution at each time-step is com-
puted through a solver based on the generalized mini-
mal residual (GMRES) method, provided by the library
MFEM. The GMRES method is designed for indefinite
non-symmetric systems of linear equations, as is the case
of FE implementation of (6a) due to the presence of the
cross-product terms. Material parameters that can differ
from subdomain to subdomain are treated as piecewise
constant functions, unless differently stated. The con-
tribution of the demagnetizing field is evaluated only on
the disconnected magnetic domain by using a hybrid ap-
proach combining the boundary element method and the
FE method [20], [21].
The weak formulation for the computation of the spin
accumulation is derived from (2). The equations are
solved for the magnetization mkat each time-step, and
the resulting spin accumulation Skis employed in (6a).
3. Weak Formulation for a Finite
Element Implementation
3.1 LLG Equation
3.2 Spin and Charge Drift-Diffusion
Equations
WSEAS TRANSACTIONS on SYSTEMS and CONTROL
DOI: 10.37394/23203.2022.17.64
S. Fiorentini, R. L. De Orio, J. Ender,
S. Selberherr, M. Bendra, N. Jørstad,
Wolfgang Goes, V. Sverdlov
E-ISSN: 2224-2856
586
For the computation of the charge potential and cur-
rent, a Laplace equation is solved for the whole domain
Ω, with the weak formulation
ZΩ
σ∇V· ∇vdx= 0 ,(7a)
ZΩ
JC·vdx=−ZΩ
σ∇V·vdx.(7b)
vrepresents a test function belonging to H1, while vis
a test function in H1.Vbelongs to H1, and equation
(7b) is employed to obtain a projection of JCin the H1
function space [18]. Dirichlet conditions are applied to
prescribe the voltage at the contacts. The Neumann
condition σ∇V·n= 0 is assumed on external bound-
aries not containing an electrode. Special treatment
of the conductivity σresults in an implementation
capable of reproducing the charge current dependence
on the tunnel magnetoresistance (TMR) and relative
direction of the magnetization vectors in the free layer
(FL) and reference layer (RL) of an MTJ [22]. The
solution of the system of equations resulting from the
FE implementation of (7) is computed through a solver
based on the conjugate gradient (CG) method, provided
by the library MFEM. The CG method is designed for
the numerical solution of systems of linear equations
whose matrix is positive-definite, as is the case of the
FE implementation of (7).
The weak formulation of the spin drift-diffusion equa-
tions presented in [6] was generalized to apply to (2)
which includes the additional spin dephasing term, re-
sulting in the following expression:
DeZΩ∇S−βσβDm⊗(∇S)Tm:∇vdx+
+DeZΩ S
λ2
sf
+S×m
λ2
J
+m×(S×m)
λ2
φ!·vdx=
µB
eβσZω
(m⊗JC) : ∇vdx+
−µB
eβσZ∂Ω∩∂ω
((m⊗JC)n)·vdx.(8)
vrepresents again a test function belonging to H1,S
also belongs to H1, and ∂Ω∩∂ω indicates the shared
external boundary of the whole domain Ω and magnetic
subdomain ω. The boundary integrals arising from
the weak formulation are put to zero by the Neumann
condition (∇S)n=0, assumed on external boundaries.
For contacting regions longer than the spin-flip length,
this condition is equivalent to an exponential decay of
Stowards the electrodes [6], [8]. The charge current
JCis the one computed from (7). The addition of
appropriate boundary conditions at the interfaces
between the tunnel barrier and the RL and FL allows
to reproduce the torque properties typical of MTJs [23].
The solution of the system of equations resulting from
the FE implementation of (8) is computed through a
solver based on the GMRES method, provided by the
Fig. 1: Anti-parallel to parallel switching of elongated
MRAM cell with composite FL. The two FL segments
switch one at a time, creating visible steps in the switch-
ing process.
library MFEM, as due to the cross product terms the
FE implementation of (8) results in a non-symmetric
system of linear equations.
The solver implementing the proposed approach can
be applied to various structures composed of ferromag-
nets, metal spacers, and tunnel barriers. It is particularly
suitable for simulating recently proposed devices pre-
senting additional layers for the reduction of switching
current and cell size. Application of the solver to elon-
gated structures with composite FL, proposed in [24],
reveals a step-by-step switching process [23], and allows
to evaluate the possibility of using such structures as
multi-bit cells. The anti-parallel to parallel switching of
a structure with two FL segments is reported in Fig.1.
The solver was also successfully applied to the switch-
ing simulation of structure presenting a double FL [25],
reproducing the observed reduction of the switching volt-
age [26]. Additional terms in the spin current expression
allow to also reproduce the spin Hall effect [8], so that
the solver can also be applied to switching simulations
of spin-orbit torque (SOT) MRAM [27].
We presented a finite element implementation of the
weak formulation of the LLG equation coupled to spin
transport in a solver based on open source software.
The developed solver has been successfully applied to
switching simulations of emerging experimental struc-
tures composed of multiple ferromagnetic layers, non-
magnetic spacers and tunnel barriers, and can be a valu-
able help to predict and investigate the switching behav-
ior of novel structures.
The financial support by the Austrian Federal Min-
istry for Digital and Economic Affairs and the National
Foundation for Research, Technology and Development
is gratefully acknowledged.
4. Applications
5. Conclusion
Acknowledgment
WSEAS TRANSACTIONS on SYSTEMS and CONTROL
DOI: 10.37394/23203.2022.17.64
S. Fiorentini, R. L. De Orio, J. Ender,
S. Selberherr, M. Bendra, N. Jørstad,
Wolfgang Goes, V. Sverdlov
E-ISSN: 2224-2856
587
[1] T. Hanyu, T. Endoh, D. Suzuki, H. Koike, Y. Ma
et al., “Standby-power-free integrated circuits using
MTJ-based VLSI computing,” Proc. IEEE, vol. 104,
no. 10, pp. 1844–1863, 2016.
[2] W. J. Gallagher, E. Chien, T. Chiang, J. Huang,
M. Shih et al., “22nm STT-MRAM for reflow
and automotive uses with high yield, reliability,
and magnetic immunity and with performance and
shielding options,” in Proc. IEDM Conf., 2019, pp.
2.7.1–2.7.4.
[3] S. H. Han, J. M. Lee, H. M. Shin, J. H. Lee,
K. S. Suh et al., “28nm 0.08 mm2/Mb embedded
MRAM for frame buffer memory,” in Proc. IEDM
Conf., 2020, pp. 11.2.1–11.2.4.
[4] Y.-C. Shih, C.-F. Lee, Y.-A. Chang, P.-H. Lee,
H.-J. Lin et al., “A reflow-capable, embedded 8Mb
STT-MRAM macro with 9ns read access time in
16nm FinFET logic CMOS process,” in Proc. IEDM
Conf., 2020, pp. 11.4.1–11.4.4.
[5] V. B. Naik, K. Yamane, T. Lee, J. Kwon,
R. Chao et al., “JEDEC-qualified highly reliable
22nm FD-SOI embedded MRAM for low-power
industrial-grade, and extended performance to-
wards automotive-grade-1 applications,” in Proc.
IEDM Conf., 2020, pp. 11.3.1–11.3.4.
[6] C. Abert, M. Ruggeri, F. Bruckner, C. Vogler,
G. Hrkac et al., “A three-dimensional spin-diffusion
model for micromagnetics,” Sci. Rep., vol. 5, no. 1,
p. 14855, 2015.
[7] C. Abert, M. Ruggeri, F. Bruckner, C. Vogler,
A. Manchon et al., “A self-consistent spin-diffusion
model for micromagnetics,” Sci. Rep., vol. 6, no. 1,
p. 16, Dec. 2016.
[8] S. Lepadatu, “Unified treatment of spin torques us-
ing a coupled magnetisation dynamics and three-
dimensional spin current solver,” Sci. Rep., vol. 7,
no. 1, p. 12937, 2017.
[9] D. Braess, Finite Elements: Theory, Fast Solvers,
and Applications in Solid Mechanics, 3rd ed. Cam-
bridge University Press, 2007.
[10] M. G. Larson and F. Bengzon, The Finite Ele-
ment Method: Theory, Implementation, and Ap-
plications. Springer Berlin Heidelberg, 2013.
[11] R. Anderson, J. Andrej, A. Barker, J. Bramwell,
J.-S. Camier et al., “MFEM: A modular finite ele-
ment library,” Comp. & Math. with Appl., 2020.
[12] F. Alouges and P. Jaisson, “Convergence of a finite
element discretization for the Landau-Lifshitz equa-
tions in micromagnetism,” Math. Models Methods
Appl. Sci., vol. 16, no. 2, p. 299 – 316, 2006.
[13] S. Bartels, J. Ko, and A. Prohl, “Numerical
analysis of an explicit approximation scheme
for the Landau-Lifshitz-Gilbert equation,” Math.
Comput., vol. 77, no. 262, p. 773 – 788, 2008.
[14] F. Alouges, “A new finite element scheme for
Landau-Lifchitz equations,” Discrete Contin. Dyn.
Syst. S, vol. 1, no. 2, pp. 187–196, 2008.
[15] F. Alouges, E. Kritsikis, and J.-C. Toussaint, “A
convergent finite element approximation for Lan-
dau–Lifschitz–Gilbert equation,” Physica B, vol.
407, no. 9, pp. 1345–1349, 2012.
[16] F. Bruckner, D. Suess, M. Feischl, T. F¨uhrer,
P. Goldenits et al., “Multiscale modeling in micro-
magnetics: Existence of solutions and numerical
integration,” Math. Models Methods Appl. Sci.,
vol. 24, no. 13, pp. 2627–2662, 2014.
[17] C. Abert, G. Hrkac, M. Page, D. Praetorius,
M. Ruggeri, and D. Suess, “Spin-polarized transport
in ferromagnetic multilayers: An unconditionally
convergent FEM integrator,” Comp. & Math. with
Appl., vol. 68, no. 6, pp. 639 – 654, 2014.
[18] C. Abert, “Micromagnetics and spintronics: Models
and numerical methods,” Eur. Phys. J. B, vol. 92,
no. 6, p. 120, June 2019.
[19] G. Hrkac, C.-M. Pfeiler, D. Praetorius, M. Ruggeri,
A. Segatti, and B. Stiftner, “Convergent tangent
plane integrators for the simulation of chiral
magnetic skyrmion dynamics,” Adv. Comput.
Math., vol. 45, no. 3, pp. 1329–1368, June 2019.
[20] J. Ender, M. Mohamedou, S. Fiorentini, R. L.
de Orio, S. Selberherr et al., “Efficient demagne-
tizing field calculation for disconnected complex ge-
ometries in STT-MRAM cells,” in Proc. SISPAD
Conf., 2020, pp. 213–216.
[21] M. Bendra, J. Ender, S. Fiorentini, T. Hadamek,
R. L. de Orio et al., “Finite element method ap-
proach to MRAM modeling,” in Proc. MIPRO
Conf., 2021, pp. 70–73.
[22] S. Fiorentini, J. Ender, S. Selberherr, R. L. de
Orio, W. Goes, and V. Sverdlov, “Coupled spin and
charge drift-diffusion approach applied to magnetic
tunnel junctions,” Solid-State Electron., vol. 186,
p. 108103, 2021.
[23] S. Fiorentini, M. Bendra, J. Ender,
R. L. de Orio, W. Goes et al., “Spin and
charge drift-diffusion in ultra-scaled MRAM cells,”
Sci. Rep., vol. 12, no. 1, p. 20958, Dec. 2022.
[24] B. Jinnai, J. Igarashi, K. Watanabe, T. Funatsu,
H. Sato et al., “High-performance shape-anisotropy
magnetic tunnel junctions down to 2.3 nm,” in Proc.
IEDM Conf., 2020, pp. 24.6.1–24.6.4.
[25] G. Hu, G. Lauer, J. Z. Sun, P. Hashemi, C. Safranski
et al., “2X reduction of STT-MRAM switching cur-
rent using double spin-torque magnetic tunnel junc-
tion,” in Proc. IEDM Conf., 2021, pp. 2.5.1–2.5.4.
[26] W. J. Loch, S. Fiorentini, N. P. Jørstad, W. Goes,
S. Selberherr, and V. Sverdlov, “Double refer-
ence layer STT-MRAM structures with improved
performance,” Solid-State Electron., vol. 194,
p. 108335, 2022.
[27] N. P. Jørstad, S. Fiorentini, W. J. Loch,
W. Goes, S. Selberherr, and V. Sverdlov, “Finite
element modeling of spin–orbit torques,” Solid-
State Electron., vol. 194, p. 108323, 2022.
References
Creative Commons Attribution License 4.0
(Attribution 4.0 International, CC BY 4.0)
This article is published under the terms of the Creative
Commons Attribution License 4.0
https://creativecommons.org/licenses/by/4.0/deed.en_US
WSEAS TRANSACTIONS on SYSTEMS and CONTROL
DOI: 10.37394/23203.2022.17.64
S. Fiorentini, R. L. De Orio, J. Ender,
S. Selberherr, M. Bendra, N. Jørstad,
Wolfgang Goes, V. Sverdlov
E-ISSN: 2224-2856
588