# Energy Barriers for Thermally Activated Magnetization Reversal in
Perpendicularly Magnetized Nanodisks in a Transverse Field

Corrado C. M. Capriata<EMAIL_ADDRESS>Division of Electronics and Embedded
Systems, KTH - Royal Institute of Technology, Stockholm, Sweden  Center for
Quantum Phenomena, Department of Physics, New York University, New York, NY
10003, USA  B. Gunnar Malm  Division of Electronics and Embedded Systems, KTH
- Royal Institute of Technology, Stockholm, Sweden  Andrew D. Kent  Center for
Quantum Phenomena, Department of Physics, New York University, New York, NY
10003, USA  Gabriel D. Chaves-O’Flynn<EMAIL_ADDRESS>Institute of
Molecular Physics, Polish Academy of Sciences, Poznań, Poland

###### Abstract

Thermally-induced transitions between bistable magnetic states of magnetic
tunnel junctions (MTJ) are of interest for generating random bitstreams and
for applications in stochastic computing. An applied field transverse to the
easy axis of a perpendicularly magnetized MTJ (pMTJ) can lower the energy
barrier ($E_{b}$) to these transitions leading to faster fluctuations. In this
study, we present analytical and numerical calculations of $E_{b}$ considering
both coherent (macrospin) reversal and non-uniform wall-mediated magnetization
reversal for a selection of nanodisk diameters and applied fields. Non-uniform
reversal processes dominate for larger diameters, and our numerical
calculations of $E_{b}$ using the String method show that the transition state
has a sigmoidal magnetization profile. The latter can be described with an
analytical expression that depends on only one spatial dimension, parallel to
the applied field, which is also the preferred direction of profile motion
during reversal. Our results provide nanodisk energy barriers as a function of
the transverse field, nanodisk diameter, and material characteristics, which
are useful for designing stochastic bitstreams.

##  I Introduction

There has been a recent surge in interest in innovative computational
approaches that mimic the flexibility of neural systems [1, 2, 3]. Proposals
based on magnetic systems include neuromorphic computing [4], reservoir [5],
and stochastic computing [6]. For these applications, it is necessary to find
physical systems with true random behavior [7, 8, 9]. Magnetic tunnel junction
devices can exhibit random two-state fluctuations. However, magnetic tunnel
junctions have primarily been used for traditional information storage, where
the state of the information bit must remain stable over time and be resistant
to thermal fluctuations for decades [10, 11]. To achieve this goal,
perpendicular magnetic tunnel junctions (pMTJs) have been extensively studied
and optimized as memory elements [12, 13]. In contrast, devices suitable for
random number generation require energy barriers that can be easily overcome
by thermal perturbations. Fortunately, pMTJs can also be designed to generate
random numbers with a high rate of switching between two equally likely
metastable states [14, 15, 16, 17, 18, 9].

Kramers’ theory for thermally activated escape over a potential barrier [19,
20] is useful in this regard. This theory states that the rate of thermally
induced switching between two metastable states, $\Gamma$, obeys an Arrhenius
law $\Gamma=\Gamma_{0}\exp(-E_{b}/(k_{B}T))$, where $\Gamma_{0}$ is an attempt
frequency, typically in the GHz range, $k_{B}$ is Boltmann’s constant and $T$
is the temperature. $E_{b}$ is the energy barrier, the difference between the
energy of the transition state and the energy of the metastable state, where
the transition state corresponds to the lowest energy saddle point between the
two metastable states.

Previous research on randomly fluctuating magnetic tunneling junctions (MTJs)
has primarily focused on easy-plane systems, which have been considered the
faster alternative [21, 22]. However, in this article, we present a
theoretical model for MTJs with uniaxial magnetic anisotropy in the presence
of a transverse field. In fact, in pMTJs, applying a transverse magnetic field
lowers the energy barrier, which in turn increases the fluctuation rates [23].

We investigated perpendicularly magnetized nanodisks, that correspond to the
so-called free layer of a perpendicular magnetic tunnel junction (pMTJ) memory
element. Our starting point is a macrospin model developed by Garanin in Ref.
[24]. However, our analysis accounts for non-uniform magnetization switching,
which is characteristic for nanodisks larger than a certain minimum size. The
theory presented in this study predicts how $E_{b}$ varies with nanodisk
diameter and the magnitude of the applied field. To validate our analytical
predictions, we compared them against numerical results obtained using the
String method [25, 26] implemented in OOMMF micromagnetic simulator [27]. This
technique has been previously used by our group to study thin films [28],
nanomagnets [29], and nanorings [30, 31].

The paper is organized as follows. In Section II, we describe the structure we
model and simulate. In Section III, we present the macrospin model and the
micromagnetic simulation results used to test it. Section IV covers the non-
uniform reversal model and the String method technique used to assess the
model. The results of the String method simulation are presented and discussed
in detail in Section V. The paper ends with a summary in Section VI.

##  II Free layer geometry

A schematic of a pMTJ free layer is shown in Fig. 1. It is modeled by a
ferromagnetic nanodisk with geometrical and magnetic properties chosen to have
net magnetic anisotropy perpendicular to the film plane. The coordinate system
has $\bm{\hat{x}}$ parallel to the applied magnetic field and $\bm{\hat{z}}$
perpendicular to the film plane. The magnetocrystalline anisotropy axis
$\mathbf{\hat{k}}$ is parallel to $\mathbf{\hat{z}}$. To describe the
magnetization, we use spherical coordinates in which $\Theta$ is the angle of
the magnetization from the $\bm{\hat{z}}$ axis. Because of the orientation of
the field, the magnetization lies in the $x-z$ plane.

We use SI units and dimensionless expressions are used to simplify the algebra
and provide generalizations. Dimensionless quantities are expressed in lower
case if there is no ambiguity ($\mathbf{m}=\mathbf{M}/M_{s}$), or with tildes
if necessary (as in $\tilde{x}=x/l_{\mathrm{ex}}$).

Figure 1: Schematic of the ferromagnetic disk under consideration. It has
radius $R$, thickness $t$ and an easy magnetic axis in the z-direction (i.e.,
perpendicular magnetic anisotropy). The field is applied in-plane along the
x-direction.

##  III Macrospin model

For magnetization reversal in very small samples coherent rotation of the
magnetization (macrospin reversal) can be favored, while in samples larger
than the exchange length there may be spatial variations of the magnetization
(non-uniform reversal). Thus, we start here with the macrospin model applied
to uniformly magnetized thin ferromagnetic disks with radius $R$, thickness
$t$, and volume $V$.

The magnetic energy density, $\mathcal{E}=E/V$, is obtained from the sum of
dipolar ($\mathcal{E}_{\mathrm{d}}$), anisotropy ($\mathcal{E}_{\mathrm{k}}$),
and Zeeman ($\mathcal{E}_{\mathrm{Z}}$) energy densities:

$\displaystyle\mathcal{E}(\mathbf{m})$
$\displaystyle=\mathcal{E}_{\mathrm{d}}+\mathcal{E}_{\mathrm{k}}+\mathcal{E}_{\mathrm{Z}}.$
(1)

The individual terms depend on $\mathbf{m}$ as follows:

$\displaystyle\mathcal{E}_{\mathrm{d}}$ $\displaystyle=\zeta
K_{d}(\mathbf{m\cdot\hat{z}})^{2}$ (2) $\displaystyle\mathcal{E}_{\mathrm{k}}$
$\displaystyle=-K(\mathbf{m\cdot\hat{k}})^{2}$ (3)
$\displaystyle\mathcal{E}_{\mathrm{Z}}$
$\displaystyle=-\mu_{0}M_{s}\mathbf{H}\cdot\mathbf{m},$ (4)

where $M_{s}$ is the saturation magnetization, $K$ is the magnetic anisotropy,
$K_{d}=\mu_{0}M_{s}^{2}/2$ is the demagnetization ﬁeld energy coeﬃcient, and
$\mathbf{H}=H_{x}\mathbf{\hat{x}}$ is the applied in-plane magnetic field.

The coefficient $\zeta$ is obtained from the disk magnetostatic tensor using
$\zeta=(3N_{zz}-1)/2$. Both are purely geometric factors that depend only on
the disk aspect ratio and are close to unity if $(t/R)\ll 1$ (see Refs. [32,
28] and Appendix A). The applied field will be normalized using the effective
perpendicular anisotropy of the layer
$\mathbf{h}\equiv\mathbf{H}/[M_{s}(Q-\zeta)]$ with $Q=K/K_{d}$.

After expressing the magnetization unit vector in polar coordinates
$\mathbf{m}=\left(\cos\Phi\sin\Theta,\sin\Phi\sin\Theta,\cos\Theta\right)$ and
energy densities in terms of $\mathcal{E}_{0}=K_{d}(Q-\zeta)$ we obtain a
rescaled energy density, $\tilde{\mathcal{E}}=\mathcal{E}/\mathcal{E}_{0}$:

$\tilde{\mathcal{E}}=-\cos^{2}\Theta-2h\sin\Theta\cos\Phi.$ (5)

Because the calculation of energy barriers is the main goal of this work, we
emphasize our conventions at this point. The total energy is obtained after
integrating over the volume of the disk and will be denoted by the unscripted
letter $E$. The problem characteristic energy, $E_{0}$, will be introduced in
Section IV. Energy barriers, $E_{b}$, will be expressed in units of $k_{B}T$
(with $T=$300\text{\,}\mathrm{K}$$). Energy density terms are denoted with
script character $\mathcal{E}$ and are of the order of $1\ \mathrm{MJ/m^{3}}$.
Dimensionless energy densities are denoted with tilde $\tilde{\mathcal{E}}$
after normalization by $\mathcal{E}_{0}$. A summary of the conversion
relations is provided in Table 2. We now proceed to find the extremals of the
energy.

Table 1: Summary of conversion relations for dimensionless units. Symbol |  Variable |  Conversion Relation | Normalization Quantity
---|---|---|---
$E$ |  Energy |  $E=E_{0}\tilde{E}$ | $E_{0}=At$
$\mathcal{E},\mathcal{E}_{d},\mathcal{E}_{k},\mathcal{E}_{Z},\mathcal{E}_{\mathrm{ex}}$ |  Energy density terms, in order: total, dipolar, anisotropy, Zeeman, exchange. |  $\mathcal{E}=\text{$\mathcal{E}_{0}$$\tilde{\mathcal{E}}$}$ | $\mathcal{E}_{0}=K_{d}(Q-\zeta)$
$\mathbf{H}=H\mathbf{x}$ |  Applied external field. |  $H=hM_{s}(Q-\zeta)$ | $M_{s}(Q-\zeta)$
$\tilde{x},\tilde{y}$ |  In plane dimensionless coordinates. |  $x=l_{\mathrm{ex}}\tilde{x}$ | $l_{ex}=\sqrt{\frac{2A}{\mu_{0}M_{s}^{2}\left(Q-\zeta\right)}}$

The transition and metastable states of the energy are magnetization
directions where the effective field,
$\frac{\delta\mathcal{E}}{\delta\mathbf{m}}$, vanishes:

$\displaystyle\frac{\partial\tilde{\mathcal{E}}}{\partial\Theta}=-\cos\Theta\left(\sin\Theta+h\cos\Phi\right)=0$
(6)
$\displaystyle\frac{1}{\sin\Theta}\frac{\partial\tilde{\mathcal{E}}}{\partial\Phi}=-h\sin\Phi=0$
(7)

From Eq. 7 we conclude that $\Phi=0$. There are four uniform solutions of Eq.
6. The first two

$\Theta_{\pm}=\pm\frac{\pi}{2};\qquad\tilde{\mathcal{E}}_{\pm}=\mp 2h,$ (8)

are the polar angles of the energy maximum and saddle. There are also two
mutually supplementary angles that share the same energy

$\Theta_{\uparrow,\downarrow}=\sin^{-1}h\qquad\tilde{\mathcal{E}}_{1}=-\left(1+h^{2}\right),$
(9)

and correspond to energy minima. In Garanin’s macrospin reversal model [24],
the energy barrier for the most likely thermally activated transition is given
by:

$\tilde{\mathcal{E}_{b}}=\tilde{\mathcal{E}}_{+}-\tilde{\mathcal{E}}_{1}=\left(1-h\right)^{2}=\frac{E_{b}}{K_{d}V(Q-\zeta)}$
(10)

For the reader’s convenience, we write $E_{b}$ explicitly in dimension-full
units.

$E_{b}=\left(1-\frac{H}{M_{s}(Q-\zeta)}\right)^{2}K_{d}V(Q-\zeta)$ (11)

Our model slightly refines Garanin’s prediction by explicitly considering
shape anisotropy effects caused by the magnetostatic interaction which are
captured here in the size-dependent factor $(Q-\zeta)$.

###  III.1 Micromagnetic Calculations

To exemplify the predictions of the macrospin we perform pairs of overdamped
micromagnetic simulations (damping factor $\alpha=1$) using
$\mathbf{m}\approx\mathbf{\hat{x}}$ as the initial configuration. For these
simulations, the initial magnetization was set marginally out of the $x-y$
plane and in opposite directions so that the magnetization relaxes to
different metastable states. The parameters for these simulations were set to
match those of our previous experimental studies [33]. They are as follows:
Heisenberg exchange constant
$A=$$4.2\text{\,}\mathrm{pJ}\text{\,}{\mathrm{m}}^{-1}$, saturation
magnetization $M_{s}=$$0.58\text{\,}\mathrm{MA}\text{\,}{\mathrm{m}}^{-1}$,
bulk magnetic anisotropy
$K=$$0.39\text{\,}\mathrm{M}\mathrm{J}\mathrm{/}\mathrm{m}^{3}$. The material
constant becomes $Q=1.84$. The cells in the simulation had dimensions:
$2.5\times 2.5\times 2.6~{}\mathrm{nm^{3}}$.

Figure 2:  (a) Two overdamped micromagnetic simulations of a
$60\text{\,}\mathrm{nm}$ device at a field of $0.6\text{\,}\mathrm{T}$ that
relax to different minima. The energy barrier is measured directly from the
difference between the initial and final states (b). The slow initial
evolution confirms that the macrospin saddle is a critical point of the energy
landscape. The two minima are symmetric and follow identical curves. (c) The
non-uniformity as defined in the text remain low and overlapping for the two
simulations.

A typical result from the overdamped simulations is shown in Fig. 2. The
difference between the two traces can only be seen in the plot of $\Theta(t)$,
Fig. 2(a), because the energy and non-uniformity results are identical. This
behavior can be taken as confirmation that the initial magnetization is indeed
a saddle point of the energy landscape. The initial magnetization dynamics are
at first very slow before decaying to either of the energy minima, as is
evident in the plateau of the energy vs time curve of Fig. 2(b).

We use the spatial average of each magnetization component to quantify the
non-uniformity of the configuration,
$\epsilon=\sqrt{1-\bar{m_{x}}^{2}-\bar{m_{y}}^{2}-\bar{m_{z}}^{2}}$. As can be
seen in Fig. 2(c), this quantity remains small during the overdamped
relaxation. This fact supports our use of overdamped dynamics to validate
macrospin predictions.

From a systematic repetition of this procedure, we determine the energy
barrier versus applied transverse field for six distinct diameters in the
range $5<D<100$ $\text{\,}\mathrm{nm}$. The energy barrier and field are
scaled as indicated in Eq. 11. In consequence, all data appear on a single
curve Fig. 3(b). This result confirms the universality of Garanin’s macrospin
approximation, provided the proper shape-induced magnetization corrections are
taken into account.

For comparison, Fig. 3(a) shows the same data in SI units. An immediate
consequence of these studies is the counter-intuitive result observable here.
It is generally expected that the energy barrier of a uniformly magnetized
system increases with volume. However, it is shown that the energy barrier for
fixed disk thickness actually decreases with the radius for any fixed field.
This is because the effective perpendicular anisotropy weakens as
$R/t\rightarrow\infty$. In Fig. 3(a), the simulation data points do not
exactly follow the theoretical lines. We did not find a clear relation between
the offset between the macrospin model and the theory and the radius of the
device. Therefore we attribute the mismatch to effects of the numerical
discretization.

The decrease in the energy barrier occurs even under the assumption of uniform
magnetization and the curves seem to reach convergence at $D\approx$
$100\text{\,}\mathrm{nm}$. Nonetheless, this limiting energy barrier is much
larger than a barrier obtained from domain wall mediated reversal, as will be
introduced in the next section.

Figure 3:  (a) Energy barrier $\mathcal{E}_{b}$ vs field $\mu_{0}H_{z}$
without the shape factor correction $(Q-\zeta)$. The black curves show the
limiting values for $D=0$ and $D\rightarrow\infty$. Points are from
micromagnetic simulation, continuous lines are obtained from Eq. 11. (b) Non-
dimensional energy barriers, $\tilde{\mathcal{E}_{b}}$, vs scaled field, $h$.
All results from overdamped simulations collapse onto the macrospin model
after appropriate normalization.

##  IV Domain Wall Mediated Reversal

While the macrospin model is a useful approximation for device sizes smaller
than the exchange length, it is natural to expect that it fails once the
system is large enough to accommodate variations of the magnetization. To
consider this scenario, we rewrite the energy as a surface integral with four
spatially varying contributions:

$\displaystyle
E=t\int\mathcal{E}(\mathbf{m})d^{2}\mathbf{r}=t\int\left(\mathcal{E}_{\mathrm{ex}}+\mathcal{E}_{\mathrm{d}}+\mathcal{E}_{\mathrm{k}}+\mathcal{E}_{\mathrm{Z}}\right)d^{2}\mathbf{r}.$
(12)

The energy now includes the exchange energy

$\displaystyle\mathcal{E}_{\mathrm{ex}}$
$\displaystyle=A\left|\nabla\mathbf{m}\right|^{2}=A\left(\left|\nabla\Theta\right|^{2}+\sin^{2}\Theta\left|\nabla\Phi\right|^{2}\right)$
(13)

where $A$ is the Heisenberg exchange constant.

Equation 12 estimates the dipolar interaction using the magnetostatic tensor
of uniformly magnetized disks. It equals the dipolar interaction between a
uniformly magnetized disk and a small differential volume, $td^{2}r$, with the
same magnetization. As such, only the surface charges are accounted for, at
the top and bottom surfaces of the disk, and at the in-plane edges. This rough
approximation neglects the role of bulk magnetic charges but is a necessary
step towards a local energy functional amenable to analysis with variational
calculus. A rigorous justification is beyond the scope of this paper, but the
leading order corrections to this approach are usually very small (see, e.g.,
Ref. [34, 35, 36, 30]).

To obtain non-dimensional expressions we rescale the in-plane spatial
coordinates ($\mathbf{r}=l_{\mathrm{ex}}\mathbf{\tilde{r}}$) by an exchange
length,

$l_{\mathrm{ex}}=\sqrt{\frac{2A}{\mu_{0}M_{s}^{2}(Q-\zeta)}};$ (14)

and the energies with a characteristic value, $E_{0}=At$;

The energy density rescaling factor can now be written as
$\mathcal{E}_{0}=E_{0}/(l_{\mathrm{ex}}^{2}t)$ and the dimensionless energy
density ($\tilde{\mathcal{E}}=\mathcal{E}/\mathcal{E}_{0}$) becomes:

$\tilde{\mathcal{E}}=\left|\tilde{\nabla}\Theta\right|^{2}+\sin^{2}\Theta\left|\tilde{\nabla}\Phi\right|^{2}-\cos^{2}\Theta-2h\sin\Theta\cos\Phi.$
(15)

Notice that, as a consequence of our local approximation for the magnetostatic
energy, the exchange length (Eq. 14) depends on the disk radius. The implicit
definition is not a problem for computation because $\zeta(R)$ is monotonic,
and therefore, $l_{\mathrm{ex}}(\zeta(R))$ is single-valued. The shape-induced
effects of the local approximation into the effective anisotropy are
incorporated in this definition to smooth the transition between the macrospin
and the wall-mediated models.

We will assume that $\Phi$ is constant ($\Phi=0$) and the profile $\Theta$
only depends on $x$.

$\tilde{E}(\Theta(x))=\int\left[\left|\frac{\partial\Theta}{\partial\tilde{x}}\right|^{2}-\cos^{2}\Theta-2h\sin\Theta\right]d^{2}\mathbf{\tilde{r}}.$
(16)

Here, we maintain a two-dimensional measure of integration to account for the
circular shape of our device, even though the $y$ coordinate is ignored when
performing variational calculus. This is the key approximation in our model.
We solve for the transversal profile of an infinite stripe and assume that it
is still valid for a circle.

The configurations $\Theta(x)$ that correspond to energy minima or saddle
points are extremes of this function which can be found from the corresponding
Euler-Lagrange equation and its accompanying boundary condition (see App. C).

$\frac{\partial^{2}\Theta}{\partial\tilde{x}^{2}}=\cos\Theta\left[\sin\Theta-h\right],\qquad\left.\frac{\partial\Theta}{\partial\tilde{x}}\right\|_{\tilde{x}=\pm\tilde{R}}=0.$
(17)

A non-uniform solution that satisfies Eq. 17 exists and corresponds to a
reversal driven by a domain wall that moves across the disk (detailed
derivation in Appendix D). The key observation is that the saddle state occurs
when the domain wall is at the center of the disk, and the magnetization
profile described by:

$\Theta_{0}(\tilde{x})\equiv\frac{\pi}{2}-2\tan^{-1}\left[\sqrt{m}\
\mathcal{B}\
\mathrm{sn}\left(\left.\frac{\tilde{x}}{\mathcal{A}}\right|m\right)\right].$
(18)

In the equation, $\mathrm{sn}(\cdot|m)$ is the Jacobi Elliptic function with
parameter $m$ chosen to satisfy the boundary conditions (not to be confused
with the magnetization magnitude $|\mathbf{m}|$=1). Appendix B explains how to
obtain the parameters used in the description of these profiles, i.e.,
$\mathcal{A}$, $\mathcal{B}$, and $m$.

Numerically integrating the energy density (Eq. 15) over the surface of the
disk,

$E=At\int_{-\tilde{R}}^{\tilde{R}}2\sqrt{\tilde{R}^{2}-\tilde{x}^{2}}\tilde{\mathcal{E}}(\Theta_{0}(x))d\tilde{x},$
(19)

gives the energy of the non-uniform saddle state. The energy barrier is
therefore:

$E_{b}=E-At\pi\tilde{R}^{2}\tilde{\mathcal{E}}_{1}.$ (20)

###  IV.1  String Method Calculations

The String method [26, 25] is a numerical procedure for calculating transition
energies and paths within the context of large fluctuations and rare events
and it is especially useful to find the minimum energy path (MEP) connecting
two metastable configurations. Practically, it is a chain-of-states algorithm
[37] that allows for precise estimates of $E_{b}$ when analytical solutions of
the saddle state cannot be obtained. To use it, we provide a guess for the
optimal escape trajectory in configuration space (called a string). Each of
the points along this path is a micromagnetic configuration of the disk. If
this path is parametrized by the normalized arc-length in magnetization space,
the location of each configuration in this path is described by a number from
0 to 1 and is referred to as the reaction coordinate.

Figure 4:  Example of energy values along the initial and relaxed string.
Configurations along the initial string are shown as snapshots in the upper
part of the figure. The string magnetization configurations, after relaxation,
are shown at the bottom of the figure. The chosen color palette is linked to
the $z$ component of the magnetization (blue for $m_{z}=-1$; white for
$m_{z}=0$; and red for $m_{z}=+1$). The saddle configuration is at the center
and the two metastable states are far left and far right.

We have used our OOMMF-based implementation of this method to find the
transition states of the nanodisk at different applied fields. An example of
this calculation is shown in Fig. 4. The initial trajectory consisted of 100
configurations that contained domains of zero width ($m_{z}$ changed abruptly
from +1 to -1 between neighboring cells). The $x$ coordinate of the wall
center varies from the left to right of the disk as the reaction coordinate
changes from 0 to 1. The top row of thumbnails contains the spin
configurations of the initial string.

As we can see in Fig. 4, the energy of the magnetic configurations along the
string starts as a staircase or steplike curve but quickly relaxes to a smooth
curve with a single maximum in $\sim 250$ iterations. At the bottom of Fig. 4,
we present the relaxed spin configuration at different reaction coordinates,
which shows the domain wall motion that mediates the switching. The individual
configurations are slightly curved Néel walls that sweep the disk from left to
right.

After the relaxation of the string, the transition between the two domains
became smoother with a profile that closely approximates our theoretical
predictions (Eq. 18), Fig. 5. The String method data is plotted for each cell
(symbols) and the scatter measures the deviation of a full micromagnetic
calculation from our one-dimensional model. In OOMMF the magnetization profile
is allowed to vary in both spatial directions.

Figure 5:  Profile of selected configurations for a $60\text{\,}\mathrm{nm}$
device. For low fields ($0.2\text{\,}\mathrm{T}$), the prediction from the
domain wall mediated model and the String method simulations match closely. At
larger fields ($0.6\text{\,}\mathrm{T}$), a deviation is visible at the edges
of the device, which can be attributed to the fact that our boundary
conditions are strictly valid only in infinitely long stripes. Points show the
magnetization of each cell obtained from OOMMF. Inset shows the magnetization
configuration for the $0.2\text{\,}\mathrm{T}$ case.

The re-scaled field $\mathbf{h}$ varies between 0 and 1 in our model, and as
such determines a critical field, $B_{c}(R)=\mu_{0}M_{s}(Q-\zeta)$, at which
the system becomes saturated in-plane and the bi-stability is lost. As the
magnet’s aspect ratio changes from a high needle-like cylinder
($\zeta\rightarrow-1/2$) to an extended film $\zeta\rightarrow 1$, the field
decreases from $B_{c}=$ $1.7\text{\,}\mathrm{T}$ to $0.61\text{\,}\mathrm{T}$.
Since film devices have small aspect ratios, the fields of interest are all
below $0.7\text{\,}\mathrm{T}$.

##  V Results and Discussion

Since the energy barrier height greatly influences the rate of fluctuation and
this barrier can be changed by applying a transverse field to the device, the
following section examines the effect of the field on the energy barrier for
several disk diameters.

###  V.1 Size and Field Dependence

In Fig. 6(a), we show the energy barrier change with the applied field for
different device sizes according to the domain wall (DW) mediated models and
the String method calculations. In general, the barrier decreases with
increasing field and decreasing diameter. There is, however, a deviation from
this trend for high fields as expected from the dependence of the critical
field $B_{c}$ on device radius.

Overall, there is a good agreement between predictions from the DW-mediated
model and the String method results although the latter yields slightly lower
energy barriers. This is expected since our model only considers a one-
dimensional dependence of the magnetization, and the String method calculation
allows for other relaxation paths. As a rule, chain-of-states algorithms reach
transition paths that are lower in energy than their analytical counterparts
since they have larger degrees of freedom and are able to access a wider
configurational space. This explains why most of our analytical values are
above their corresponding simulation results in Fig. 6.

Figure 6:  (a) Energy barrier ($E_{b}$) vs size at different fields in a DW-
mediated transition. We notice a non-monotonicity vs. field for high fields
(i.e. $>$ $0.6\text{\,}\mathrm{T}$). (b) Scaled energy barrier vs. scaled
field. All results obtained from the macrospin model are characterized by the
top curve in red.

Another result is a non-monotonicity of the energy barrier with respect to
nanodisk size, Fig. 6(a). This is evident for the curves at large fields
($H=0.6$ and $0.7\text{\,}\mathrm{T}$) since the macrospin’s energy barrier
vanishes for $H>\mu_{0}M_{s}(Q-1)$ as $R\rightarrow\infty$.

Figure 6(b) shows the scaled energy density barrier vs. the scaled applied
field. The dashed lines capture String method calculations at different
diameters. The data are the same as those in panel (a) but this new
representation helps to easily compare the behavior of nanodisks (of different
dimensions) under transverse fields. As disk size decreases towards 10 nm
diameter, the energy barrier approaches the value predicted by the macrospin
model (Sec. III). At $10\text{\,}\mathrm{nm}$ we have reached the critical
dimension for uniform switching, in fact, the curve overlaps perfectly with
the macrospin even at zero applied field. For larger disk diameters, a
stronger field needs to be employed to have an agreement ($h=0.4$ for a
$20\text{\,}\mathrm{nm}$ device). The larger disks considered here ($D=$ 40
and $60\text{\,}\mathrm{nm}$), do not show macrospin behavior. They exhibit a
lower energy barrier for all the applied fields.

###  V.2 Tunneling Magnetoresistance Dependence on Field

An expected effect of the application of a transverse field to an MTJ is a
decrease in tunneling magnetoresistance (TMR). To predict the TMR vs applied
field behavior, we have assumed a fixed reference layer (blue layer in the
inset of Fig. 7).

The tunneling conductance is assumed to depend on $\theta$ as follows (e.g.,
[38]):

$G(\theta)=\frac{1}{2}(G_{P}+G_{AP})+\frac{1}{2}(G_{P}-G_{AP})\cos\theta$ (21)

and the tunneling magnetoresistance is the fraction:

$T=\frac{G_{P}-G_{AP}}{G_{AP}}$ (22)

In our case, from Eq. 9, the angle is defined as
$\Theta_{\uparrow}=\sin^{-1}h$ and
$\Theta_{\downarrow}=\pi-\Theta_{\uparrow}$.

A nice way to see the reduction in TMR is to relate $\mathcal{T}(h)$ to $T$,
the TMR as a function of the transverse field to the zero field value (i.e.,
that given by Eq. 22). Then the maximum device TMR is given by:

$\mathcal{T}(h)=\frac{2T\sqrt{1-h^{2}}}{T(1-\sqrt{1-h^{2}})+2}.$ (23)

Figure 7 shows a plot of the TMR as a function of the transverse field $h$
with $T=1.0$, an initial TMR of 100 %.

Figure 7:  Tunneling magnetoresistance $\mathcal{T}(h)$ vs field $h$. The
magnetoresistance is reasonably large even for large transverse fields, for
fields at which the energy barrier can be reduced significantly. In the inset,
for the reader’s convenience, parallel and antiparallel states schematic in a
pMTJ.

###  V.3 Intermediate Metastable State

Figure 8:  Energy values along relaxed strings for different device sizes with
an applied field $0.7\text{\,}\mathrm{T}$. From $D=$ $60\text{\,}\mathrm{nm}$
a new metastable state appears at the center of the string and it cannot be
predicted by our model. Inset shows the magnetization configuration for $D=$
$100\text{\,}\mathrm{nm}$ and at a reaction coordinate of 0.5.

In our investigation, a new intermediate state appears for large devices
($D\geq$ $60\text{\,}\mathrm{nm}$), Fig. 8. This new metastable state is
higher in energy with respect to the endpoint reaction coordinates. The
configuration at this new local energy minimum is highly non-uniform and
cannot be predicted from our one-dimensional model that features a narrow
domain wall only. We will just provide a heuristic explanation about what may
stabilize the metastable state by considering the preferred magnetization
direction of an element of volume $td^{2}r$ depending on its position inside
the disk.

The inset of Fig. 8 shows the magnetization configuration for $D=$
$100\text{\,}\mathrm{nm}$ and at reaction coordinate of 0.5. To understand
this case, let us first consider the magnetization at the rim of the disk,
under the applied field $H\mathbf{\hat{x}}$. At position $x=\pm R$, there is
an out-plane buckling of the magnetization as if the effective perpendicular
anisotropy became enhanced in the vicinity of the in-plane edge. This out-of-
plane orientation decreases the crystalline anisotropy energy at the expense
of increasing the Zeeman and magnetostatic energies.

Next, we consider the area closer to the disk center. The penalty for out-of-
plane magnetization orientations is stronger in the interior of the disk than
at the edges. In the presence of a strong field along $\mathbf{x}$ direction,
this magnetization direction becomes favorable. By enlarging the area with
close to uniform magnetization, the exchange energy penalty is reduced.

Putting together all these observations, we now provide an explanation for the
stability of the newly found configuration. The field is strong, so both
metastable directions are slightly out of plane. The perpendicular anisotropy
is stronger at the disk edges than in the interior. As a result, the
magnetization aligns with the field for most points except within a narrow rim
near the disk edge where the contribution to the exchange energy is
concentrated.

##  VI Conclusions

In summary, we have compared two possible transition pathways for thermally
activated switching, in perpendicularly magnetized nanodisks. One pathway is a
coherent rotation of the magnetization (macrospin) and another transition is
domain wall-mediated switching with spatially non-uniform magnetization
states. These transition paths were tested using overdamped micromagnetic and
String method simulations respectively. Specifically, we have compared energy
barriers and the spatial profiles of analytical models to those obtained
numerically resulting in minimal discrepancies, and only in extreme cases.
Moreover, a new metastable state in the energy barrier behavior was found in
large devices ($D\geq$ $60\text{\,}\mathrm{nm}$) for very high fields.

From an application perspective, the aim of this work was to prove the
possibility of using a transverse field to lower the energy barrier of pMTJ
devices with different nanodisk diameters. Both the models and the
simulations, show that the energy barrier of a perpendicularly magnetized
nanodisk can be tuned with the applied field. A quick estimate of the
switching rate can be done assuming $\Gamma_{0}=$ $1\text{\,}\mathrm{GHz}$ as
a characteristic frequency of magnetic excitations. In this case, a barrier
$E_{b}=4.6~{}\mathrm{k_{B}T}$ results in switching times in the order of
$100\text{\,}\mathrm{ns}$. For a device with a zero field barrier of
$E_{b}=20~{}\mathrm{k_{B}T}$ (i.e. a $12\text{\,}\mathrm{nm}$ diameter device,
with a saturation field of $B_{\mathrm{sat}}=$ $0.98\text{\,}\mathrm{T}$) an
applied field of $0.51\text{\,}\mathrm{T}$ ($h=0.52$) would produce this
fluctuation rate. Similarly, a $40~{}\mathrm{k_{B}T}$ device
($18\text{\,}\mathrm{nm}$ diameter,
$B_{\mathrm{sat}}=$$0.9\text{\,}\mathrm{T}$) achieves the above rate with a
field of $0.60\text{\,}\mathrm{T}$ ($h=0.67$). Thus, as we have illustrated in
Fig. 7, the TMR would still be significant, enabling state determination. The
findings confirm that pMTJ devices, whose fluctuation rates are determined by
a transverse magnetic field, can produce rapid and measurable fluctuations,
presenting exciting possibilities for their application in generating genuine
random numbers.

What still remains to be explored are the optimal device dimension, the
testing of our theory, and simulations against experimental device
measurements. Moreover, an accurate estimation of the prefactor in the attempt
rate expression would be needed, to fully assess the stochasticity of the
device.

## Acknowledgments

The research at NYU was supported by the DOE Office of Science (ASCR/BES)
Microelectronics Co-Design project COINFLIPS. This work was partially
supported by the Swedish Research Council (VR), project Fundamental
Fluctuations in Spintronics, 2017-04196. We also thank the funding agencies
Nils and Hans Backmark Foundation (J-2021-2437) and Karl Engvers Foundation
for supporting the project. The computations were partially enabled by
resources provided by the Swedish National Infrastructure for Computing (SNIC)
at High Performance Computing Center North (HPC2N) partially funded by the
Swedish Research Council through grant agreement no. 2018-05973. Gabriel D.
Chaves-O’Flynn was supported in part by the National Science Centre Poland
under OPUS funding Grant No. 2019/33/B/ST5/02013.

##  Appendix A Demagnetizing factor for a nanodisk.

The demagnetization tensor coefficient $N_{zz}$ for a uniformly magnetized
thin disk can be calculated [32] using hypergeometric functions [39, 40]

$N_{zz}=1+\frac{8R}{3\pi
t}-{}_{2}F{{}_{1}}\left[-\frac{1}{2},\frac{1}{2};2,-\left(\frac{2R}{t}\right)^{2}\right];$
(24)

or obtained from any standard micromagnetic solver from the magnetostatic
energies of mutually orthogonal saturated states, as was done in [28]. It is
easily checked that the values obtained with both approaches are identical.

##  Appendix B Profile parameters

The coefficients $\mathcal{A}$ and $\mathcal{B}$ can be obtained from $m$ and
$h$ (ignoring nonphysical solutions) as follows:

$\displaystyle\mathcal{A}(m,h)$
$\displaystyle=\frac{1-m}{\sqrt{1+m-\sqrt{h^{2}(1-m)^{2}+4m}}}$ (25)
$\displaystyle\sin(\Theta_{R})$
$\displaystyle=h+\frac{m+1}{m-1}+\sqrt{h^{2}+\frac{4m}{(m-1)^{2}}}$ (26)
$\displaystyle\beta$
$\displaystyle\equiv\frac{1}{2}\left(\frac{\pi}{2}-\Theta\right);\quad\beta_{R}=\frac{1}{2}\left(\frac{\pi}{2}-\Theta_{R}\right)$
(27) $\displaystyle\mathcal{B}(h,\Theta_{R})$
$\displaystyle=\tan\left(\beta_{R}\right)$ $\displaystyle\mathcal{B}(m,h)$
$\displaystyle=\sqrt{\frac{2-h(1-m)+\sqrt{h^{2}(m-1)^{2}+4m}}{h(1-m)-2m-\sqrt{h^{2}(m-1)^{2}+4m}}}$
(28)

The configuration described by Eq. 18 has a wavelength equal to
$\lambda=4\mathbf{K}(m)$; to satisfy the boundary conditions, the disk
diameter has to be a multiple of half this length $D=n\lambda/2$. We then have
an additional relation for the non-uniform saddle
$\tilde{R}=n\mathcal{A}(m,h)\mathbf{K}(m)$ for which $m$ can be found
numerically. Although large disks will be able to sustain a hierarchy of
states with multiple domain walls, we will restrict our analysis to the
simplest case $(n=1)$. As part of the supplementary materials, we provide a
Mathematica notebook to obtain $m,\mathcal{A}$ and $\mathcal{B}$ for specific
geometric and material parameters as well as a description of the necessary
algebraic steps to obtain the above explicit dependence of $\mathcal{A}$ and
$\mathcal{B}$ on $m$ and $h$ [supplementary material - S1].

##  Appendix C Derivation of the Boundary Conditions

This is obtained from variational calculus:

$\delta\tilde{E}=\iint\delta\left[\left|\frac{\partial\Theta}{\partial\tilde{x}}\right|^{2}-\cos^{2}\Theta-2h\sin\Theta\right]d\tilde{x}d\tilde{y}.$
(29)

To simplify our work, in the text we solve the equation at $y=0$ and assume
the same profile is valid for the rest of the disk. In reality, this is the
solution for a square of side $2R$. We assume the discrepancy with the
solution for a disk will be small

$2\iint\left[\frac{\partial\Theta}{\partial\tilde{x}}\frac{\partial\delta\Theta}{\partial\tilde{x}}+\left[\sin\Theta\cos\Theta-h\cos\Theta\right]\delta\Theta\right]d\tilde{x}d\tilde{y}.$
(30)

Using integration by parts in the first term we obtain:

$2\iint\left[-\frac{\partial^{2}\Theta}{\partial\tilde{x}^{2}}+\sin\Theta\cos\Theta-h\cos\Theta\right]\delta\Theta
d\tilde{x}d\tilde{y}.$ (31)
$+2\int^{+\tilde{R}}_{-\tilde{R}}\left[\frac{\partial\Theta}{\partial\tilde{x}}\right]^{\tilde{R}}_{-\tilde{R}}\delta\Theta
d\tilde{y}$ (32)

An extremal satisfies $\delta E=0$ for arbitrary $\delta\Theta$; therefore,
the first of these integrals provides the differential equation to be solved;
the second integral, the boundary conditions to be satisfied
$\left.\frac{\partial\Theta}{\partial\tilde{x}}\right\|_{\tilde{x}=\pm\tilde{R}}=0$.

##  Appendix D Derivation of the non-constant saddle profile.

We postulate that, in addition to the Neumann boundary conditions in Eq.17,
the profiles should also satisfy the following symmetry:

$\Theta(\tilde{x}=-\tilde{R})=\Theta_{R}\leq\frac{\pi}{2},\quad\Theta(\tilde{x}=+\tilde{R})=\pi-\Theta_{R}.$
(33)

Which is motivated on physical grounds. It is reasonable to expect that the
non-uniform extremal is made of two energetically equivalent orientations
separated by a transition region in between. We expect that for large disks
the magnetization fully rotates between $\Theta_{\uparrow}$ and
$\Theta_{\downarrow}$, i.e.
$\lim_{R\rightarrow\infty}\Theta_{R}=\Theta_{\uparrow}$.

Multiplying the differential equation in 17 by $d\Theta/d\tilde{x}$ and
integrating in $x$ results in the equivalent equation:

$\frac{1}{2}\left(\frac{d\Theta}{d\tilde{x}}\right)^{2}=\frac{1}{2}\sin^{2}\Theta-h\sin\Theta+C_{1}.$
(34)

We solve for the first constant of integration, $C_{1}$, after substitution of
the boundary conditions.

$C_{1}=\sin\Theta_{R}\left[h-\frac{1}{2}\sin\Theta_{R}\right]\qquad\lim_{R\rightarrow\infty}C_{1}=\frac{h^{2}}{2}.$
(35)

We use $m_{x}=\sin(\Theta)$ to rewrite Eq.34 as an elliptical integral:

$d\tilde{x}=\frac{dm_{x}}{\sqrt{1-m_{x}^{2}}\sqrt{m_{x}^{2}-2hm_{x}+2C_{1}}}.$
(36)

After integrating we obtain an expression involving, $F(\tilde{x}|m)$, the
Elliptic Integral of the first kind:

$\displaystyle\tilde{x}+C_{2}$ $\displaystyle=\mathcal{A}\cdot
F\left(\left.\sin^{-1}\left(\mathcal{B}\cot\left(\beta\right)\right)\right)|m\right)$
(37)

Solving Eq.(37) for $\Theta$, we get an expression involving a Jacobi elliptic
function, $\mathrm{sn}(x|m)$:

$\Theta=\frac{\pi}{2}-2\cot^{-1}\left[\frac{1}{\mathcal{B}}\mathrm{sn}\left(\left.\frac{\tilde{x}+C_{2}}{\mathcal{A}}\right|m\right)\right].$
(38)

Since we expect $\Theta(x=0)=\pi/2$, we can find $C_{2}$ that satisfies:

$\mathrm{sn}\left(\left.\frac{C_{2}}{\mathcal{A}}\right|m\right)=\infty.$ (39)

Since the first pole of $\mathrm{sn}(x|m)$ is located at $i\mathbf{K}(1-m)$,
we find $C_{2}=i\mathcal{A}\mathbf{K}(1-m)$. This allows us to use a change of
variables and other identities of Jacobi elliptic functions and simplify this
to Eq. 18

In the large disk limit ($m$=1), the above parameters can be simplified:

$\displaystyle\lim_{R\rightarrow\infty}\mathcal{A}$
$\displaystyle=\sqrt{\frac{4}{1-h^{2}}}$ (40)
$\displaystyle\lim_{R\rightarrow\infty}\mathcal{B}$
$\displaystyle=\sqrt{\frac{1-h}{1+h}}=\frac{1}{\cot{\beta_{\mathrm{max}}}}$
(41) $\displaystyle\beta_{\mathrm{max}}$
$\displaystyle\equiv\frac{1}{2}\left(\frac{\pi}{2}-\Theta_{\uparrow}\right)$
(42)

and, when $R\rightarrow\infty$,

$\Theta_{0}=\frac{\pi}{2}-2\tan^{-1}\left[\tan\beta_{\mathrm{max}}\tanh\left(x\sqrt{\frac{1-h^{2}}{4}}\right)\right].$
(43)

In this form, the solution provides for a quick intuitive interpretation of
this result. Since $\tanh(x)$ changes sigmoidally from -1 to 1 as x grows from
$-\infty$ to $+\infty$, the angle $\Theta$ acquires the same sigmoidal
dependence but varies between $\Theta_{\uparrow}$ and $\Theta_{\downarrow}$.

## References

  * Darwish  _et al._ [2020] A. Darwish, A. E. Hassanien, and S. Das, A survey of swarm and evolutionary computing approaches for deep learning, Artificial Intelligence Review 53, 1767 (2020).
  * Shastri  _et al._ [2021] B. J. Shastri, A. N. Tait, T. Ferreira de Lima, W. H. P. Pernice, H. Bhaskaran, C. D. Wright, and P. R. Prucnal, Photonics for artificial intelligence and neuromorphic computing, Nature Photonics 15, 102 (2021), number: 2 Publisher: Nature Publishing Group.
  * Hoffmann  _et al._ [2022] A. Hoffmann, S. Ramanathan, J. Grollier, A. D. Kent, M. J. Rozenberg, I. K. Schuller, O. G. Shpyrko, R. C. Dynes, Y. Fainman, A. Frano, E. E. Fullerton, G. Galli, V. Lomakin, S. P. Ong, A. K. Petford-Long, J. A. Schuller, M. D. Stiles, Y. Takamura, and Y. Zhu, Quantum materials for energy-efficient neuromorphic computing: Opportunities and challenges, APL Materials 10, 070904 (2022).
  * Grollier  _et al._ [2020] J. Grollier, D. Querlioz, K. Y. Camsari, K. Everschor-Sitte, S. Fukami, and M. D. Stiles, Neuromorphic spintronics, Nature Electronics 3, 360 (2020), number: 7 Publisher: Nature Publishing Group.
  * Tanaka  _et al._ [2019] G. Tanaka, T. Yamane, J. B. Héroux, R. Nakane, N. Kanazawa, S. Takeda, H. Numata, D. Nakano, and A. Hirose, Recent advances in physical reservoir computing: A review, Neural Networks 115, 100 (2019).
  * Liu  _et al._ [2021] Y. Liu, S. Liu, Y. Wang, F. Lombardi, and J. Han, A survey of stochastic computing neural networks for machine learning applications, IEEE Trans. on Neural Net. and Learning Sys. 32, 2809 (2021).
  * Qu  _et al._ [2018] Y. Qu, B. F. Cockburn, Z. Huang, H. Cai, Y. Zhang, W. Zhao, and J. Han, Variation-Resilient True Random Number Generators Based on Multiple STT-MTJs, IEEE Transactions on Nanotechnology 17, 1270 (2018).
  * Bassham  _et al._ [2010] L. E. Bassham, A. L. Rukhin, J. Soto, J. R. Nechvatal, M. E. Smid, S. D. Leigh, M. Levenson, M. Vangel, N. A. Heckert, and D. L. Banks, A Statistical Test Suite for Random and Pseudorandom Number Generators for Cryptographic Applications, NIST (2010).
  * Rehm  _et al._ [2023] L. Rehm, C. C. M. Capriata, S. Misra, J. D. Smith, M. Pinarbasi, B. G. Malm, and A. D. Kent, Stochastic Magnetic Actuated Random Transducer Devices Based on Perpendicular Magnetic Tunnel Junctions, Physical Review Applied 19, 024035 (2023).
  * Augustine  _et al._ [2010] C. Augustine, A. Raychowdhury, D. Somasekhar, J. Tschanz, K. Roy, and V. K. De, Numerical analysis of typical stt-mtj stacks for 1t-1r memory arrays, in  _2010 International Electron Devices Meeting_ (2010) pp. 22.7.1–22.7.4.
  * Pajouhi  _et al._ [2015] Z. Pajouhi, X. Fong, and K. Roy, Device/circuit/architecture co-design of reliable stt-mram, in  _2015 Proc. Des. Autom. Test Europe Conf. Exhib._ (2015) pp. 1437–1442.
  * Puebla  _et al._ [2020] J. Puebla, J. Kim, K. Kondou, and Y. Otani, Spintronic devices for energy-efficient data storage and energy harvesting, Communications Materials 1, 1 (2020).
  * Sun and Safranski [2022] J. Z. Sun and C. Safranski, Metrology and metrics for spin-transfer-torque switched magnetic tunnel junctions in memory applications, J. of Magn. and Magn. Mat. 563, 169878 (2022).
  * Vodenicarevic  _et al._ [2017] D. Vodenicarevic, N. Locatelli, A. Mizrahi, J. S. Friedman, A. F. Vincent, M. Romera, A. Fukushima, K. Yakushiji, H. Kubota, S. Yuasa, S. Tiwari, J. Grollier, and D. Querlioz, Low-Energy Truly Random Number Generation with Superparamagnetic Tunnel Junctions for Unconventional Computing, Physical Review Applied 8, 054045 (2017).
  * Parks  _et al._ [2018] B. Parks, M. Bapna, J. Igbokwe, H. Almasi, W. Wang, and S. A. Majetich, Superparamagnetic perpendicular magnetic tunnel junctions for true random number generators, AIP Advances 8, 055903 (2018).
  * Kaiser  _et al._ [2019] J. Kaiser, A. Rustagi, K. Y. Camsari, J. Z. Sun, S. Datta, and P. Upadhyaya, Subnanosecond Fluctuations in Low-Barrier Nanomagnets, Physical Review Applied 12, 054056 (2019).
  * Hayakawa  _et al._ [2021] K. Hayakawa, S. Kanai, T. Funatsu, J. Igarashi, B. Jinnai, W. A. Borders, H. Ohno, and S. Fukami, Nanosecond Random Telegraph Noise in In-Plane Magnetic Tunnel Junctions, Physical Review Letters 126, 117202 (2021).
  * Safranski  _et al._ [2021] C. Safranski, J. Kaiser, P. Trouilloud, P. Hashemi, G. Hu, and J. Z. Sun, Demonstration of Nanosecond Operation in Stochastic Magnetic Tunnel Junctions, Nano Letters 21, 2040 (2021).
  * Hänggi  _et al._ [1990] P. Hänggi, P. Talkner, and M. Borkovec, Reaction-rate theory: fifty years after Kramers, Reviews of Modern Physics 62, 251 (1990).
  * Braun [1994] H.-B. Braun, Fluctuations and instabilities of ferromagnetic domain-wall pairs in an external magnetic field, Phys. Rev. B 50, 16485 (1994).
  * Borders  _et al._ [2019] W. A. Borders, A. Z. Pervaiz, S. Fukami, K. Y. Camsari, H. Ohno, and S. Datta, Integer factorization using stochastic magnetic tunnel junctions, Nature 573, 390 (2019).
  * Camsari  _et al._ [2021] K. Y. Camsari, M. M. Torunbalci, W. A. Borders, H. Ohno, and S. Fukami, Double-Free-Layer Magnetic Tunnel Junctions for Probabilistic Bits, Physical Review Applied 15, 044049 (2021).
  * Chen  _et al._ [2022] X. Chen, J. Zhang, and J. Xiao, Magnetic-tunnel-junction-based true random-number generator with enhanced generation rate, Phys. Rev. Applied 18, L021002 (2022).
  * Garanin  _et al._ [1999] D. A. Garanin, E. C. Kennedy, D. S. F. Crothers, and W. T. Coffey, Thermally activated escape rates of uniaxial spin systems with transverse field: Uniaxial crossovers, Physical Review E 60, 6499 (1999).
  * Weinan  _et al._ [2002] E. Weinan, W. Ren, and E. Vanden-Eijnden, String method for the study of rare events, Physical Review B 66, 052301 (2002).
  * Weinan  _et al._ [2007] E. Weinan, W. Ren, and E. Vanden-Eijnden, Simplified and improved string method for computing the minimum energy paths in barrier-crossing events, The Journal of Chemical Physics 126, 164103 (2007).
  * Donahue [1999] M. J. Donahue, OOMMF User’s Guide, Version 1.0, NIST (1999), last Modified: 2018-11-10T10:11-05:00 Publisher: Michael J. Donahue.
  * Chaves-O’Flynn  _et al._ [2015] G. D. Chaves-O’Flynn, G. Wolf, J. Z. Sun, and A. D. Kent, Thermal Stability of Magnetic States in Circular Thin-Film Nanomagnets with Large Perpendicular Magnetic Anisotropy, Physical Review Applied 4, 024010 (2015).
  * Chaves-O’Flynn  _et al._ [2013] G. D. Chaves-O’Flynn, E. Vanden-Eijnden, D. L. Stein, and A. D. Kent, Energy barriers to magnetization reversal in perpendicularly magnetized thin film nanomagnets, Journal of Applied Physics 113, 023912 (2013).
  * Martens  _et al._ [2006] K. Martens, D. L. Stein, and A. D. Kent, Magnetic reversal in nanoscopic ferromagnetic rings, Physical Review B 73, 054413 (2006).
  * Chaves-O’Flynn  _et al._ [2009] G. D. Chaves-O’Flynn, A. D. Kent, and D. L. Stein, Micromagnetic study of magnetization reversal in ferromagnetic nanorings, Physical Review B 79, 184421 (2009).
  * Beleggia  _et al._ [2006] M. Beleggia, M. D. Graef, and Y. T. Millev, The equivalent ellipsoid of a magnetized body, Journal of Physics D: Applied Physics 39, 891 (2006).
  * Rehm  _et al._ [2019] L. Rehm, G. Wolf, B. Kardasz, M. Pinarbasi, and A. D. Kent, Sub-nanosecond spin-torque switching of perpendicular magnetic tunnel junction nanopillars at cryogenic temperatures, Applied Physics Letters 115, 182404 (2019).
  * [34] R. V. Kohn and V. V. Slastikov, Another thin-film limit of micromagnetics, Archive for Rational Mechanics and Analysis 178, 227.
  * [35] G. Gioia and R. D. James, Micromagnetics of very thin films, Proc. of the Royal Society of London. Series A 453, 213.
  * [36] C. B. Muratov and V. V. Slastikov, Domain structure of ultrathin ferromagnetic elements in the presence of dzyaloshinskii–moriya interaction, Proceedings of the Royal Society A 473, 20160666.
  * [37] P. F. Bessarab, V. M. Uzdin, and H. Jónsson, Method for finding mechanism and activation energy of magnetic transitions, applied to skyrmion and antivortex annihilation, Computer Physics Communications 196, 335.
  * Zhu and Park [2006] J.-G. J. Zhu and C. Park, Magnetic tunnel junctions, Materials Today 9, 36 (2006).
  * Abramowitz and Stegun [1965] M. Abramowitz and I. A. Stegun, eds.,  _Handbook of Mathematical Functions: With Formulas, Graphs, and Mathematical Tables_ , 0009th ed. (Dover Publications, New York, NY, 1965).
  * Armitage and Eberlein [2006] J. V. Armitage and W. F. Eberlein,  _Elliptic Functions_, London Mathematical Society Student Texts (Cambridge University Press, Cambridge, 2006).