# Stability of the Néel quantum critical point in the presence of Dirac
fermions

Huanzhi Hu  London Centre for Nanotechnology, University College London,
Gordon St., London, WC1H 0AH, United Kingdom  Jennifer Lin  London Centre for
Nanotechnology, University College London, Gordon St., London, WC1H 0AH,
United Kingdom  Mikolaj D. Uryszek  London Centre for Nanotechnology,
University College London, Gordon St., London, WC1H 0AH, United Kingdom  Frank
Krüger  London Centre for Nanotechnology, University College London, Gordon
St., London, WC1H 0AH, United Kingdom  ISIS Facility, Rutherford Appleton
Laboratory, Chilton, Didcot, Oxfordshire OX11 0QX, United Kingdom

###### Abstract

We investigate the stability of the Néel quantum critical point of two-
dimensional quantum antiferromagnets, described by a non-linear $\sigma$ model
(NL$\sigma$M), in the presence of a Kondo coupling to $N_{f}$ flavours of two-
component Dirac fermion fields. The long-wavelength order parameter
fluctuations are subject to Landau damping by electronic particle-hole
fluctuations. Using momentum-shell RG, we demonstrate that the Landau damping
is weakly irrelevant at the Néel quantum critical point, despite the fact that
the corresponding self-energy correction dominates over the quadratic gradient
terms in the IR limit. In the ordered phase, the Landau damping increases
under the RG, indicative of damped spin-wave excitations. Although the Kondo
coupling is weakly relevant, sufficiently strong Landau damping renders the
Néel quantum critical point quasi-stable for $N_{f}\geq 4$ and
thermodynamically stable for $N_{f}<4$. In the latter case, we identify a new
multi-critical point which describes the transition between the Néel critical
and Kondo run-away regimes. The symmetry breaking at this fixed point results
in the opening of a gap in the Dirac fermion spectrum. Approaching the multi-
critical point from the disordered phase, the fermionic quasiparticle residue
vanishes, giving rise to non-Fermi-liquid behavior.

##  I Introduction

The discovery of topological insulators [1, 2] has initiated an explosion of
research into Dirac or Weyl semimetals and topological aspects of electronic
bandstructures [3, 4, 5, 6]. Dirac fermions with relativistic dispersion
around point-like Fermi surfaces can arise as low-energy excitations of weakly
interacting electron systems. The most prominent example is graphene [7],
which can be described by a tight-binding model of electrons on the half-
filled honeycomb lattice.

Because of their point-like Fermi surfaces nodal semi- metals provide the
simplest setting to study fermionic quantum criticality. While Dirac
semimetals are stable against weak repulsive interactions, a consequence of
the vanishing density of states at the Fermi level, sufficiently strong short-
range interactions can give rise to a range of competing instabilities. For
the extended Hubbard model on the half-filled honeycomb lattice rich phase
diagrams were established [8, 9, 10, 11, 12, 13, 14, 15, 16, 17], showing
antiferromagnetic, charge ordered, Kekule and topological Haldane phases.
Sufficiently strong on-site Hubbard repulsion gives rise to a transition to an
antiferromagnetic state with a gap in the electron spectrum that is
proportional to the staggered magnetization.

Since the fermionic particle-hole excitations are gapless at such quantum
phase transitions, the critical behavior falls outside the Landau-Ginzburg-
Wilson paradigm of a pure order parameter description [18]. Instead, the
nature of the transitions can be studied within a field theory that describes
the coupling of the bosonic order parameter field, which is introduced through
a Hubbard-Stratonovich decoupling of the interaction vertex, to the gapless
Dirac fermions [19, 20, 21]. In the high-energy community, such field theories
are known as Gross-Neveu-Yukawa (GNY) theories [22, 23]. For the
antiferromagnetic transition driven by local Hubbard repulsion the staggered
magnetization is described by an $O(3)$ order-parameter field, and the field
theory usually referred to as the Heisenberg-GNY model. At the fermion-induced
GNY fixed point, the fermions acquire an anomalous dimension, resulting in the
fermion spectral functions with branch cuts rather than quasi-particle poles
[20]. Such non-Fermi liquid behavior is the hallmark of fermionic quantum
criticality.

Interesting criticality is also expected if Dirac fermions are coupled to
local magnetic moments. In the case of graphene, local moments can be
introduced by adatoms [24] or defects [25], as evidenced by the experimental
observation of a Kondo effect [26]. Moreover, the high Kondo temperatures [26]
reflect the strong coupling between the local moments and the conduction
electrons in Dirac materials [27]. In subsequent work [28] the unusual Kondo
quantum criticality of magnetic adatoms in graphene and the very fast power-
law decay of the RKKY interaction between them was established.

In this paper we consider two-dimensional quantum antiferromagnets with Kondo
coupling between the local moments and Dirac electrons. While two-dimensional
Kondo lattice models with Dirac points close to the Fermi level might be rare,
there is prospect of engineering such models in heterostructures of graphene
and two-dimensional van der Waals magnets [29], such as the honeycomb
Heisenberg antiferromagnets MnPS3 and NiPS3 [30, 31, 32, 33, 34, 35, 36].

The aim of this work is to investigate the stability of the Néel quantum
critical point of a local moment antiferromagnet against the Kondo coupling to
Dirac fermions. We will demonstrate that the universal critical behavior of
this model is different from that of the Heisenberg GNY theory which describes
the antiferromagnetic ordering transition of a purely electronic model with
local Hubbard repulsion between Dirac electrons.

The outline of this paper is as follows. In Sec. II we introduce the non-
linear $\sigma$ model (NL$\sigma$M), describing the long-wavelength behavior
of the two-dimensional quantum antiferromagnet, with Kondo coupling to $N_{f}$
Dirac fermion pairs. We discuss the importance of Landau damping of the Néel
order parameter fluctuations by low-energy electronic particle-hole
fluctuations. In Sec. III we analyze the Landau-damped NL$\sigma$M. Using
momentum-shell RG, we demonstrate that Landau damping is weakly irrelevant at
the Néel quantum critical point but increases in the ordered state, indicating
that spin-wave excitations are damped. The full set of RG equations, including
the Kondo coupling, are derived in Sec. IV and analyzed in Sec. V. We show
that while the Kondo coupling is weakly relevant at the Néel quantum critical
point, sufficiently strong Landau damping renders the critical point quasi-
stable for any realistic system size for $N_{f}\geq 4$ and thermodynamically
stable for $N_{f}<4$. In the latter case, a new multi-critical point captures
the transition between Néel critical and Kondo run-away regimes. We analyze
the universal critical behavior associated with this fixed point. As
demonstrated in Sec. VI, the behavior in $D=3$ space-time dimensions is not
accessible within an $\epsilon$-expansion above the lower critical dimension,
$D=2+\epsilon$. Finally, in Sec. VII we summarize and discuss our results.

##  II Model

Our starting model is a non-linear $\sigma$ model (NL$\sigma$M) which
describes the Néel transition of a two-dimensional quantum antiferromagnet
[37, 38]. This model is coupled to $N_{f}$ copies of two component Dirac
electrons via the conventional Kondo coupling. On a microscopic level this
model could be realized in the low-energy limit of a quantum antiferromagnet
on the honeycomb lattice with Kondo coupling to noninteracting electrons that
move on either the same or adjacent honeycomb lattice at half-filling. For
this realization we would have $N_{f}=4$ due to twofold spin and valley
degeneracies. The effective continuum field theory at zero temperature is
given by the imaginary-time path integral over the action
$S=S_{f}+S_{N}+S_{K}$, with contributions

$\displaystyle S_{f}$ $\displaystyle=$
$\displaystyle\int_{\mathbf{q},\omega}\overline{\bm{\psi}}(\mathbf{q},\omega)\left(-i\frac{\omega}{v_{F}}+q_{x}\tau_{x}+q_{y}\tau_{y}\right)\bm{\psi}(\mathbf{q},\omega),$
$\displaystyle S_{N}$ $\displaystyle=$ $\displaystyle\frac{1}{2g}\int
d^{2}\mathbf{r}\int_{0}^{\infty}d\tau\left\\{(\nabla\vec{N})^{2}+\frac{1}{c^{2}}(\partial_{\tau}\vec{N})^{2}\right\\},$
$\displaystyle S_{K}$ $\displaystyle=$
$\displaystyle\frac{\lambda}{\sqrt{N_{f}}}\int
d^{2}\mathbf{r}\int_{0}^{\infty}d\tau\overline{\bm{\psi}}\left(\vec{N}\cdot\vec{\sigma}\otimes\tau_{z}\right)\bm{\psi},$
(1)

where $S_{f}$ describes two-dimensional Dirac fermions with Fermi velocity
$v_{F}$, written in terms of fermionic Grassmann fields $\bm{\psi}$. The term
$S_{N}$ is the conventional NL$\sigma$M in terms of the staggered three-
component Néel order parameter field $\vec{N}(\mathbf{r},\tau)$ which
satisfies the constraint $\vec{N}^{2}(\mathbf{r},\tau)=1$. Here $c$ denotes
the spin-wave velocity and the coupling constant $g$ is inversely proportional
to the spin stiffness. The last contribution $S_{K}$ is the Kondo coupling
between the local moments and Dirac electrons. Here $\vec{\sigma}$ is the
vector of spin Pauli matrices, while the Pauli matrices $\tau_{\alpha}$ act on
sub-lattice space. Note that since $\vec{N}$ describes the staggered
magnetization the coupling has opposite sign on the two sublattices, resulting
in the additional $\tau_{z}$. The low-energy continuum field theory is subject
to a UV momentum cut-off, $|\mathbf{q}|\leq\Lambda$.

As pointed out in the context of GNY theories [39, 40], the Landau damping of
the bosonic order parameter fluctuations by electronic particle-hole
fluctuations gives rise to a self-energy contribution

$\Pi(\mathbf{q},\omega)=\gamma\sqrt{\mathbf{q}^{2}+\omega^{2}/v_{F}^{2}}$ (2)

to the inverse boson propagator in two spatial dimensions. This non-analytic
self-energy correction arrises from the diagram in Fig. 1(a) from integration
of fermion modes near zero momenta and frequency. It is therefore not
generated within the momentum shell RG, but needs to be included to correctly
capture the universal critical behavior of GNY theories [39, 40]. The form of
the Landau damping does not depend on the number of order parameter components
and is not affected by the fixed-length constraint of the Néel order parameter
field. Note that although the bare Landau damping parameter $\gamma_{0}$ is
determined by the square of the bare Kondo coupling,
$\gamma_{0}\sim\lambda_{0}^{2}$, this relation is not preserved under the RG.
We therefore treat $\gamma$ and $\lambda$ as independent coupling constants.

Under the RG there will be a non-trivial flow of the velocities $c$ and
$v_{F}$. For simplicity, we will focus on the case $v_{F}=c$, which is
preserved under the RG. For convenience, we rescale to dimensionless momenta
$\mathbf{k}=\mathbf{q}/\Lambda$ and frequencies $k_{0}=\omega/(c\Lambda)$ and
absorb the additional prefactors in a redefinition of the coupling constants.
Since both the order parameter and fermion sectors are relativistic and
frequency and momenta enter the zero-temperature field theory in the same way,
the quantum critical behaviour will be described by a dynamical exponent
$z=1$. We will therefore treat frequency and momenta on an equal footing and
impose an isotropic cut-off in 2+1 dimensions,
$\sqrt{\mathbf{k}^{2}+k_{0}^{2}}\leq 1$. Note that the universal critical
behavior is independent of the choice of the UV cut-off scheme.

##  III Landau damped NL$\sigma$M

We start by investigating the effects of Landau damping on the Néel transition
in the case of vanishing Kondo coupling, $\lambda=0$. Our starting point is
the NL$\sigma$M

$S_{N}=\frac{1}{2g}\int_{k}\Omega^{-1}(k)\vec{N}(k)\cdot\vec{N}(-k),$ (3)

where we have defined $k=(\mathbf{k},k_{0})$ and
$\int_{k}=\int\frac{dk_{0}}{2\pi}\int\frac{d^{2}\mathbf{k}}{(2\pi)^{2}}$,
subject to the cut-off $|k|\leq 1$, for brevity, and include the Landau
damping $\gamma$ in the inverse propagator,

$\Omega^{-1}(k)=k^{2}+\gamma|k|=\mathbf{k}^{2}+k_{0}^{2}+\gamma\sqrt{\mathbf{k}^{2}+k_{0}^{2}}.$
(4)

We follow the conventional treatment [37, 38, 41] and decompose
$\vec{N}=(\vec{\pi},\sigma)$ and use the constraint
$\sigma(\mathbf{r},\tau)=\sqrt{1-\vec{\pi}^{2}(\mathbf{r},\tau)}$ to eliminate
$\sigma$ and derive an effective action in terms of the transverse fields
$\vec{\pi}$. In the presence of Landau damping, we need to apply the
constraint in momentum space,
$\sigma(k)=\delta(k)-\frac{1}{2}\int_{q}\vec{\pi}(q)\vec{\pi}(k-q)$. This
results in the effective action

$\displaystyle S_{N}$ $\displaystyle=$
$\displaystyle\frac{1}{2g}\int_{k}\Omega^{-1}(k)\vec{\pi}(k)\cdot\vec{\pi}(-k)-\frac{\rho}{2}\int_{k}\vec{\pi}(k)\cdot\vec{\pi}(-k)$
(5)
$\displaystyle+\frac{1}{16g}\int_{k_{1},\ldots,k_{4}}\delta(k_{1}+k_{2}+k_{3}+k_{4})\left[\Omega^{-1}(k_{1}+k_{2})\right.$
$\displaystyle\left.+\Omega^{-1}(k_{3}+k_{4})\right]\left[\vec{\pi}(k_{1})\cdot\vec{\pi}(k_{2})\right]\left[\vec{\pi}(k_{3})\cdot\vec{\pi}(k_{4})\right],$

where $\rho$ is the density and the corresponding term arises from
exponentiation and expansion of $1/(2\sqrt{1-\vec{\pi}^{2}(\mathbf{r},\tau)})$
from the path-integral measure [41].

Figure 1: One-loop diagrams relevant to our RG calculation. Solid lines
represent fermionic degrees of freedom, wiggly lines the Néel order parameter
fields. (a) The fermionic bubble diagram integrated over small momenta and
frequencies gives rise to the non-analytic Landau damping of long-wavelength
order parameter fluctuations. The momentum-shell contribution of this diagram
contributes to the renormalization of the coupling constant $g$ of the
NL$\sigma$M. (b)-(d) Diagrams relevant for the RG of the Landau-damped
NL$\sigma$M. The diagram (b) is identical to zero, (c) renormalizes the
quadratic gradient terms and hence the coupling constant $g$. The unphysical
mass term generated by (c) is cancelled by the contribution (d) from the
functional integral measure. (e) The fermionic self-energy diagram
renormalizes the overall prefactor of the free fermion action $S_{f}$. The
scaling dimension of the fermion fields is determined such that the prefactor
remains constant. (f) Diagram contributing to the renormalization of the Kondo
coupling $\lambda$.

We integrate out modes with momenta and frequencies from an infinitesimal
shell near the cut-off, $e^{-d\ell}\leq\sqrt{\mathbf{k}^{2}+k_{0}^{2}}\leq 1$,
followed by a rescaling of momenta, $\mathbf{k}\to\mathbf{k}e^{d\ell}$, and
frequencies, $k_{0}\to k_{0}e^{zd\ell}$, with dynamical exponent $z=1$. In
addition, we rescale the transverse spin-fluctuation fields as
$\vec{\pi}(k)\to\vec{\pi}(k)e^{-\Delta_{\pi}d\ell}$.

At one-loop order, the contraction of two order parameter fields,

$\langle\pi_{\alpha}(k)\pi_{\beta}(k^{\prime})\rangle_{0}=g\delta_{\alpha\beta}\delta(k+k^{\prime})\Omega(k),$
(6)

gives rise to the renormalization of the quadratic action by the quartic
vertex. The diagram in Fig. 1(b) vanishes because of $\Omega^{-1}(0)=0$. The
diagram shown in Fig. 1(c) gives rise to a term $\sim
k^{2}\vec{\pi}(k)\cdot\vec{\pi}(-k)$ and hence a renormalization of the
coupling constant $g$. In addition, it produces a mass term
$\sim\vec{\pi}(k)\cdot\vec{\pi}(-k)$ which cancels exactly with the trivial
term from the reduction of the density $\rho$ by the shell contribution [Fig.
1(d)]. Evaluating the momentum shell and frequency integrals and combining
with the rescaling contributions, we obtain the RG equations

$\displaystyle\frac{d}{d\ell}\left(\frac{1}{2g}\right)$ $\displaystyle=$
$\displaystyle(-5+2\Delta_{\pi})\frac{1}{2g}+\frac{1}{(2\pi)^{2}}\frac{1+\gamma/3}{1+\gamma},$
(7) $\displaystyle\frac{d}{d\ell}\left(\frac{\gamma}{2g}\right)$
$\displaystyle=$ $\displaystyle(-4+2\Delta_{\pi})\frac{\gamma}{2g}.$ (8)

The scaling dimension $\Delta_{\pi}$ of the transverse spin-fluctuation fields
needs to be determined such that the constraint of the NL$\sigma$M is
satisfied on all scales. This is only the case if the coupling constant $g$ in
front of the quartic vertex renormalizes in exactly the same way as the $g$ in
front of the quadratic action. Instead of evaluating the second-order, one-
loop diagrams that renormalize the vertex, we employ a trick invented by
Nelson and Pelcovits [41] to include a staggered magnetic field term
$-\frac{h}{2g}\int_{\mathbf{r},\tau}\sigma(\mathbf{r},\tau)$ in the action.
Since the scaling of the magnetic field should not depend on the field
direction, and since the staggered magnetic field couples linearly to the Néel
order parameter field, the scaling dimension of the applied field is equal to
that of the order-parameter field itself,

$\frac{d}{d\ell}\left(\frac{h}{2g}\right)=\Delta_{\pi}\frac{h}{2g}.$ (9)

On the other hand, we can use the constraint to expand
$\sigma(\mathbf{r},\tau)$ in terms of the $\vec{\pi}$ fields and explicitly
compute the one-loop renormalization of the applied field,

$\frac{d}{d\ell}\left(\frac{h}{2g}\right)=(-3+2\Delta_{\pi})\frac{h}{2g}+\frac{1}{(2\pi)^{2}}\frac{h}{1+\gamma}.$
(10)

Equating Eqs. (9) and (10), we obtain

$\Delta_{\pi}=3-\frac{2}{(2\pi)^{2}}\frac{g}{1+\gamma},$ (11)

which results in the coupled RG equations

$\displaystyle\frac{d\tilde{g}}{d\ell}$ $\displaystyle=$
$\displaystyle-\tilde{g}+\frac{1-\gamma/3}{1+\gamma}\tilde{g}^{2},$ (12)
$\displaystyle\frac{d\gamma}{d\ell}$ $\displaystyle=$
$\displaystyle\gamma\left[1-\frac{1+\gamma/3}{1+\gamma}\tilde{g}\right],$ (13)

for $\tilde{g}=\frac{2}{(2\pi)^{2}}g$ and the Landau damping $\gamma$.

For $\gamma=0$ we recover the RG equation of the conventional NL$\sigma$M in
2+1 space-time dimensions. This RG equation exhibits two fixed points: the
attractive, Néel ordered fixed point at $\tilde{g}=0$ and the critical fixed
point at $\tilde{g}=\tilde{g}_{c}=1$. For $\tilde{g}(0)<1$ the RG flow is
towards $\tilde{g}=0$, corresponding to a freezing of transverse spin-
fluctuations on larger and larger scales. On the other hand, for
$\tilde{g}(0)>1$, $\tilde{g}(\ell)\to\infty$, corresponding to a vanishing of
the spin stiffness and indicative of the destruction of long-range order by
spatial and temporal fluctuations.

The coupled RG equations (12) and (13) do not exhibit any additional fixed
points at finite $\gamma$. The RG flow in the $\tilde{g}$-$\gamma$ plane is
shown in Fig. 2. In the antiferromagnetically ordered phase, $\gamma(\ell)$
increases, indicative of damped spin-wave excitations. At the critical fixed
point of the Néel transition the Landau damping $\gamma$ is weakly irrelevant.
The separatrix between the Néel antiferromagnet and the quantum disordered
phase is given by $\tilde{g}\approx
1+\frac{4}{3}\gamma-\frac{4}{9}\gamma^{2}$. Along the separatrix and for an
initial value $\gamma_{0}=\gamma(0)\ll 1$, the Landau damping vanishes as
$\gamma(\ell)=\gamma_{0}/(1+\frac{2}{3}\gamma_{0}\ell)$.

Figure 2: RG flow of the Landau damped NL$\sigma$M as a function of the
rescaled inverse spin stiffness $\tilde{g}=g/g_{c}$ and Landau damping
$\gamma$. The red line separates the Néel AFM from the paramagnet. Along this
separatrix, $\gamma$ renormalizes to zero, demonstrating that the Néel quantum
critical point $P_{c}$ is stable against Landau damping. The increase of
$\gamma$ in the ordered phase indicates that spin-wave excitations are damped.

##  IV Including the Kondo coupling to Dirac fermions

We now include the Kondo coupling $S_{K}$ between the order parameter field to
$N_{f}$ copies of two-component Dirac fermion fields, as given in Eq. (II).
The momentum-shell contribution of the diagram Fig. 1(a) will give rise to an
additional correction to the NL$\sigma$M,

$\displaystyle\delta S_{N}$ $\displaystyle=$
$\displaystyle-\frac{1}{2}\frac{\lambda^{2}}{N_{f}}\int_{q}^{<}\vec{N}(q)\cdot\vec{N}(-q)$
(14)
$\displaystyle\times\int_{k}^{>}\textrm{Tr}\left\\{\tau_{z}G_{\psi}(k)\tau_{z}G_{\psi}(k+q)\right\\},$

where $\int_{q}^{<}$ and $\int_{k}^{>}$ denote frequency-momentum integrals
over $|q|\leq e^{-d\ell}$ and $e^{-d\ell}\leq|k|\leq 1$, respectively. The
fermionic Green function in each of the $N_{f}$ copies is given by

$G_{\psi}(k)=\frac{ik_{0}+k_{x}\tau_{x}+k_{y}\tau_{y}}{k_{0}^{2}+\mathbf{k}^{2}}.$
(15)

Note that the trace in Eq. (14) results in an additional factor of $N_{f}$.
Expanding external momenta/frequencies $q=(\mathbf{q},q_{0})$ to quadratic
order, we obtain

$\delta
S_{N}=-\frac{1}{3}\lambda^{2}\frac{2}{(2\pi)^{2}}d\ell\int_{q}^{<}q^{2}\vec{N}(q)\cdot\vec{N}(-q),$
(16)

resulting in an additional contribution
$d\left(\frac{1}{2g}\right)=-\frac{1}{3}\lambda^{2}\frac{2}{(2\pi)^{2}}d\ell$
to the renormalization of the coupling constant. This changes the RG equations
for $\tilde{g}$ and $\gamma$ to

$\displaystyle\frac{d\tilde{g}}{d\ell}$ $\displaystyle=$
$\displaystyle-\tilde{g}+\left[\frac{1-\gamma/3}{1+\gamma}+\frac{2}{3}\lambda^{2}\right]\tilde{g}^{2},$
(17) $\displaystyle\frac{d\gamma}{d\ell}$ $\displaystyle=$
$\displaystyle\gamma\left[1-\frac{1+\gamma/3}{1+\gamma}\tilde{g}+\frac{1}{3}\lambda^{2}\tilde{g}\right].$
(18)

In order to determine the renormalization of the Kondo coupling constant
$\lambda$, we first need to determine the scaling dimension $\Delta_{\psi}$ of
the fermion fields. The diagram in Fig. 1(e) results in a correction

$\displaystyle\delta S_{f}$ $\displaystyle=$
$\displaystyle-\frac{2\lambda^{2}g}{N_{f}}\int_{k}^{<}\overline{\psi}(k)\left(\int_{q}^{>}D(q)\tau_{z}G_{\psi}(k+q)\tau_{z}\right)\psi(k)$
(19) $\displaystyle=$
$\displaystyle\frac{2}{3N_{f}}\frac{1}{1+\gamma}\lambda^{2}\tilde{g}\,d\ell\int_{k}^{<}\overline{\psi}(k)G_{\psi}^{-1}(k)\psi(k),$

where the factor of two arises from the number of components of the transverse
spin-fluctuation field $\vec{\pi}$, $N_{\pi}=2$.

After rescaling frequency and momenta as before and fermion fields as
$\psi(k)\to\psi(k)e^{-\Delta_{\psi}d\ell}$, we demand that that the prefactor
of $S_{f}$ remains scale invariant, which results in

$\Delta_{\psi}=2-\frac{1}{3N_{f}}\frac{1}{1+\gamma}\lambda^{2}\tilde{g}.$ (20)

The diagram that contributes to the renormalization of the Kondo vertex is
shown in Fig. 1(f) and equals

$\delta
S_{K}=\frac{g\lambda^{3}}{\sqrt{N_{f}}^{3}}\sum_{i}\int_{k_{1},k_{2}}^{<}\pi_{i}(k_{1}-k_{2})\overline{\psi}(k_{1})\Omega_{i}\psi(k_{2}),$
(21)

with coupling matrices

$\displaystyle\Omega_{i}$ $\displaystyle=$
$\displaystyle\sum_{j}\int_{q}^{>}D(q)(\sigma_{j}\otimes\tau_{z})G_{\psi}(q)(\sigma_{i}\otimes\tau_{z})$
(22) $\displaystyle\times G_{\psi}(q)(\sigma_{j}\otimes\tau_{z}).$

Since $G_{\psi}$ is independent of spin, we can evaluate the products of spin
Pauli matrices and carry out the sum over $j$,
$\sum_{j}\sigma_{j}\sigma_{i}\sigma_{j}=(2-N_{\pi})\sigma_{i}$. The momentum-
shell integral is trivial and we indeed find that $\Omega_{i}$ is proportional
to the original Kondo coupling matrix $\sigma_{i}\otimes\tau_{z}$,

$\Omega_{i}=\frac{2}{(2\pi)^{2}}(N_{\pi}-2)\frac{1}{1+\gamma}d\ell(\sigma_{i}\otimes\tau_{z}).$
(23)

However, the result crucially depends upon the number $N_{\pi}$ of order-
parameter components, as discussed in the literature [42, 43]. While the
results for $N_{\pi}=1$ and $N_{\pi}=3$ are equal but of opposite sign, the
diagram vanishes in the relevant case of $N_{\pi}=2$ components, $\delta
S_{Y}=0$.

The rescaling of momenta, frequencies and fields gives rise to the RG equation

$\displaystyle\frac{d\lambda}{d\ell}$ $\displaystyle=$
$\displaystyle(-6+\Delta_{\pi}+2\Delta_{\psi})\lambda$ (24) $\displaystyle=$
$\displaystyle\lambda\left[1-\frac{\tilde{g}}{1+\gamma}-\frac{2}{3N_{f}}\frac{\lambda^{2}\tilde{g}}{1+\gamma}\right]$

for the Kondo coupling $\lambda$.

##  V Analysis of RG flow

We will now discuss the coupled RG equations for the inverse spin stiffness
$\tilde{g}$ (17), the Landau damping $\gamma$ of the Néel order parameter (18)
and the Kondo coupling $\lambda$ to the Dirac fermions (24). In Sec. III we
found that in the absence of Kondo coupling ($\lambda=0$), the Landau damping
$\gamma$ is weakly irrelevant at the Néel quantum critical point $P_{c}$.

Let us first investigate the stability of $P_{c}$ against Kondo coupling in
the absence of Landau damping ($\gamma=0$). In this case the RG equations
reduce to
$d\tilde{g}/d\ell=-\tilde{g}+\left(1+\frac{2}{3}\lambda^{2}\right)\tilde{g}^{2}$
and
$d\lambda/d\ell=\lambda\left(1-\tilde{g}-\frac{2}{3N_{f}}\lambda^{2}\tilde{g}\right)$.
In this case we find a separatrix $\tilde{g}=1-\frac{2}{3}\lambda^{2}$, along
which the flow of the Kondo coupling increases according to
$d\lambda/d\ell=\frac{2}{3}(1-1/N_{f})\lambda^{3}$, resulting in
$\lambda(\ell)=\lambda_{0}/\sqrt{1-\frac{4}{3}\lambda_{0}^{2}(1-1/N_{f})\ell}$.
The Néel quantum critical point is therefore very weakly unstable against the
Kondo coupling to Dirac fermions. The RG flow in the $\tilde{g}$-$\lambda$
plane is shown in Fig. 3.

Figure 3: RG flow as a function of inverse spin stiffness $\tilde{g}=g/g_{c}$
and Kondo coupling $\lambda$ to Dirac fermions with $N_{f}=4$ flavors. At the
Néel quantum critical point $P_{c}$, the Kondo coupling is a weakly relevant
perturbation, indicated by the increase of $\lambda$ along the separatrix
shown in red.

To determine the critical surface in the three-dimensional parameter space of
$\tilde{g}$, $\gamma$ and $\lambda$, we insert a polynomial ansatz
$\tilde{g}=f(\gamma,\lambda)$ into the RG equations (17), (18) and (24). To
second order we obtain

$\tilde{g}=1+\frac{4}{3}\gamma-\frac{4}{9}\gamma^{2}-\frac{2}{3}\lambda^{2}.$
(25)

The critical surface is shown in Fig. 4. As expected, the critical surface
contains the separatrices in the $\lambda=0$ and $\gamma=0$ planes. For
initial values of the coupling constants slightly outside the surface, the RG
flow is away from the surface: the inverse spin stiffness $\tilde{g}$
renormalizes to zero on one side, indicative of a freezing of spin-wave
fluctuations, and to infinity on the other side, corresponding to a quantum
disordered state. The Landau damping $\gamma$ has a stabilizing effect on the
Néel order, while the Kondo coupling $\lambda$ has a destabilizing effect.

Figure 4: RG flow within the critical surface for $N_{f}=4$, relevant to Dirac
electrons on the honeycomb lattice. The Néel quantum critical point $P_{c}$ is
stable against Landau damping $\gamma$ but unstable against Kondo coupling
$\lambda$. For sufficiently strong Landau damping, the RG flow is towards
$P_{c}$ until the trajectories turn to hit the magenta lines, which are given
by $\gamma^{\prime}(\ell)=0$ and $\lambda^{\prime}(\ell)=0$, respectively, and
closely track each other. At this point the RG flow becomes extremely slow and
the parameters acquire small metastable values.

To analyze the competition between $\gamma$ and $\lambda$ within the critical
surface we replace $\tilde{g}$ in the corresponding RG equations, using Eq.
(25),

$\displaystyle\frac{d\gamma}{d\ell}$ $\displaystyle=$
$\displaystyle-\frac{2}{3}\gamma^{2}+\frac{2}{3}\gamma^{3}+\gamma\lambda^{2}$
(26) $\displaystyle\frac{d\lambda}{d\ell}$ $\displaystyle=$
$\displaystyle-\frac{1}{3}\gamma\lambda+\frac{7}{9}\gamma^{2}\lambda+\frac{2}{3}(1-1/N_{f})\lambda^{3},$
(27)

where we have expanded up to cubic order in the coupling constants. For
$N_{f}\geq 4$, the RG equations only exhibit a single fixed point at
$\lambda=0$ and $\gamma=0$, corresponding to the Néel quantum critical point
$P_{c}$. The RG flow in the critical surface and several trajectories obtained
from numerical integration of the RG equations (26) and (27) are shown in Fig.
4 for the case $N_{f}=4$.

The RG flow is best understood in terms of the lines along which (i)
$\gamma^{\prime}(\ell)=0$ and (ii) $\lambda^{\prime}(\ell)=0$, shown in
magenta in Fig. 4, and given by (i)
$\lambda^{2}=\frac{2}{3}(\gamma-\gamma^{2})$ and (ii)
$\lambda^{2}=\frac{1}{2(1-1/N_{f})}(\gamma-\frac{7}{3}\gamma^{2})$,
respectively. These lines merge at $P_{c}$ and because of they exhibit the
same asymptotic functional form, $\lambda\sim\sqrt{\gamma}$, they closely
track each other. As a result, the RG flow becomes very slow in the vicinity
of this pair of lines and it is not possible for trajectories to cross them on
scales relevant to any realistic system size.

The case $N_{f}=4$, relevant to Dirac electrons on the honeycomb lattice, is
the most extreme since in this case the coefficients of the leading
$\sqrt{\gamma}$ terms are identical. For weak Landau damping,
$\lambda^{2}>\frac{2}{3}\gamma$, corresponding to points above the magenta
lines, the flow is towards the regime of strong Kondo coupling. This indicates
that the Néel quantum critical $P_{c}$ point becomes unstable toward Kondo
physics, which falls outside the validity of our analysis.

Figure 5: RG flow of the Landau damping $\gamma$ and the Kondo coupling
$\lambda$ within the critical surface for $N_{f}=2$. The Néel quantum critical
point $P_{c}$ is thermodynamically stable in the blue region. In the purple
region the flow is towards increasing $\lambda$, indicative of strong coupling
Kondo physics. The transition between the two regimes is controlled by a new
multi-critical point $\tilde{P}_{c}$.

On the other hand, if the Landau damping is sufficiently strong,
$\lambda^{2}<\frac{2}{3}\gamma$, both $\lambda(\ell)$ and $\gamma(\ell)$
decrease under the RG. The corresponding trajectories approach $P_{c}$ until
they eventually turn to hit the magenta lines. Here the RG flow practically
comes to a standstill and $\lambda(\ell)$ and $\gamma(\ell)$ reach metastable
plateaux values $\gamma_{*}$ and
$\lambda_{*}^{2}\approx\frac{2}{3}\gamma_{*}$. The non-zero Kondo coupling
leads to the opening of a small electronic gap
$\Delta\sim\lambda_{*}|\langle\vec{N}(\mathbf{r},\tau)\rangle|$ in the Néel
ordered phase where the spin-rotational symmetry is broken. On the critical
surface the finite values $\gamma_{*}$ and $\lambda_{*}$ result in an
anomalous contribution to the scaling dimension $\Delta_{\psi}$ of the fermion
fields, giving rise to non-Fermi liquid behavior. However, the corresponding
critical exponents are non-universal since the behavior is not associated with
a true fixed point.

For $N_{f}<4$ the RG equations (26) and (27) exhibit an additional fixed point
$\tilde{P}_{c}$ at

$\tilde{\gamma}_{c}=\frac{4-N_{f}}{4+3N_{f}},\quad\tilde{\lambda}^{2}_{c}=\frac{8}{3}\frac{N_{f}(4-N_{f})}{(4+3N_{f})^{2}},$
(28)

which merges with the Néel quantum critical point $P_{c}$ as $N_{f}\to 4$,
showing again that the case $N_{f}=4$ is marginal.

In Fig. 5 the RG flow of $\gamma$ and $\lambda$ within the critical surface is
shown for the representative case $N_{f}=2$. In the blue region the RG flow is
towards $P_{c}$, demonstrating that the Néel fixed point is thermodynamically
stable rather than metastable. In the regime of small Kondo coupling
$\lambda$, this stability is achieved by finite but very small Landau damping
$\gamma$. In the purple region the RG flow is towards large values of $\gamma$
and $\lambda$, beyond the validity of our RG equations. The transition between
this Kondo run-away regime and the Néel critical region is described by the
multi-critical fixed point $\tilde{P}_{c}$.

We proceed to analyze the universal critical behavior of the multi-critical
point $\tilde{P}_{c}$ for general $N_{f}<4$. The correlation length exponent
$\tilde{\nu}$ can be obtained from linearizing the RG equation for the inverse
spin stiffness (17) around the critical value $\tilde{g}_{c}\approx
1+\frac{4}{3}\tilde{\gamma}_{c}-\frac{2}{3}\tilde{\lambda}_{c}^{2}$. The
resulting RG equation is of the general form
$d(\tilde{g}-\tilde{g_{c}})/d\ell=\tilde{\nu}^{-1}(\tilde{g}-\tilde{g_{c}})$.
A short calculation gives $\tilde{\nu}=1$ which is identical to the
correlation-length exponent $\nu=1$ at the Néel quantum critical point
$P_{c}$.

From the scaling dimension $\Delta_{\pi}$ (11) of the transverse spin
fluctuations fields $\vec{\pi}$ we obtain the anomalous dimension
$\eta_{\pi}=1$ at the Néel critical point $P_{c}$ and

$\tilde{\eta}_{\pi}=\frac{\tilde{g}_{c}}{1+\tilde{\gamma}_{c}}\approx
1+\frac{1}{9}\frac{(4-N_{f})(12-7N_{f})}{(4+3N_{f})^{2}}$ (29)

at the new multi-critical point $\tilde{P}_{c}$. The additional contribution
to $\tilde{\eta}_{\pi}$ results in a slightly different exponent of the
algebraic order parameter correlations at criticality,
$\langle\vec{\pi}(\mathbf{r})\vec{\pi}(\mathbf{0})\rangle\sim
r^{-D+2-\tilde{\eta}_{\pi}}$, and corrections to other critical exponents,
which can be obtained from the conventional scaling and hyper-scaling
relations.

Due to the finite value $\tilde{\lambda}_{c}$ (28) of the Kondo coupling at
$\tilde{P}_{c}$, the symmetry breaking transition will be accompanied with the
opening of a gap [20]

$\Delta\sim(\tilde{g}_{c}-\tilde{g})^{z\tilde{\nu}}=(\tilde{g}_{c}-\tilde{g})$
(30)

in the Dirac fermion spectrum for $\tilde{g}<\tilde{g}_{c}$, in the Néel
ordered phase. Moreover, at $\tilde{P}_{c}$ the fermions acquire a small
anomalous dimension [see Eq. (20)],

$\tilde{\eta}_{\psi}=\frac{1}{3N_{f}}\frac{\tilde{\lambda}_{c}^{2}\tilde{g}_{c}}{1+\tilde{\gamma}_{c}}\approx\frac{8}{9}\frac{4-N_{f}}{(4+3N_{f})^{2}},$
(31)

which implies that the fermion Green’s function has branch cuts rather than
quasiparticle poles. The multicritical point $\tilde{P}_{c}$ is therefore
associated with non-Fermi liquid behavior. From a scaling analysis of the
fermionic spectral function [20] we find that the quasiparticle pole strength
vanishes as

$Z\sim(\tilde{g}-\tilde{g}_{c})^{(z-1+\tilde{\eta}_{\psi})\tilde{\nu}}=(\tilde{g}-\tilde{g}_{c})^{\tilde{\eta}_{\psi}}$
(32)

as the critical point is approached from the semi-metallic, non-magnetic phase
($\tilde{g}>\tilde{g}_{c}$).

##  VI Comparison with the $\epsilon$-expansion

We now address the question whether the same qualitative behavior can be found
within an $\epsilon$-expansion above the lower critical dimension,
$D=2+\epsilon$. Although such an expansion gives analytic control of the
criticality of the NL$\sigma$M [41, 44], the Padé-Borel extrapolation to
$\epsilon=1$ may be problematic due to a lack of sign oscillations in the
coefficients of the $\epsilon$-expansion of $1/\nu$ [45].

An additional problem arrises when the NL$\sigma$M is coupled to Dirac
fermions since the form of the resulting Landau damping of the Néel order
parameter field explicitly depends on the dimension $D$,
$\Pi(\mathbf{k})=\gamma|\mathbf{k}|^{D-2}$ [40]. At one-loop order the RG
equations for the inverse spin stiffness $\tilde{g}=g/(2\pi)$ and the Landau
damping $\gamma$ in $D=2+\epsilon$ are given by

$\displaystyle\frac{d\tilde{g}}{d\ell}$ $\displaystyle=$
$\displaystyle-\epsilon\tilde{g}+\frac{1-\gamma\epsilon^{2}/4}{1+\gamma}\tilde{g}^{2},$
(33) $\displaystyle\frac{d\gamma}{d\ell}$ $\displaystyle=$
$\displaystyle\gamma\left[2-\epsilon-\frac{1+\gamma\epsilon^{2}/4}{1+\gamma}\tilde{g}\right],$
(34)

where we have determined the scaling dimension of the order parameter field
$\Delta_{\pi}=2+\epsilon-\tilde{g}/(1+\gamma)$ from the renormalization of an
auxiliary magnetic field, as before.

Without Landau damping, $\gamma=0$, we obtain the Néel quantum critical point
at $\tilde{g}_{c}=\epsilon$. At the critical spin stiffness the linearized RG
equation for $\gamma$ is equal to $d\gamma/d\ell=2(1-\epsilon)\gamma$, showing
that near the lower critical dimension the Landau damping is a relevant
perturbation.

Interestingly, the shell contribution of the diagram in Fig. 1(a) is equal to
zero in $D=2$ due to a vanishing angular integral. As a result, the Kondo
coupling $\lambda$ does not contribute to the renormalization of $\tilde{g}$
and $\gamma$, unlike in $D=3$.

For similar reasons, the angular integration over the $D=2$ dimensional shell
causes the fermionic self-energy diagram, shown in Fig. 1(e), to vanish. The
fermion field does therefore not acquire an anomalous dimension and the
scaling dimension is trivial, $\Delta_{\psi}=(3+\epsilon)/2$. From the scaling
dimensions $\Delta_{\pi}$ and $\Delta_{\psi}$ we obtain the renormaization of
the Kondo coupling,

$\frac{d\lambda}{d\ell}=\lambda\left[1-\frac{\tilde{g}}{1+\gamma}\right].$
(35)

At $\tilde{g}_{c}=\epsilon$ and $\gamma=0$ the RG equation reduces to
$d\lambda/d\ell=(1-\epsilon)\lambda$, demonstrating that the Kondo coupling is
a relevant perturbation at the Néel quantum critical point for $\epsilon<1$.
Note that for $\epsilon=1$ both the the Landau damping $\gamma$ and the Kondo
coupling $\lambda$ become marginal, consistent with our calculation in $D=3$.

##  VII Discussion

We have investigated the stability of the Néel quantum critical point of a
two-dimensional quantum antiferromagnet with a Kondo coupling to $N_{f}$
flavors of two-component Dirac fermion fields. For $N_{f}=4$ this would
describe Dirac electrons on the honeycomb lattice with two-fold spin and
valley degeneracies.

The resulting long-wavelength field theory is given by a NL$\sigma$M coupled
to the Dirac fermion fields. It is crucial to account for the Landau damping
of the Néel order parameter field. From simple scaling arguments, the
resulting self-energy correction to the order-parameter propagator is expected
to dominate the IR physics.

At first glance the field theory seems very similar to the Heisenberg GNY
theory, which describes the criticality in a purely electronic model with
strong local repulsions between the Dirac electrons. There are crucial
differences, however. While in the GNY theory the quantum phase transition is
tuned by the mass of the order parameter field, the NL$\sigma$M only contains
gradient terms and the criticality occurs as a function of the inverse spin-
stiffness. It is therefore essential to follow the scale dependence of the
order parameter propagator with both the quadratic gradient terms and the non-
analytic self energy correction from Landau damping. In the Heisenberg GNY
theory on the other hand, the quadratic gradient terms can be discarded.

Another important difference is that the scaling dimension of the transverse
spin-fluctuation field of the NL$\sigma$M is fixed by the requirement that the
constraint $\vec{N}^{2}=1$ is satisfied on all length scales. As a result, the
boson scaling dimension cannot be used to enforce scale invariance of the
Kondo coupling, in contrast to the large-$N_{f}$ Heisenberg GNY theory.

We have employed momentum-shell RG to analyze the scale dependence of the
inverse spin stiffness $g$, the Landau damping $\gamma$ and the Kondo coupling
$\lambda$. Although $\gamma$ and $\lambda$ are initially linked to each other
via the fermionic polarization diagram, the two parameters flow independently
under the RG. At the Néel quantum critical point the scaling dimensions of
both $\gamma$ and $\lambda$ vanish and a bifurcation analysis is required. We
have investigated the coupled RG flow of the two perturbations within the
critical surface $g=f(\gamma,\lambda)$ which contains the unperturbed Néel
quantum critical point and separates the regions where transverse spin
fluctuations freeze or diverge, respectively.

The flow within the critical surface shows that while the Landau damping
$\gamma$ is weakly irrelevant at the Néel critical point, the Kondo coupling
$\lambda$ is a weakly relevant perturbation. Interestingly, the interplay
between the two parameters crucially depends on the number $N_{f}$ of Dirac
fermion flavors. For $N_{f}\geq 4$, sufficiently strong Landau damping renders
the Néel quantum critical point metastable. This is evident from an RG flow
towards the Néel critical point up to scales larger than those relevant to
experiments. This behavior is most pronounced for the marginal case $N_{f}=4$,
representing Dirac electrons on the honeycomb lattice.

For $N_{f}<4$ the Néel critical point becomes thermodynamically stable over a
region where the Landau damping dominates over the Kondo coupling. We have
established a new multicritical point on the critical surface which controls
the transition between the Néel-critical and Kondo-runaway regimes. The finite
values of $\gamma$ and $\lambda$ lead to distinct critical exponents and an
anomalous dimension of the fermion fields, resulting in non-Fermi-liquid
behavior. It would be interesting to investigate if our results are robust at
higher-loop order and if the change of behavior still occurs at $N_{f}=4$ or
at a different number of Dirac fermion flavors.

Finally, we have investigated the problem in $D=2+\epsilon$ space-time
dimensions where the RG calculation of the NL$\sigma$M is controlled. We found
that near the lower critical dimension both the Landau damping and the Kondo
coupling are relevant perturbations at the Néel quantum critical point,
resulting in a run-away flow. Unfortunately, the rich behaviour in 2+1
dimensions is not accessible within an $\epsilon$ expansion above the lower-
critical dimension, and the large-$N_{f}$ limit does not provide analytic
control as in the Heisenberg GNY theory. Nevertheless, our results point
towards interesting novel critical behavior that could potentially be further
investigated using quantum Monte Carlo.

## References

  * Hasan and Kane [2010] M. Z. Hasan and C. L. Kane, Rev. Mod. Phys. 82, 3045 (2010).
  * Qi and Zhang [2011] X.-L. Qi and S.-C. Zhang, Rev. Mod. Phys. 83, 1057 (2011).
  * Vafek and Vishwanath [2014] O. Vafek and A. Vishwanath, Annual Review of Condensed Matter Physics 5, 83 (2014).
  * Yan and Felser [2017] B. Yan and C. Felser, Annual Review of Condensed Matter Physics 8, 337 (2017).
  * Armitage  _et al._ [2018] N. P. Armitage, E. J. Mele, and A. Vishwanath, Rev. Mod. Phys. 90, 015001 (2018).
  * Hasan  _et al._ [2021] M. Z. Hasan, G. Chang, I. Belopolski, G. Bian, S.-Y. Xu, and J.-X. Yin, Nature Reviews Materials 6, 784 (2021).
  * Castro Neto  _et al._ [2009] A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, and A. K. Geim, Rev. Mod. Phys. 81, 109 (2009).
  * Grushin  _et al._ [2013] A. G. Grushin, E. V. Castro, A. Cortijo, F. de Juan, M. A. H. Vozmediano, and B. Valenzuela, Phys. Rev. B 87, 085136 (2013).
  * García-Martínez  _et al._ [2013] N. A. García-Martínez, A. G. Grushin, T. Neupert, B. Valenzuela, and E. V. Castro, Phys. Rev. B 88, 245123 (2013).
  * Daghofer and Hohenadler [2014] M. Daghofer and M. Hohenadler, Phys. Rev. B 89, 035103 (2014).
  * Capponi and Läuchli [2015] S. Capponi and A. M. Läuchli, Phys. Rev. B 92, 085146 (2015).
  * Motruk  _et al._ [2015] J. Motruk, A. G. Grushin, F. de Juan, and F. Pollmann, Phys. Rev. B 92, 085147 (2015).
  * Scherer  _et al._ [2015] D. D. Scherer, M. M. Scherer, and C. Honerkamp, Phys. Rev. B 92, 155137 (2015).
  * Volpez  _et al._ [2016] Y. Volpez, D. D. Scherer, and M. M. Scherer, Phys. Rev. B 94, 165107 (2016).
  * Kurita  _et al._ [2016] M. Kurita, Y. Yamaji, and M. Imada, Phys. Rev. B 94, 125131 (2016).
  * de la Peña  _et al._ [2017] D. S. de la Peña, J. Lichtenstein, and C. Honerkamp, Phys. Rev. B 95, 085143 (2017).
  * Christou  _et al._ [2018] E. Christou, B. Uchoa, and F. Krüger, Phys. Rev. B 98, 161120 (2018).
  * Li  _et al._ [2017] Z.-X. Li, Y.-F. Jiang, S.-K. Jian, and H. Yao, Nature Communications 8, 314 (2017).
  * Herbut [2006] I. F. Herbut, Phys. Rev. Lett. 97, 146401 (2006).
  * Herbut  _et al._ [2009] I. F. Herbut, V. Juricic, and B. Roy, Phys. Rev. B 79, 085116 (2009).
  * Assaad and Herbut [2013] F. F. Assaad and I. F. Herbut, Phys. Rev. X 3, 031010 (2013).
  * Gross and Neveu [1974] D. J. Gross and A. Neveu, Phys. Rev. D 10, 3235 (1974).
  * Zinn-Justin [1991] J. Zinn-Justin, Nuclear Physics B 367, 105 (1991).
  * Lehtinen  _et al._ [2003] P. O. Lehtinen, A. S. Foster, A. Ayuela, A. Krasheninnikov, K. Nordlund, and R. M. Nieminen, Phys. Rev. Lett. 91, 017202 (2003).
  * Chen  _et al._ [2009] J.-H. Chen, W. G. Cullen, C. Jang, M. S. Fuhrer, and E. D. Williams, Phys. Rev. Lett. 102, 236805 (2009).
  * Chen  _et al._ [2011] J.-H. Chen, L. Li, W. G. Cullen, E. D. Williams, and M. S. Fuhrer, Nature Physics 7, 535 (2011).
  * Hentschel and Guinea [2007] M. Hentschel and F. Guinea, Phys. Rev. B 76, 115407 (2007).
  * Uchoa  _et al._ [2011] B. Uchoa, T. G. Rappoport, and A. H. Castro Neto, Phys. Rev. Lett. 106, 016801 (2011).
  * Burch  _et al._ [2018] K. S. Burch, D. Mandrus, and J.-G. Park, Nature 563, 47 (2018).
  * Joy and Vasudevan [1992] P. A. Joy and S. Vasudevan, Phys. Rev. B 46, 5425 (1992).
  * Chandrasekharan and Vasudevan [1994] N. Chandrasekharan and S. Vasudevan, Journal of Physics: Condensed Matter 6, 4569 (1994).
  * Wildes  _et al._ [1998] A. R. Wildes, B. Roessli, B. Lebech, and K. W. Godfrey, Journal of Physics: Condensed Matter 10, 6417 (1998).
  * Wildes  _et al._ [2006] A. R. Wildes, H. M. Rønnow, B. Roessli, M. J. Harris, and K. W. Godfrey, Phys. Rev. B 74, 094422 (2006).
  * Wildes  _et al._ [2015] A. R. Wildes, V. Simonet, E. Ressouche, G. J. McIntyre, M. Avdeev, E. Suard, S. A. J. Kimber, D. Lançon, G. Pepe, B. Moubaraki, and T. J. Hicks, Phys. Rev. B 92, 224408 (2015).
  * Lançon  _et al._ [2018] D. Lançon, R. A. Ewings, T. Guidi, F. Formisano, and A. R. Wildes, Phys. Rev. B 98, 134414 (2018).
  * Kim  _et al._ [2019] K. Kim, S. Y. Lim, J. Kim, J.-U. Lee, S. Lee, P. Kim, K. Park, S. Son, C.-H. Park, J.-G. Park, and H. Cheong, 2D Materials 6, 041001 (2019).
  * Chakravarty  _et al._ [1988] S. Chakravarty, B. I. Halperin, and D. R. Nelson, Phys. Rev. Lett. 60, 1057 (1988).
  * Chakravarty  _et al._ [1989] S. Chakravarty, B. I. Halperin, and D. R. Nelson, Phys. Rev. B 39, 2344 (1989).
  * Isobe  _et al._ [2016] H. Isobe, B.-J. Yang, A. Chubukov, J. Schmalian, and N. Nagaosa, Phys. Rev. Lett. 116, 076803 (2016).
  * Uryszek  _et al._ [2020] M. D. Uryszek, F. Krüger, and E. Christou, Phys. Rev. Research 2, 043265 (2020).
  * Nelson and Pelcovits [1977] D. R. Nelson and R. A. Pelcovits, Phys. Rev. B 16, 2191 (1977).
  * Sur and Roy [2019] S. Sur and B. Roy, Phys. Rev. Lett. 123, 207601 (2019).
  * Uryszek  _et al._ [2019] M. D. Uryszek, E. Christou, A. Jaefari, F. Krüger, and B. Uchoa, Phys. Rev. B 100, 155101 (2019).
  * Polyakov [1975] A. Polyakov, Physics Letters B 59, 79 (1975).
  * Hikami and Brezin [1978] S. Hikami and E. Brezin, Journal of Physics A: Mathematical and General 11, 1141 (1978).


# Switching time of spin-torque-driven magnetization in biaxial ferromagnets

Ankit Shukla<EMAIL_ADDRESS>Holonyak Micro and Nanotechnology
Laboratory, University of Illinois at Urbana-Champaign, Urbana, IL 61801  Arun
Parthasarathy<EMAIL_ADDRESS>Department of Electrical and Computer
Engineering, New York University, Brooklyn 11201, USA  Shaloo Rakheja
<EMAIL_ADDRESS>Holonyak Micro and Nanotechnology Laboratory, University
of Illinois at Urbana-Champaign, Urbana, IL 61801

###### Abstract

We analytically model the magnetization switching time of a biaxial
ferromagnet driven by an antidamping-like spin torque. The macrospin
magnetization dynamics is mapped to an energy-flow equation, wherein a
rational-function approximation of the elliptic integrals for moderate spin
current and small damping results in a closed-form expression of the switching
time. Randomness in the initial angle of magnetization gives the distribution
function of the switching time. The analytic model conforms to the results
obtained from Monte Carlo simulation for a broad range of material parameters.
Our results can ameliorate design and benchmarking of in-plane spin torque
magnetic memory by obviating expensive numerical computation.

Spin torque, switching time, biaxial anisotropy, constant-energy approximation

††preprint: AIP/123-QED

##  I Introduction

Current-induced spin phenomena, such as spin-transfer torque (STT) and spin-
orbit torque (SOT), allow for electrical manipulation of the magnetic order,
and form the basis of emerging spintronic technologies such has non-volatile
memory, Kent and Worledge (2015) magnonic interconnects, Sadovnikov  _et al._
(2017) and radio-frequency oscillators. Chen  _et al._ (2016) Spin current can
transfer angular momentum to a magnetic layer and reorient its magnetization,
similar to how electric current can transfer charge to a capacitor and
modulate its voltage. As shown in Fig. 1(a), an electric current flowing
orthogonal to the plane of the spin valve becomes spin polarized in a
direction parallel to the magnetization of the fixed layer. This spin-
polarized current affects the magnetization of the free layer due to STT.
Slonczewski (1996); Berger (1996) On the other hand, an in-plane electric
current flowing through a nonmagnetic material with spin-orbit coupling is
spin polarized in the plane of the nonmagnetic material, but transverse to the
electric current due to the spin-Hall effect. Sinova  _et al._ (2015) This in-
plane polarized spin current can exert SOT Liu  _et al._ (2012a) on the free-
layer magnetization of the spin valve as shown in Fig. 1(b).

Thin-film magnets—the path toward miniaturized spintronics—are subject to
epitaxial strain from substrate and finite-size effects, which can elicit an
in-plane or a perpendicular spin orientation. Sander (2004) The symmetries in
the energy landscape of thin films up to quadratic (lowest-order) terms of
magnetization components are characterized by a biaxial anisotropy, consisting
of an axis of minimum energy that is ‘easy’ for spins to orient along and an
orthogonal axis of maximum energy which is ‘hard’. The uniaxial anisotropy is
a special case of the biaxial anisotropy where the hard axis is absent. For
perpendicularly magnetized films, a perpendicular easy axis approximates the
anisotropy by assuming symmetry in the plane. For in-plane magnetized films,
an in-plane easy axis and a perpendicular hard axis offer the correct
description.

In STT memory, the perpendicular free-layer configuration is superior to the
in-plane configuration due to its lower switching current, faster speeds, and
higher density. Kent and Worledge (2015) A key problem in writing STT memory
is its vulnerability to dielectric breakdown of the tunnel barrier. This is
addressed in SOT memory, where the writing occurs with an in-plane current
that need not traverse the tunnel barrier. The three-terminal SOT memory
separates the read and write paths, improving memory endurance at the cost of
cell size. However, deterministic switching of the perpendicular free layer in
SOT memory requires either a biasing magnetic field Liu  _et al._ (2012b) or
additional layers in the device stack adding to its fabrication complexity. Oh
_et al._ (2016) The in-plane free-layer configuration of the SOT memory is
preferred due to its fabrication simplicity, magnetic-field-free switching,
and lower switching currents, Fukami  _et al._ (2016) although its writing
speed is inferior to that of the perpendicular SOT memory.

When the spin polarization of the injected spin current is antiparallel to the
stable orientation of the free-layer magnetization, the spin torque is
antidamping-like—it competes with the intrinsic damping to raise the macrospin
energy—until halfway in the magnetization reversal process when it becomes
damping-like—it contributes to damping to cause dissipation of macrospin
energy. The switching process is characterized by the time to reverse the
orientation of the free layer as a function of the input spin current. A
closed-form expression of the switching time is useful to design and optimize
the performance of spin-torque memory. Previous works(Sun, 2000; Liu  _et al._
, 2012a; Pinna, Kent, and Stein, 2013a; d’Aquino, Scalera, and Serpico, 2019)
have derived expressions of the switching time for the uniaxial anisotropy,
but not for the more general biaxial anisotropy presented in this work.

There are three equivalent approaches to analyze the magnetization dynamics:
perturbative approach, (Sun, 2000) constructing a Fokker-Planck
representation, (Apalkov and Visscher, 2005; Taniguchi  _et al._ , 2013) or
using a constant-energy stochastic equation. (Mayergoyz, Bertotti, and
Serpico, 2009; Pinna, Kent, and Stein, 2013b; Newhall and Vanden-Eijnden,
2013; Pinna, Stein, and Kent, 2014) In this work, we adopt the constant-energy
orbit averaging (CEOA) to study magnetization reversal in a biaxial magnetic
system because this approach simplifies a coupled three-dimensional (3D) (two-
dimensional if the magnetization has a constant magnitude) stochastic problem
into a tractable one-dimensional (1D) problem at low temperature and for low-
to-moderate applied spin current. (Pinna, Kent, and Stein, 2013b; Newhall and
Vanden-Eijnden, 2013; Pinna, Stein, and Kent, 2014)

The switching dynamics is modeled as a slow perturbation of the rapid
constant-energy gyration around the easy axis (Sec. II). In the deterministic
limit, a closed-form expression of the switching time as a function of input
spin current, initial magnetization energy, and material parameters is
obtained. Average switching time and the probability distribution of the
switching time follow for an initial Boltzmann distributed ensemble of spins
(Secs. III and IV).

Figure 1: Injection of spin current $I_{s}$ (dashed red line) into a free
magnetic layer produces (a) STT in a spin valve and (b) SOT in a spin-Hall
structure. $\mathbf{m}$ represents the direction of free-layer magnetization,
$\mathbf{\hat{n}}_{p}$ represents the direction of spin polarization, and
$\mathbf{\hat{n}}_{e}$, $\mathbf{\hat{n}}_{i}$, and $\mathbf{\hat{n}}_{h}$ are
the unit vectors along the easy, intermediate, and hard anisotropy axes,
respectively. The flow of electron in the nonmagnetic (NM) layer is opposite
to that of the injected charge current $I$.

##  II Theory

The dynamics of magnetization subject to an effective magnetic field,
intrinsic damping, and spin torque is described by the Landau-Lifshitz-Gilbert
(LLG) equation. Using dimensionless form of physical parameters listed in
Table 1, the LLG equation is (Pinna, Kent, and Stein, 2013b; Yu, Muratov, and
Moore, 2019)

$\begin{split}\frac{\partial{\mathbf{m}}}{\partial\tau}&=-\left(\mathbf{m}\times\mathbf{h}_{\mathrm{eff}}\right)-\alpha\mathbf{m}\times\left(\mathbf{m}\times\mathbf{h}_{\mathrm{eff}}\right)\\\
&-I_{s}\mathbf{m}\times\left(\mathbf{m}\times\mathbf{\hat{n}}_{p}\right)+\alpha
I_{s}\left(\mathbf{m}\times\mathbf{\hat{n}}_{p}\right),\end{split}$ (1)

where $\mathbf{m}$ is the normalized free-layer magnetization, $\alpha$ is the
Gilbert damping constant, $I_{s}$ is the input spin current, and
$\mathbf{\hat{n}}_{p}$ is the unit vector along spin polarization.

Table 1:  Symbols and Definitions in Eq. (1). $M_{s}$, $H_{k}$, and $V$ are the saturation magnetization, Stoner-Wohlfarth field, and volume, respectively, of the free layer. $\mu_{0}$ is the free-space permeability, $\gamma$ is the gyromagnetic ratio, and $\alpha$ is the Gilbert damping constant. Physical parameters | Dimensionless form
---|---
Magnetization $\mathbf{M}$ | $\mathbf{m}=\mathbf{M}/M_{s}$
Time $t$ | $\tau=\gamma\mu_{0}H_{k}t/(1+\alpha^{2})$
Effective magnetic field $\mathbf{H}_{\mathrm{eff}}$ | $\mathbf{h}_{\mathrm{eff}}=\mathbf{H}_{\mathrm{eff}}/H_{k}$
Spin current $\beta_{s}$ | $I_{s}=\beta_{s}/(\mu_{0}M_{s}H_{k}V)$

The effective magnetic field, $\mathbf{h}_{\mathrm{eff}}$, includes
contributions from an internal field produced by the magnetic anisotropy,
externally applied magnetic fields ($\mathbf{H}_{a}$), and a thermal field.
The sum of the internal and external magnetic fields, normalized to the
Stoner-Wohlfarth field, $H_{k}$, is given as

$\mathbf{h}=-\frac{1}{2}\nabla_{\mathbf{m}}\left[g_{L}(\mathbf{m},\mathbf{H}_{a})\right],$
(2)

where $g_{L}$ is the free energy of the macrospin normalized to its uniaxial
energy $K_{u}V$, where $K_{u}=\mu_{0}M_{s}H_{k}/2$ ($\mu_{0}$ is the free-
space permeability and $M_{s}$ is the saturation magnetization of the free
layer). Neglecting higher order anisotropy terms,

$\begin{split}g_{L}(\mathbf{m},\mathbf{H}_{a})&=D_{e}\left(\mathbf{m}\cdot\mathbf{\hat{n}}_{e}\right)^{2}+D_{h}\left(\mathbf{m}\cdot\mathbf{\hat{n}}_{h}\right)^{2}-\frac{2}{H_{k}}\mathbf{m}\cdot\mathbf{H}_{a},\end{split}$
(3)

where $D_{e}$ and $D_{h}$ are the effective anisotropy coefficients along the
easy $(\mathbf{\hat{n}}_{e})$ and hard axes $(\mathbf{\hat{n}}_{h})$,
respectively. For thin-film magnets, the energy landscape is characterized by
a biaxial anisotropy with $D_{e}=-1$ and $D_{h}=M_{s}/H_{k}$. Without loss of
generality, $\mathbf{\hat{n}}_{e}$ and $\mathbf{\hat{n}}_{h}$ are assumed to
coincide with $\mathbf{\hat{x}}$ and $\mathbf{\hat{z}}$ axes, respectively.
Using Eq. (3) in Eq. (2), we obtain
$\mathbf{h}=m_{x}\mathbf{\hat{x}}-Rm_{z}\mathbf{\hat{z}}+\mathbf{H}_{a}/H_{k}$
where $R=D_{h}=M_{s}/H_{k}$.

The thermal field is a Langevin field (Kani, 2017) that is spatially isotropic
and uncorrelated in space and time,

$\displaystyle\begin{split}\langle\mathbf{h}_{T}(t)\rangle&=0,\\\
\langle\mathbf{h}_{T,p}(t_{1})\mathbf{h}_{T,r}(t_{2})\rangle&=D\delta_{p,r}\delta(t_{1}-t_{2}),\end{split}$
(4)

where $p,r$ represent Cartesian coordinates, $\delta_{p,r}$ is the Kronecker
delta, and $\delta(t)$ is the Dirac delta. Assuming the macrospin to be in
thermal equilibrium with the thermal bath and neglecting Joule heating, the
diffusion coefficient

$D=\frac{\alpha
k_{B}T}{\left(1+\alpha^{2}\right)K_{u}V}=\frac{\alpha}{\left(1+\alpha^{2}\right)\Delta_{0}},$
(5)

where $k_{B}$ is the Boltzmann constant, $T$ is the temperature of the bath,
and $\Delta_{0}$ measures the thermal stability of the macrospin (Kani, 2017)
and is referred to as the barrier height of a uniaxial anisotropy magnet.

###  II.1 Constant-Energy Orbit Averaging (CEOA)

The energy of a macrospin is conserved when damping, thermal field (noise),
external magnetic field, and spin torque are absent. In this case,
$\mathbf{m}$ precesses around the easy axis on the unit magnetization sphere
with a fixed macrospin energy $g_{L}$ $(<0)$. Trajectories of conserved
motion, illustrated in Fig. 2(a), are obtained by solving Eqs. (1) and (3)
with $\alpha$, $\mathbf{H}_{a}$, and $I_{s}$ set equal to zero, and
$\norm{\mathbf{m}}=1$. However, with finite damping, $\mathbf{m}$ loses
energy, eventually relaxing to a stable equilibrium state ($m_{x}=\pm 1$). A
non-zero spin torque can pump energy to the macrospin and act against its
inherent damping, causing the magnetization to deviate from its equilibrium
position.

In the case of a small to moderate input spin torque, small damping, and low
temperature, two distinct time-scales of magnetization dynamics emerge: (i) a
fast time-scale associated with constant-energy gyration around the easy axis
and (ii) a slow time-scale corresponding to perpendicular diffusion of
magnetization from one constant-energy orbit to another as a result of
damping, spin torque, and themal field. Figure 2(b) shows one such trajectory
where the magnetization switches from an anti-parallel well to a parallel well
under the simultaneous effects of damping, input spin current, and thermal
field

Figure 2: (a) Constant-energy curves for $g_{L}<0$ in both anti-parallel and
parallel wells for $R=15$. (b) Switching trajectory of a magnetization
initially in the anti-parallel well for $R=15$ and a large input spin current.
External field is assumed to be absent. The arrows shown in the figure mark
the initial (magenta), intermediate (yellow) and final (green) positions of
the magnetization, in the counter-clockwise sense. The dashed red curves in
both figures correspond to a zero energy separatrix.

The rate of change of macrospin energy due to the non-conservative torques is
given as

$\frac{\partial
g_{L}}{\partial\tau}=\nabla_{\mathbf{m}}g_{L}\cdot\frac{\partial\mathbf{{m}}}{\partial\tau}=2\left(Rm_{z}\frac{\partial
m_{z}}{\partial\tau}-m_{x}\frac{\partial m_{x}}{\partial\tau}\right).$ (6)

Averaging the above equation over one time period of the undamped motion
reduces the coupled stochastic dynamics of Eq. (1) to a 1D stochastic dynamics
as (Pinna, Kent, and Stein, 2013b; Pinna, Stein, and Kent, 2014) (see Appendix
A)

$\displaystyle\begin{split}&\left\langle\frac{\partial
g_{L}}{\partial\tau}\right\rangle=\frac{\pi\alpha}{K\left(R,g_{L}\right)}\sqrt{\frac{R-g_{L}}{1+R}}\left[\frac{I_{s}}{\alpha}\left(1+g_{L}\right)\right.\\\
&-\left.\frac{2}{\pi}\sqrt{\left(1+R\right)\left(R-g_{L}\right)}\left\\{E\left(R,g_{L}\right)+g_{L}K\left(R,g_{L}\right)\right\\}\right]\\\
&+2\sqrt{\frac{\alpha}{\Delta_{0}}}\sqrt{\frac{\left(R-g_{L}\right)}{K\left(R,g_{L}\right)}}\sqrt{E\left(R,g_{L}\right)+g_{L}K\left(R,g_{L}\right)}\circ\dot{W}_{g_{L}}.\end{split}$
(7)

Here, $\circ$ denotes multiplication of thermal noise in the Stratonovich
sense (Pinna, Kent, and Stein, 2013b; Kani, 2017), while
$K\left(R,g_{L}\right)$ and $E\left(R,g_{L}\right)$ are the complete elliptic
integrals of the first and second kind, respectively. Time averaging
$\frac{\partial g_{L}}{\partial\tau}$ over a period of precessional motion
enables us to study the dynamics due to slow diffusion of energy with respect
to fast periodic oscillations (Pinna, Kent, and Stein, 2013b; Pinna, Stein,
and Kent, 2014).

Figure 3: (a) Threshold spin current versus $R$. The asterisk denotes
$R=R_{c}=5.094$ where $I_{s}^{\mathrm{th0}}=I_{s}^{\mathrm{th1}}$. For
$I_{s}>\mathrm{max}\left(I_{s}^{\mathrm{th0}},I_{s}^{\mathrm{th1}}\right)$,
deterministic switching (DS) occurs. Thermal switching region is denoted as
TS, while limit cycles are denoted as LC. (b) For $R=3.0$ and
$I_{s}>I_{s}^{\mathrm{th1}}$, the rate of change of macrospin energy is
positive for $g_{L}\in(-1,0)$. For
$I_{s}^{\mathrm{th0}}<I_{s}<I_{s}^{\mathrm{th1}}$, thermal assistance is
necessary to cause magnetization reversal. For $I_{s}<I_{s}^{\mathrm{th0}}$,
switching occurs due to thermal noise. (c) For $R=15.0$ and
$I_{s}^{\mathrm{th0}}<I_{s}$, reversal is deterministic as the rate of energy
change is always positive. For $I_{s}<I_{s}^{\mathrm{th1}}$, thermal switching
dominates, while in the intermediate region, thermal assistance is required to
cause switching.

Assuming deterministic dynamics and using Eq. (7) to evaluate for zero energy
flow at $g_{L}=0$ and $g_{L}=-1$ leads to two different threshold currents
(Sun, 2000; Pinna, Kent, and Stein, 2013b; Taniguchi  _et al._ , 2013): the
minimum current required to push the magnetization over the energy barrier
into the adjoining basin,
$I_{s}^{\mathrm{th0}}=\alpha\left[\frac{2}{\pi}\sqrt{R\left(1+R\right)}\right]$
and the minimum current required to move the magnetization away from stable
equilibrium, $I_{s}^{\mathrm{th1}}=\alpha[R/2+1]$. $R_{c}$ denotes the
critical value of $R$ for which $I_{s}^{\mathrm{th0}}$ equals
$I_{s}^{\mathrm{th1}}$. The threshold currents demarcate regions of
deterministic switching from those that require thermal assistance as shown in
Fig. 3(a). Figures 3(b) and 3(c) show the rate of change of energy for
deterministic dynamics. It is positive for the complete range of macrospin
energy only for
$I_{s}>\mathrm{max}\left(I_{s}^{\mathrm{th0}},I_{s}^{\mathrm{th1}}\right)$
which is consistent with Fig. 3(a). For other values of current, thermal
assistance is required for switching. The CEOA is valid when the variation in
macrospin energy over one precessional cycle is small, i.e.
$|T\langle\frac{\partial g_{L}}{\partial\tau}\rangle|\ll\max[|g_{L}|]=1$. To
satisfy this constraint on the variation of macrospin energy, the maximum spin
current is given as
$I_{s}^{\mathrm{thM}}=I_{s}^{\mathrm{th0}}\left[1+1/\left(8\alpha\sqrt{R}\right)\right]$.

##  III Analytic switching time model

The switching time due to spin torque is defined as the time required for the
macrospin energy to change from an initial value $g_{L_{i}}$ to its final
value $g_{L_{f}}$. Analytic solutions of the switching time are obtained by
approximating the elliptic integrals in Eq. (7) as follows

$\frac{2}{\pi}\left[E\left(R,g_{L}\right)+g_{L}K\left(R,g_{L}\right)\right]=A(R)g_{L}^{2}+B(R)g_{L}+C(R),$
(8a) $K(x)=\frac{\pi}{2}\left[\frac{x-4}{2x-4}\right],$ (8b)

where $x=R\frac{\left(1+g_{L}\right)}{\left(R-g_{L}\right)}$. Defining
$\tilde{I}_{s}=I_{s}/\alpha$ and using the approximate representation of the
rate of change of magnetization energy, we obtain the magnetization switching
time (see Appendix B for derivation)

$\displaystyle\begin{split}\tau_{s}&=\frac{1}{4\alpha}\int_{g_{L_{i}}}^{g_{L_{f}}}\left[\frac{3R-g_{L}(R+4)}{R-g_{L}(R+2)}\right]\sqrt{\frac{1+R}{R-g_{L}}}\times\\\
&\frac{\partial
g_{L}}{\left[\tilde{I}_{s}\left(1+g_{L}\right)-\sqrt{1+R}\sqrt{R-g_{L}}\left(Ag_{L}^{2}+Bg_{L}+C\right)\right]}.\end{split}$
(9)

In the above equation, the parameters $A$, $B$, and $C$ are functions of $R$
given as $k_{1}+k_{2}R^{k_{3}}$, where the values of $k_{1}$, $k_{2}$ and
$k_{3}$ are chosen for different intervals of $R$ to reduce the error in
approximating the elliptic integrals. See Table 2 for details. Note that the
values of $A$, $B$, and $C$ are independent of the device geometry and depend
only on $R=M_{s}/H_{k}$. Equation (9) is simplified and integrated using
partial fractions to arrive at the closed-form expression of switching time:

$\displaystyle\begin{split}\tau_{s}&=\frac{1}{2\alpha\left(R+2\right)}\bBigg<EMAIL_ADDRESS>\end{split}$
$\displaystyle\begin{split}&+\left.\sqrt{\frac{R(1+R)}{R+2}}\left\\{\frac{\log{\left[\frac{\sqrt{(R-g_{L_{f}})(R+2)}-\sqrt{R(1+R)}}{\sqrt{(R-g_{L_{i}})(R+2)}-\sqrt{R(1+R)}}\right]}}{A\prod\limits_{n=1}^{5}\left(\frac{R(1+R)}{R+2}-w_{i}\right)}\right.\right.\\\
&+\bBigg@{5}.\left.\frac{\log{\left[\frac{\sqrt{(R-g_{L_{f}})(R+2)}+\sqrt{R(1+R)}}{\sqrt{(R-g_{L_{i}})(R+2)}+\sqrt{R(1+R)}}\right]}}{A\prod\limits_{n=1}^{5}\left(\frac{R(1+R)}{R+2}+w_{i}\right)}\right\\}\bBigg@{5}],\end{split}$
(10)

where $w_{i}$’s are the roots of a fifth-degree polynomial-
$x^{5}-(\frac{B}{A}+2R)x^{3}+\frac{\tilde{I}_{s}}{A\sqrt{1+R}}x^{2}+(R^{2}+\frac{B}{A}R+\frac{C}{A})x-\frac{\tilde{I}_{s}\sqrt{1+R}}{A}$
with $N=(R+4)+\frac{2R(1+R)}{(R+2)w_{i}^{2}-R(1+R)}$, and
$D=5Aw_{i}^{4}-3(B+2AR)w_{i}^{2}+2\frac{\tilde{I}_{s}}{\sqrt{1+R}}w_{i}+(AR^{2}+BR+C)$.

Table 2:  Parameters $A$, $B$, and $C$ in Eq. (9) are approximated as $k_{1}+k_{2}R^{k_{3}}$, where $k_{1}$, $k_{2}$, and $k_{3}$ depend on $R$. $R$ |  | $k_{1}$ | $k_{2}$ | $k_{3}$
---|---|---|---|---
[1, 3] | $A$ | 0.35661 | -0.51244 | -0.38689
$B$ | 1.05148 | -0.55504 | -0.28598
$C$ | 0.61670 | 0.03018 | -1.00153
[3, 50] | $A$ | 0.20223 | -0.38439 | -0.68424
$B$ | 0.81746 | -0.34729 | -0.63939
$C$ | 0.61765 | 0.02994 | -1.08243
[50, 100] | $A$ | 0.17370 | -0.51992 | -0.97986
$B$ | 0.78501 | -0.48295 | -0.97726
$C$ | 0.61755 | 0.02625 | -1.01442

A major advantage of this analytic result is that the approximations of
elliptic integrals are independent of the input spin current and Gilbert
damping. Additionally, the results obtained in this work are valid for a broad
range of $R$ as opposed to prior works (Sun, 2000; Vincent  _et al._ , 2014)
that are valid only for $R<R_{c}$ (= 5.09).

Considering $R\rightarrow\infty$ in Eq. (7), the elliptic integrals are
approximated as (Hollos, 2012)

$\displaystyle\begin{split}K(x)=\frac{\pi}{2}\left[\frac{x-4}{2x-4}\right];E(x)=\frac{\pi}{2}\left[1-\frac{x}{4}\frac{x^{2}-28x+64}{4x^{2}-40x+64}\right],\end{split}$

where $x=R\frac{1+g_{L}}{R-g_{L}}$. These approximations of elliptic integrals
simplify the energy flow equation and the switching time is (see Appendix C
for details)

$\displaystyle\begin{split}\tau_{s}&=\frac{1}{2\alpha\left[\tilde{I}_{s}-0.5R\right]}\left[\log{\left[\frac{1+g_{L_{f}}}{1+g_{L_{i}}}\right]}\right.\\\
&+\left.\frac{b(a-4)(a-8)}{32(a-b)}\log{\left[\frac{1+g_{L_{f}}-a}{1+g_{L_{i}}-a}\right]}\right.\\\
&-\left.\frac{a(b-4)(b-8)}{32(a-b)}\log{\left[\frac{1+g_{L_{f}}-b}{1+g_{L_{i}}-b}\right]}\right],\end{split}$
(11)

where $a$ and $b$ are the roots of the quadratic equation
$x^{2}-\left(\frac{160\tilde{I}_{s}-60R}{16\tilde{I}_{s}-7R}\right)x+\left(\frac{256\tilde{I}_{s}-128R}{16\tilde{I}_{s}-7R}\right)=0$.

In the limit $R\to 0$, the energy landscape becomes uniaxial and the switching
time is (Pinna, Kent, and Stein, 2013a, b; Sun, 2000; d’Aquino, Scalera, and
Serpico, 2019; Liu  _et al._ , 2014) (see Appendix D for details)

$\displaystyle\begin{split}\tau_{s}&=\frac{1}{2\alpha\left(\tilde{I}_{s}^{2}-1\right)}\left[\tilde{I}_{s}\left\\{\log{\left[\frac{1+\sqrt{-g_{L_{i}}}}{1+\sqrt{-g_{L_{f}}}}\right]}\right.\right.\\\
&\left.\left.-\log{\left[\frac{1-\sqrt{-g_{L_{i}}}}{1-\sqrt{-g_{L_{f}}}}\right]}\right\\}-\log{\left[\frac{1+g_{L_{i}}}{1+g_{L_{f}}}\right]}\right.\\\
&+\left.2\log{\left[\frac{\tilde{I}_{s}-\sqrt{-g_{L_{i}}}}{\tilde{I}_{s}-\sqrt{-g_{L_{f}}}}\right]}\right].\end{split}$
(12)

Figure 4: Histogram distribution of $g_{L}$ using $\alpha=0.03$ and
$\Delta_{0}=102.4$ for (a) $R$ = 3.0 and (b) $R$ = 30.0. Inset compares the
cumulative distribution function obtained analytically (Eq. (16)) and
numerically.

###  III.1 Equilibrium Distribution

In the absence of input spin current, the magnetization is considered to be in
thermal equilibrium in its stable energy well. An average switching time,
$\langle\tau_{s}\rangle$, is obtained by averaging $\tau_{s}$ over the
equilibrium energy distribution, which in the case of a large energy barrier
is the Boltzmann distribution given as

$\displaystyle\begin{split}w_{eq}\left(\mathbf{m}\right)&=\frac{1}{Z(\Delta_{0},R)}\exp\left(-\Delta_{0}g_{L}(\mathbf{m},0)\right),\\\
\end{split}$ (13)

where $Z(\Delta_{0},R)$ is the partition function. The Boltzmann distribution
function is evaluated in terms of the random variable $G_{L}=-M_{x}^{2}$ (see
Appendix E for details of the transformation). Accordingly the probability
density function

$\rho_{G_{L}}\left(g_{L}\right)=\frac{1}{Z(\Delta_{0},R)}\frac{\exp\left(-\Delta_{0}g_{L}\right)}{\sqrt{-g_{L}}},$
(14)

the partition function

$\displaystyle\begin{split}Z(\Delta_{0},R)&=\int_{-1}^{0}\rho_{G_{L}}[x]dx=\frac{2\exp\left(\Delta_{0}\right)F[\sqrt{\Delta_{0}}]}{\sqrt{\Delta_{0}}},\end{split}$
(15)

and the cumulative distribution function

$\displaystyle\begin{split}P_{G_{L}}\left(g_{L}\right)&=\int_{-1}^{g_{L}}\rho_{G_{L}}[x]dx\\\
&=1-\exp\left(-\Delta_{0}\left(1+g_{L}\right)\right)\frac{F\left[\sqrt{-\Delta_{0}g_{L}}\right]}{F\left[\sqrt{\Delta_{0}}\right]},\end{split}$
(16)

where $F[x]=\exp(-x^{2})\int_{0}^{x}\exp(y^{2})dy$ is the Dawson’s integral.

###  III.2 Average Switching Time and Model Validation

The magnetization is considered to have switched successfully when it crosses
the separatrix $\left(g_{L_{f}}=0\right)$ and consequently moves into the
adjoining energy well. Once the magnetization moves into the target energy
well, the spin current could be switched off. The magnetization would
eventually settle into its stable well due to its intrinsic damping.
Therefore, $\langle\tau_{s}\rangle$ is given as

$\langle\tau_{s}\rangle=\int_{-1}^{0}\rho_{G_{L}}\left(g_{L_{i}}\right)\
\tau_{s}\left(g_{L_{i}},g_{L_{f}}=0\right)\ dg_{L_{i}}.$ (17)

Figure 5: Average switching time $\langle
t_{s}\rangle=\frac{\left(1+\alpha^{2}\right)}{\gamma\mu_{0}}\frac{\langle\tau_{s}\rangle}{H_{k}}$
as a function of injected spin current density, $J_{s}$, for different values
of $R$. The lower horizontal dashed black line in each figure is the threshold
current demarcating region of deterministic switching from that of thermal
activation. The upper horizontal dashed cyan line is the maximum current for
which CEOA is valid. The numerical data is obtained for an ensemble of
$10^{4}$ macrospins.

To benchmark our analytic results, we solve Eq. (1) numerically using the Heun
integration scheme implemented in CUDA and run in parallel on GPUs. Numerical
simulations were calibrated against published results to ensure their
accuracy. (Pinna, Kent, and Stein, 2013b, a; Kani, Rakheja, and Naeemi, 2016)
For all simulations, the time step of integration was set as 0.3 ps. Other
material parameters are $\alpha=0.03$, $K_{u}=0.6\times 10^{6}$ J/m3, and
$V=\pi\times 15\times 7.5\times 2\times 10^{-27}$ m-3, while the value of
$M_{s}$ was varied. Simulations were also conducted for other parameter values
but are not reported here. Generally, the key features and trends of switching
dynamics remain the same as material parameters are varied. For each set of
parameters, simulations were performed on an equilibrium ensemble of $10^{4}$
independent macrospins. Figure 4 shows that Eqs. (14) and (16) describe the
numerical distribution very well.

Figure 5 shows the average switching time $\langle
t_{s}\rangle=\frac{\left(1+\alpha^{2}\right)}{\gamma\mu_{0}}\frac{\langle\tau_{s}\rangle}{H_{k}}$
as a function of the input spin current density,
$J_{s}=(4eK_{u}V/\hbar)(I_{s}/Ar)$, where $e$ is the electron charge, $\hbar$
is the reduced Planck’s constant and $Ar$ is the cross-section area of the
magnet. For $R=0$ (uniaxial anisotropy), there is an excellent agreement
between numerical data and closed-form solutions given in Eqs. (10) and (12)
for moderate to large current levels. However, for current levels approaching
the threshold switching currents, numerical results predict a lower average
switching time as the presence of thermal noise aids the switching process.

For $R>0$ (biaxial anisotropy), $\langle t_{s}\rangle$ obtained using Eq. (11)
predicts a larger average switching time compared to the numerical results for
current levels comparable to the threshold value. This is expected since the
analytic solutions neglect the effect of thermal noise during the switching
process. As the input spin current density ($J_{s}$) increases beyond the
threshold value, $J_{s}^{\mathrm{th0}}$
($=(4eK_{u}V/\hbar)(I_{s}^{\mathrm{th0}}/Ar)$), the agreement between analytic
and numerical results improves. Near the threshold current level, the average
switching time obtained from Eq. (10) is slightly lower than that obtained
from numerical results. This slight deviation is due to the quadratic
approximation of elliptic integrals in Eq. (7). Note that for
$J_{s}>J_{s}^{\mathrm{thM}}$, the validity of CEOA is not fully justified,
which could be the source of mismatch between analytic and numerical models.
Here, $J_{s}^{\mathrm{thM}}$ is $(4eK_{u}V/\hbar)(I_{s}^{\mathrm{thM}}/Ar)$.

##  IV Switching time distribution and Error Rate

###  IV.1 Probability and Cumulative Distribution Functions

Defining a random variable switching time $T_{s}$, the fraction of macrospins
in an ensemble that have switched from an anti-parallel to a parallel state at
time $\tau_{s}$ is

$\displaystyle\begin{split}&P_{T_{s}}\left[\tau_{s}\right]=1-\int_{-1}^{g_{L_{i}}(\tau_{s})}\frac{1}{Z(\Delta_{0},R)}\frac{\exp\left(-\Delta_{0}x\right)}{\sqrt{-x}}dx\\\
&=\exp\left(-\Delta_{0}\left(1+g_{L_{i}}(\tau_{s})\right)\right)\frac{F\left[\sqrt{-\Delta_{0}g_{L_{i}}(\tau_{s}})\right]}{F\left[\sqrt{\Delta_{0}}\right]},\end{split}$
(18)

where $g_{L_{i}}(\tau_{s})$ corresponds to the initial energy for switching
time $\tau_{s}$. An exact analytic expression for $g_{L_{i}}(\tau_{s})$ is not
feasible; however, the monotonically decreasing nature of Eqs. (10)-(12) makes
it possible to numerically invert (inv, 2016) them.
$P_{T_{s}}\left[\tau_{s}\right]$ includes all macrospins in an ensemble with
$g_{L}\geq g_{L_{i}}(\tau_{s})$ as we consider only deterministic switching in
this work. $P_{T_{s}}\left[\tau_{s}\right]$ also refers to the probability of
switching of a macrospin with $\tau_{s}\geq T_{s}$; therefore,
$P_{T_{s}}\left[\tau_{s}\right]=Pr\left[T_{s}\leq\tau_{s}\right]$. Finally,
the probability density function of the switching time is

$\displaystyle\begin{split}\rho_{T_{s}}(\tau_{s})&=\frac{dP_{T_{s}}\left[\tau_{s}\right]}{d\tau_{s}}\\\
&=\frac{\sqrt{\Delta_{0}}}{2F[\sqrt{\Delta_{0}}]}\frac{\exp\left(-\Delta_{0}(1+g_{L_{i}}(\tau_{s}))\right)}{\sqrt{-g_{L_{i}}(\tau_{s})}}\
\left|\frac{dg_{L_{i}}(\tau_{s})}{d\tau_{s}}\right|.\\\ \end{split}$ (19)

###  IV.2 Write-Error Rate

The write-error rate (WER) quantifies the probability of unsuccessful spin
torque switching of the magnet. Using Eq. (18), the WER is

$\displaystyle\begin{split}\text{WER}&=1-P_{T_{s}}[\tau]\\\
&=1-\exp\left(-\Delta_{0}\left(1+g_{L_{i}}(\tau_{s})\right)\right)\frac{F\left[\sqrt{-\Delta_{0}g_{L_{i}}(\tau_{s})}\right]}{F\left[\sqrt{\Delta_{0}}\right]}.\end{split}$
(20)

The probability distribution function (pdf) and WER obtained using analytic
results of Eqs. (19) and (20) are compared against numerical solution in Figs.
6-9 for an ensemble of $10^{6}$ macrospins. Simulation parameters are noted in
Table 3. It is observed that for $R=0$, the accuracy of the pdf and WER
improves as the applied spin current increases. (Kani, Rakheja, and Naeemi,
2016) For $R\geq 15$, the accuracy of analytic solutions also increases as
current increases from $J_{s}^{\mathrm{th0}}$ toward $J_{s}^{\mathrm{thM}}$.
However, for spin currents larger than $J_{s}^{\mathrm{thM}}$, the accuracy of
analytic results drops again as the validity of CEOA becomes questionable. To
arrive at analytic results reported in Figs. 7-9, we numerically invert Eq.
(11) due to its simplicity.

Table 3:  Area and threshold spin current to switch the magnetization of free layer for various values of $R$. $R$ | Area $\left(\mathrm{nm}^{2}\right)$ | $J_{s}^{\mathrm{th}}\ \mathrm{(A/cm^{2})}$ | $J_{s}^{\mathrm{thM}}\ \mathrm{(A/cm^{2})}$
---|---|---|---
0 | $\pi\times 15\times 15$ | $5.5\times 10^{6}$ | $2.9\times 10^{7}$
15 | $\pi\times 15\times 7.5$ | $2.15\times 10^{8}$ | $4.47\times 10^{8}$
50 | $\pi\times 15\times 7.5$ | $7.02\times 10^{8}$ | $1.12\times 10^{9}$
100 | $\pi\times 15\times 7.5$ | $1.40\times 10^{9}$ | $1.98\times 10^{9}$

##  V Conclusion

Analytic models of average switching time, probability distribution function
of switching times, and the write-error rate developed in this paper for thin-
film magnets with biaxial anisotropy show good agreement against numerical
results for moderate to large spin current densities. In the vicinity of the
threshold spin current density, the error between analytic and numerical data
is significant due to thermal noise. For very large spin current densities,
the constant energy orbit averaging approach adopted in this work becomes
inadequate, even though the error between numerical and analytic results is
well under a tolerance limit. The models of this paper should complement
experimental results and aid the analysis, design and development of non-
volatile memory driven by both spin-transfer and spin-orbit torques.

Figure 6: All results reported for $R$ = 0 (uniaxial anisotropy) and
$J_{s}^{\mathrm{th1}}=5.5\times 10^{6}$ A/cm2. The top panel shows the pdf,
while the bottom panel shows the write-error rate (WER). The accuracy of
analytical results improves as the input spin current density increases with
respect to $J_{s}^{\mathrm{th1}}$. Figure 7: All results are reported for $R$
= 15 and $J_{s}^{\mathrm{th0}}=2.15\times 10^{8}$ A/cm2. The top panel shows
the pdf, while the bottom panel shows the write-error rate (WER). The accuracy
of analytical solutions improves for spin current densities larger than the
threshold current density. However, CEOA is strictly valid for
$J_{s}<J_{s}^{\mathrm{thM}}=2.1J_{s}^{\mathrm{th0}}$. Figure 8: All results
are reported for $R$ = 50 and $J_{s}^{\mathrm{th0}}=7.02\times 10^{8}$ A/cm2.
The top panel shows the pdf, while the bottom panel shows the write-error rate
(WER). The accuracy of analytic results improves for spin current density
larger than the threshold current density but drops again for
$J_{s}>J_{s}^{\mathrm{thM}}=1.6J_{s}^{\mathrm{th0}}$. Figure 9: All results
are reported for $R$ = 100 and $J_{s}^{\mathrm{th0}}=1.4\times 10^{9}$ A/cm2.
The top panel shows the pdf, while the bottom panel shows the write-error rate
(WER). The accuracy of analytical results improves for as $J_{s}$ increases
above $J_{s}^{\mathrm{th0}}$ but decreases when
$J_{s}>J_{s}^{\mathrm{thM}}=1.42J_{s}^{\mathrm{th0}}$.

##  Appendix A Simplifying Eq. (6)

Without any loss of generality, we consider the easy and hard axes to coincide
with the $\mathbf{\hat{x}}$ and $\mathbf{\hat{z}}$ axes, respectively. Next,
we consider spin polarization of the fixed layer to be in the plane of the
magnet at an angle $\phi$ with the easy axis, therefore,
$\mathbf{\hat{n}}_{p}=\cos{\phi}\ \mathbf{\hat{x}}+\sin{\phi}\
\mathbf{\hat{y}}$. Therefore, for zero external magnetic field
($\mathbf{H}_{a}=0$), the effective magnetic field
$\mathbf{h}_{\mathrm{eff}}=m_{x}\mathbf{\hat{x}}-Rm_{z}\mathbf{\hat{z}}+\mathbf{h}_{T}$,
where $\mathbf{h}_{T}$ is the thermal field. As a result the three components
of Eq. (1) become

$\displaystyle\begin{split}\frac{\partial{m}_{x}}{\partial\tau}&=Rm_{y}m_{z}+\alpha
m_{x}\left(1-m_{x}^{2}+Rm_{z}^{2}\right)\\\
&+I_{s}\left(\cos{\phi}-m_{x}\left(\cos{\phi}\ m_{x}+\sin{\phi}\
m_{y}\right)\right)\\\ &-\alpha I_{s}\ m_{z}\sin{\phi}\\\
&+n_{S,x},\end{split}$ (21)

$\displaystyle\begin{split}\frac{\partial{m}_{y}}{\partial\tau}&=-\left(R+1\right)m_{x}m_{z}-\alpha
m_{y}\left(m_{x}^{2}-Rm_{z}^{2}\right)\\\
&+I_{s}\left(\sin{\phi}-m_{y}\left(\cos{\phi}\ m_{x}+\sin{\phi}\
m_{y}\right)\right)\\\ &+\alpha I_{s}\ m_{z}\cos{\phi}\\\
&+n_{S,y},\end{split}$ (22)

$\displaystyle\begin{split}\frac{\partial{m}_{z}}{\partial\tau}&=\ m_{x}m_{y}\
-\ \alpha m_{z}\left(R+m_{x}^{2}\ -\ Rm_{z}^{2}\right)\\\
&-I_{s}m_{z}\left(\cos{\phi}\ m_{x}\ +\ \sin{\phi}\ m_{y}\right)\\\ &+\alpha
I_{s}\left(\sin{\phi}\ m_{x}\ -\ \cos{\phi}\ m_{y}\right)\\\
&+n_{S,z},\end{split}$ (23)

where each of $n_{S,p}$ denote thermal noise component in the Stratonovich
sense. Assuming the thermal field $\mathbf{h}_{T}$ as

$\mathbf{h}_{T}\left(t\right)=\sqrt{\frac{\alpha}{\left(1+\alpha^{2}\right)\Delta_{0}}}\left(\frac{\partial\mathbf{W}}{\partial\tau}\right)=\sqrt{D}\mathbf{\dot{W}},$

we get

$\displaystyle\begin{split}&\begin{pmatrix}n_{S,x}\\\ n_{S,y}\\\
n_{S,z}\end{pmatrix}=\mathcal{D}\circ\begin{pmatrix}h_{T,x}\\\ h_{T,y}\\\
h_{T,z}\end{pmatrix}=\sqrt{D}\ \mathcal{D}\circ\begin{pmatrix}\dot{W}_{x}\\\
\dot{W}_{y}\\\ \dot{W}_{z}\end{pmatrix},\end{split}$ (24)

where $\mathbf{\dot{W}}$ represents a 3D stochastic Wiener process whose each
component is a Gaussian random variable with zero mean and unit standard
deviation. $\mathcal{D}$ is referred to as the diffusion matrix and is given
as

$\begin{pmatrix}\alpha\left(1-m_{x}^{2}\right)&m_{z}-\alpha
m_{x}m_{y}&-m_{y}-\alpha m_{x}m_{z}\\\ -m_{z}-\alpha
m_{x}m_{y}&\alpha\left(1-m_{y}^{2}\right)&m_{x}-\alpha m_{y}m_{z}\\\
m_{y}-\alpha m_{x}m_{z}&-m_{x}-\alpha
m_{y}m_{z}&\alpha\left(1-m_{z}^{2}\right)\end{pmatrix}.$

Substituting Eqs. (21) and (23) into Eq. (6) leads to

$\displaystyle\begin{split}\frac{\partial
g_{L}}{\partial\tau}&=2Rm_{z}\left[m_{x}m_{y}-\alpha
m_{z}\left(R-g_{L}\right)\right.\\\ &-I_{s}m_{z}\left(\cos{\phi}\
m_{x}+\sin{\phi}\ m_{y}\right)\\\ &+\left.\alpha I_{s}\left(\sin{\phi}\
m_{x}-\cos{\phi}\ m_{y}\right)\right]\\\ &-2m_{x}\left[Rm_{y}m_{z}+\alpha
m_{x}\left(1+g_{L}\right)\right.\\\
&+I_{s}\left(\cos{\phi}-m_{x}\left(\cos{\phi}\ m_{x}+\sin{\phi}\
m_{y}\right)\right)\\\ &-\left.\alpha I_{s}\sin{\phi}\ m_{y}\right]\\\
&+2\sqrt{D}\begin{pmatrix}d_{x}&d_{y}&d_{z}\end{pmatrix}\circ\begin{pmatrix}\dot{W}_{x}\\\
\dot{W}_{y}\\\ \dot{W}_{z}\end{pmatrix};\\\
\mathrm{where}&\begin{pmatrix}d_{x}\\\ d_{y}\\\
d_{z}\end{pmatrix}=\begin{pmatrix}Rm_{z}\mathcal{D}_{31}-m_{x}\mathcal{D}_{11}\\\
Rm_{z}\mathcal{D}_{32}-m_{x}\mathcal{D}_{12}\\\
Rm_{z}\mathcal{D}_{33}-m_{x}\mathcal{D}_{13}\end{pmatrix}.\end{split}$ (25)

It can be observed from Eq. (25) that each of damping, spin current and
thermal noise contribute to the rate of energy change. In addition, the last
term in Eq. (25) is simplified to
$\sqrt{d_{x}^{2}+d_{y}^{2}+d_{z}^{2}}\sqrt{\dot{W}_{x}^{2}+\dot{W}_{y}^{2}+\dot{W}_{Z}^{2}}$
which leads to

$\displaystyle\begin{split}&\frac{\partial
g_{L}}{\partial\tau}=-2\alpha\left[\left(1+R\right)m_{x}^{2}+g_{L}\left(R-g_{L}\right)\right.\\\
&+\frac{I_{s}}{\alpha}\left(\left(1+g_{L}\right)\cos{\phi}\
m_{x}+g_{L}\sin{\phi}\ m_{y}\right)\\\ \end{split}$
$\displaystyle\begin{split}&\left.+I_{s}\left(\sin{\phi}\
m_{x}m_{y}-R\left(\sin{\phi}\ m_{x}m_{z}-\cos{\phi}\
m_{y}m_{z}\right)\right)\right]\\\
&+2\sqrt{\frac{\alpha}{\Delta_{0}}}\sqrt{\left(1+R\right)m_{x}^{2}+g_{L}\left(R-g_{L}\right)}\circ\dot{W}_{g_{L}},\end{split}$
(26)

where $\dot{W}_{g_{L}}$ is 1D white Gaussian noise and it acts away from the
constant-energy orbit along a normal. (Pinna, Stein, and Kent, 2014)

We now average $\frac{\partial g_{L}}{\partial\tau}$ over a constant-energy
orbit since the non-conservative effects act over long time-scales and we are
interested in studying their effects on short time-scale rapid periodic
motion. In the limit of zero damping, zero spin current at absolute zero
temperature, the constant energy solutions to Eqs. (21)-(23) in the anti-
parallel well are given as (Mayergoyz, Bertotti, and Serpico, 2009; Pinna,
Kent, and Stein, 2013b; Pinna, Stein, and Kent, 2014)

$\displaystyle\begin{split}&m_{x}^{c}(t)=-\sqrt{\frac{R-g_{L}}{R+1}}\
dn\left[\sqrt{R-g_{L}}\ t,k^{2}\right],\\\ &m_{y}^{c}(t)=\sqrt{1+g_{L}}\
sn\left[\sqrt{R-g_{L}}\ t,k^{2}\right],\\\
&m_{z}^{c}(t)=\sqrt{\frac{1+g_{L}}{R+1}}\ cn\left[\sqrt{R-g_{L}}\
t,k^{2}\right],\end{split}$

where $k^{2}\equiv R\frac{1+g_{L}}{R-g_{L}}$, $0<k^{2}<1$ and $sn[\cdot]$,
$cn[\cdot]$, $dn[\cdot]$ are Jacobi elliptic functions. Magnetization
oscillates around the easy axis with a time period,
$T\left(g_{L}\right)=\frac{4}{\sqrt{R-g_{L}}}K\left(k^{2}\right)$, where
$K\left(k^{2}\right)$ is the complete elliptic integral of first kind. In Eq.
(A) only the time averages of $m_{x}$ and $m_{x}^{2}$ are non-zero while those
of $m_{y}$, $m_{x}m_{y}$, $m_{x}m_{z}$, and $m_{y}m_{z}$ are zero due to their
periodic nature. To evaluate these averages we geometrically parametrize
(Pinna, Kent, and Stein, 2013b; Pinna, Stein, and Kent, 2014) the constant-
energy orbit in the anti-parallel well as

$\displaystyle\begin{split}\cosh^{2}{(w)}-\sinh^{2}{(w)}&=1,\\\
\frac{1}{-g_{L}}m_{x}^{2}-\frac{R}{-g_{L}}m_{z}^{2}&=1\end{split}$

which leads to the form as given below

$\displaystyle\begin{split}m_{x}^{c}&=-\sqrt{-g_{L}}\cosh{(w)},\\\
m_{z}^{c}&=\pm\sqrt{\frac{-g_{L}}{R}}\sinh{(w)},\\\
m_{y}^{c}&=\pm\sqrt{\frac{R-g_{L}}{R}}\sqrt{1-\zeta^{2}\cosh^{2}{(w)}},\\\
\zeta^{2}&=-\frac{g_{L}\left(R+1\right)}{R-g_{L}},\\\ \end{split}$ (27)

where $\zeta^{2}=1-k^{2}$, and
$-\cosh^{-1}{(1/\zeta)}<w<\cosh^{-1}{(1/\zeta)}$. Defining time average of any
function $p\left(t\right)$ as

$\displaystyle\begin{split}\langle
p\left(t\right)\rangle&=\frac{4}{T\left(g_{L}\right)}\int_{0}^{t=m_{y}:\sqrt{1+g_{L}}\rightarrow
0}p\left(t\right)\ dt,\\\ \end{split}$ (28)

the time averages of $\langle m_{x}\rangle$ and $\langle m_{x}^{2}\rangle$ is
evaluated as

$\displaystyle\begin{split}\langle
m_{x}\rangle&=\frac{4}{T\left(g_{L}\right)}\int\limits_{0}^{\cosh^{-1}{(1/\zeta)}}m_{x}\left(w\right)\Bigg{|}\frac{\sqrt{R}m_{z}}{Rm_{y}m_{z}}\Bigg{|}\
dw\\\
&=-\frac{4}{T\left(g_{L}\right)}\frac{\sqrt{-g_{L}}}{\sqrt{R-g_{L}}}\int\limits_{\text{\hbox
to0.0pt{\hss$0$\hss}}}^{\text{\hbox
to0.0pt{\hss$\cosh^{-1}{(1/\zeta)}$\hss}}}\frac{\cosh{(w)}\
dw}{\sqrt{1-\zeta^{2}\cosh^{2}{(w)}}}\\\
&=-\frac{4}{T\left(g_{L}\right)}\frac{\sqrt{-g_{L}}}{\zeta\sqrt{R-g_{L}}}\int\limits_{\text{\hbox
to0.0pt{\hss$0$\hss}}}^{\text{\hbox
to0.0pt{\hss$\cosh^{-1}{(1/\zeta)}$\hss}}}\frac{\cosh{(w)}\
dw}{\sqrt{\frac{1-\zeta^{2}}{\zeta^{2}}-\sinh^{2}{(w)}}}\\\
&=-\frac{\pi}{2}\frac{{\sqrt{R+-g_{L}}}}{\sqrt{1+R}}\frac{1}{K\left(\zeta\right)}=-\frac{\pi}{2}\sqrt{\frac{R-g_{L}}{R+1}}\frac{1}{K\left(R,g_{L}\right)},\end{split}$
(29)

and

$\displaystyle\begin{split}\langle
m_{x}^{2}\rangle&=\frac{4}{T\left(g_{L}\right)}\int\limits_{0}^{\cosh^{-1}{(1/\zeta)}}m_{x}^{2}\left(w\right)\Bigg{|}\frac{m_{x}\sqrt{R}m_{z}}{m_{x}Rm_{y}m_{z}}\Bigg{|}\
dw\\\
&=\frac{-4}{T\left(g_{L}\right)}\frac{g_{L}}{\sqrt{R-g_{L}}}\int\limits_{0}^{\cosh^{-1}{(1/\zeta)}}\frac{\cosh^{2}{(w)}\
dw}{\sqrt{1-\zeta^{2}\cosh^{2}{(w)}}}\\\
&=\frac{-1}{K\left(\zeta\right)}\frac{g_{L}}{\sqrt{1-\zeta^{2}}}\int\limits_{0}^{\cosh^{-1}{(1/\zeta)}}\frac{\cosh^{2}{(w)}\
dw}{\sqrt{1-\frac{\zeta^{2}}{1-\zeta^{2}}\sinh^{2}{(w)}}}\\\
&=\frac{-1}{K\left(\zeta\right)}\frac{g_{L}}{\zeta}\int\limits_{0}^{1}\frac{\sqrt{1-\left(1-\frac{1}{\zeta^{2}}\right)u^{2}}\
du}{\sqrt{1-u^{2}}}\\\
&=\frac{-1}{K\left(\zeta\right)}\frac{g_{L}}{\zeta^{2}}E\left(1-\zeta^{2}\right)=\frac{R-g_{L}}{R+1}\frac{E\left(R,g_{L}\right)}{K\left(R,g_{L}\right)},\end{split}$
(30)

where we have used $dt=\frac{dp\left(w\right)/dw}{dp\left(t\right)/dt}dw$.
Substituting $\langle m_{x}\rangle$ and $\langle m_{x}^{2}\rangle$ in Eq. (A)
results in the average rate of energy flow as

$\displaystyle\begin{split}&\left\langle\frac{\partial
g_{L}}{\partial\tau}\right\rangle=-\frac{\pi\alpha}{K\left(R,g_{L}\right)}\sqrt{\frac{R-g_{L}}{1+R}}\left[-\frac{I_{s}}{\alpha}\cos{\phi}\left(1+g_{L}\right)\right.\\\
&+\left.\frac{2}{\pi}\sqrt{1+R}\sqrt{R-g_{L}}\left(E\left(R,g_{L}\right)+g_{L}K\left(R,g_{L}\right)\right)\right]\\\
&-2\sqrt{\frac{\alpha}{\Delta_{0}}}\sqrt{\frac{\left(R-g_{L}\right)}{K\left(R,g_{L}\right)}}\sqrt{E\left(R,g_{L}\right)+g_{L}K\left(R,g_{L}\right)}\circ\dot{W}_{g_{L}}.\end{split}$
(31)

In this paper, we have considered the spin polarization $\mathbf{\hat{n}}_{p}$
to be collinear to the easy axis $\mathbf{\hat{n}}_{e}$ so $\phi=0$ which
leads to Eq. (7).

##  Appendix B Solving the integral in Eq. (9)

In Eq. (9), substitute $R-g_{L}=u^{2}$, so that we have

$\displaystyle\begin{split}\tau_{s}&=\frac{1}{2A\alpha}\frac{R+4}{R+2}\int\limits_{\sqrt{R-g_{L_{i}}}}^{\sqrt{R-g_{L_{f}}}}\frac{u^{2}-\frac{R(1+R)}{R+4}}{u^{2}-\frac{R(1+R)}{R+2}}\left[\frac{du}{P(u)}\right],\end{split}$
(32)

where
$P(u)=u^{5}-(B/A+2R)u^{3}+\frac{\tilde{I}_{s}}{A\sqrt{1+R}}u^{2}+(R^{2}+(B/A)R+C/A)u-\tilde{I}_{s}\sqrt{1+R}/A$
and $w_{i}$ are the roots of the polynomial $P(u)$ which are evaluated
numerically. Using partial fractions to resolve the denominator of the
integral we finally have

$\displaystyle\begin{split}\tau_{s}&=\frac{1}{2\alpha\left(R+2\right)}\bBigg<EMAIL_ADDRESS>\end{split}$
$\displaystyle\begin{split}&+\left.\sqrt{\frac{R(1+R)}{R+2}}\left\\{\frac{\log{\left[\frac{\sqrt{(R-g_{L_{f}})(R+2)}-\sqrt{R(1+R)}}{\sqrt{(R-g_{L_{i}})(R+2)}-\sqrt{R(1+R)}}\right]}}{A\prod\limits_{n=1}^{5}\left(\frac{R(1+R)}{R+2}-w_{i}\right)}\right.\right.\\\
&+\bBigg<EMAIL_ADDRESS>(33)

##  Appendix C $R\to\infty$

For large values of R, Eq. (7) in the deterministic domain along with rational
approximations for the elliptic integrals can be simplified as

$\displaystyle\begin{split}\frac{\partial
g_{L}}{\partial\tau}&=4\alpha\left[\frac{1-g_{L}}{3-g_{L}}\right]\left[\tilde{I}_{s}\left(1+g_{L}\right)-R\left\\{g_{L}\frac{3-g_{L}}{2-2g_{L}}\right.\right.\\\
&\left.\left.+1-\frac{(1+g_{L})}{4}\frac{(1+g_{L})^{2}-28(1+g_{L})+64}{4(1+g_{L})^{2}-40(1+g_{L})+64}\right\\}\right].\end{split}$
(34)

Now substituting $1+g_{L}=x$ in the previous equation leads us to

$\displaystyle\begin{split}\frac{dx}{d\tau}&=4\alpha\left[\frac{x-2}{x-4}\right]\left[\tilde{I}_{s}x-R\left\\{\left(x-1\right)\frac{x-4}{2x-4}\right.\right.\\\
&\left.\left.+1-\frac{x}{4}\frac{x^{2}-28x+64}{4x^{2}-40x+64}\right\\}\right],\end{split}$
(35)

which can then be simplified and rearranged to an integral of the form

$\displaystyle\begin{split}\int\limits_{1+g_{L_{i}}}^{1+g_{L_{f}}}\frac{(x-4)(x-8)dx}{x\left[x^{2}-Ex+F\right]}&=\frac{\alpha}{4}\left(16\tilde{I}_{s}-7R\right)\tau_{s},\end{split}$
(36)

where $E=\frac{160\tilde{I}_{s}-60R}{16\tilde{I}_{s}-7R}$ and
$F=\frac{256\tilde{I}_{s}-128R}{16\tilde{I}_{s}-7R}$. If the roots of the
quadratic equation is $x^{2}-Ex+F=0$ are a and b then the switching time
$\tau_{s}$ can be evaluated as

$\displaystyle\begin{split}\tau_{s}&=\frac{1}{2\alpha\left[\tilde{I}_{s}-0.5R\right]}\left[\log{\left[\frac{1+g_{L_{f}}}{1+g_{L_{i}}}\right]}\right.\\\
&+\left.\frac{b(a-4)(a-8)}{32(a-b)}\log{\left[\frac{1+g_{L_{f}}-a}{1+g_{L_{i}}-a}\right]}\right.\\\
&-\left.\frac{a(b-4)(b-8)}{32(a-b)}\log{\left[\frac{1+g_{L_{f}}-b}{1+g_{L_{i}}-b}\right]}\right].\end{split}$
(37)

##  Appendix D Uniaxial Limit

In the uniaxial limit Eq. (7) in the deterministic domain reduces to

$\displaystyle\begin{split}\left\langle\frac{\partial
g_{L}}{\partial\tau}\right\rangle&=2\alpha\sqrt{-g_{L}}\left(1+g_{L}\right)\left(\tilde{I}_{s}-\sqrt{-g_{L}}\right),\end{split}$
(38)

which can then be integrated as

$\int\limits_{g_{L_{i}}}^{g_{L_{f}}}\frac{dg_{L}}{\sqrt{-g_{L}}\left(1+g_{L}\right)\left(\tilde{I}_{s}-\sqrt{-g_{L}}\right)}=2\alpha\tau_{s}.$

Substituting $g_{L}=-u^{2}$, and using partial fractions to separate the terms
in the denominator and integrating with proper limits leads us to

$\begin{split}\tau_{s}&=\frac{1}{2\alpha\left(\tilde{I}_{s}^{2}-1\right)}\left[\tilde{I}_{s}\left\\{\log{\left[\frac{1+\sqrt{-g_{L_{i}}}}{1+\sqrt{-g_{L_{f}}}}\right]}\right.\right.\\\
&\left.\left.-\log{\left[\frac{1-\sqrt{-g_{L_{i}}}}{1-\sqrt{-g_{L_{f}}}}\right]}\right\\}-\log{\left[\frac{1+g_{L_{i}}}{1+g_{L_{f}}}\right]}\right.\\\
&+\left.2\log{\left[\frac{\frac{\tilde{I}_{s}}{\alpha}-\sqrt{-g_{L_{i}}}}{\frac{\tilde{I}_{s}}{\alpha}-\sqrt{-g_{L_{f}}}}\right]}\right].\end{split}$
(39)

##  Appendix E Eqs. (14) and (16)

Using $R=0$ in Eq. (13), the probability for a random variable
$G_{L}=-M_{x}^{2}$ is obtained as

$\displaystyle\begin{split}P_{G_{L}}\left[g_{L}\right]&=\mathbb{P}\left[G_{L}\leq
g_{L}\right]=\mathbb{P}\left[-M_{x}^{2}\leq g_{L}\right]\\\
&=1-\mathbb{P}\left[|M_{x}|\leq\sqrt{-g_{L}}\right]\\\
&=1-\mathbb{P}\left[-\sqrt{-g_{L}}\leq M_{x}\leq\sqrt{-g_{L}}\right]\\\
&=1-\frac{2}{Z(\Delta_{0},0)}\int\limits_{0}^{\sqrt{-g_{L}}}dm_{x}\
\exp\left(\Delta_{0}m_{x}^{2}\right)\\\ \end{split}$
$\displaystyle\begin{split}&=1-\frac{2}{Z(\Delta_{0},0)\sqrt{\Delta_{0}}}\int\limits_{0}^{\sqrt{-g_{L}}}dm_{x}\
\sqrt{\Delta_{0}}\ \exp\left(\Delta_{0}m_{x}^{2}\right)\\\
&=1-\frac{2\exp{(\Delta_{0}g_{L})}}{Z(\Delta_{0},0)\sqrt{\Delta_{0}}\exp(\Delta_{0}g_{L})}\int\limits_{0}^{\sqrt{-\Delta_{0}g_{L}}}du\
\exp\left(u^{2}\right)\\\
&=1-\frac{2\exp(-\Delta_{0}g_{L})}{Z(\Delta_{0},0)\sqrt{\Delta_{0}}}F[\sqrt{-\Delta_{0}g_{L}}].\end{split}$
(40)

In step 6 above we have substituted $\sqrt{\Delta_{0}}m_{x}=u$, and used
$F[x]=\exp(-x^{2})\int_{0}^{x}dy\exp(y^{2})$, the Dawson’s integral in step 7.
Also, $P_{G_{L}}[-1]=0$, therefore,

$Z(\Delta_{0},0)=\frac{2\exp(\Delta_{0})}{\sqrt{\Delta_{0}}}F[\sqrt{\Delta_{0}}].$
(41)

Eqs. (40) and (41) together lead to Eq. (16).

In order to arrive at the probability distribution function (pdf) we now
differentiate Eq. (40) with respect to $g_{L}$ at step 4 to get

$\displaystyle\begin{split}\rho_{G_{L}}[g_{L}]&=\frac{1}{Z(\Delta_{0},0)}\frac{\exp\left(-\Delta_{0}g_{L}\right)}{\sqrt{-g_{L}}}\\\
&=\frac{\sqrt{\Delta_{0}}}{2F[\sqrt{\Delta_{0}}]}\frac{\exp\left(-\Delta_{0}\left(1+g_{L}\right)\right)}{\sqrt{-g_{L}}}.\end{split}$
(42)

We have used $Z(\Delta_{0},R)=Z(\Delta_{0},0)$ and have found that all the
results match well.

## References

  * Kent and Worledge (2015) A. D. Kent and D. C. Worledge, Nature nanotechnology 10, 187 (2015).
  * Sadovnikov  _et al._ (2017) A. V. Sadovnikov, S. A. Odintsov, E. N. Beginin, S. E. Sheshukova, Y. P. Sharaevskii, and S. A. Nikitov, Physical Review B 96, 144428 (2017).
  * Chen  _et al._ (2016) T. Chen, R. K. Dumas, A. Eklund, P. K. Muduli, A. Houshang, A. A. Awad, P. Dürrenfeld, B. G. Malm, A. Rusu, and J. Åkerman, Proceedings of the IEEE 104, 1919 (2016).
  * Slonczewski (1996) J. C. Slonczewski, Journal of Magnetism and Magnetic Materials 159, L1 (1996).
  * Berger (1996) L. Berger, Physical Review B 54, 9353 (1996).
  * Sinova  _et al._ (2015) J. Sinova, S. O. Valenzuela, J. Wunderlich, C. H. Back, and T. Jungwirth, Reviews of Modern Physics 87, 1213 (2015).
  * Liu  _et al._ (2012a) L. Liu, C.-F. Pai, Y. Li, H. Tseng, D. Ralph, and R. Buhrman, Science 336, 555 (2012a).
  * Sander (2004) D. Sander, Journal of Physics: Condensed Matter 16, R603 (2004).
  * Liu  _et al._ (2012b) L. Liu, O. J. Lee, T. J. Gudmundsen, D. C. Ralph, and R. A. Buhrman, Physical review letters 109, 096602 (2012b).
  * Oh  _et al._ (2016) Y.-W. Oh, S.-h. C. Baek, Y. Kim, H. Y. Lee, K.-D. Lee, C.-G. Yang, E.-S. Park, K.-S. Lee, K.-W. Kim, G. Go, _et al._ , Nature nanotechnology 11, 878 (2016).
  * Fukami  _et al._ (2016) S. Fukami, T. Anekawa, C. Zhang, and H. Ohno, Nature Nanotechnology 11, 621 (2016).
  * Sun (2000) J. Z. Sun, Physical Review B 62, 570 (2000).
  * Pinna, Kent, and Stein (2013a) D. Pinna, A. Kent, and D. Stein, Journal of Applied Physics 114, 033901 (2013a).
  * d’Aquino, Scalera, and Serpico (2019) M. d’Aquino, V. Scalera, and C. Serpico, Journal of Magnetism and Magnetic Materials 475, 652 (2019).
  * Apalkov and Visscher (2005) D. M. Apalkov and P. B. Visscher, Physical Review B 72, 180405 (R) (2005).
  * Taniguchi  _et al._ (2013) T. Taniguchi, Y. Utsumi, M. Marthaler, D. S. Golubev, and H. Imamura, Physical Review B 87, 054406 (2013).
  * Mayergoyz, Bertotti, and Serpico (2009) I. D. Mayergoyz, G. Bertotti, and C. Serpico,  _Nonlinear magnetization dynamics in nanosystems_ (Elsevier, 2009).
  * Pinna, Kent, and Stein (2013b) D. Pinna, A. D. Kent, and D. L. Stein, Physical Review B 88, 104405 (2013b).
  * Newhall and Vanden-Eijnden (2013) K. A. Newhall and E. Vanden-Eijnden, Journal of Applied Physics 113, 184105 (2013).
  * Pinna, Stein, and Kent (2014) D. Pinna, D. L. Stein, and A. D. Kent, Physical Review B 90, 174405 (2014).
  * Yu, Muratov, and Moore (2019) Y. Yu, C. B. Muratov, and R. O. Moore, IEEE Transactions on Magnetics (2019).
  * Kani (2017) N. Kani,  _Modeling of magnetization dynamics and applications to spin-based logic and memory devices_ , Ph.D. thesis, Georgia Institute of Technology (2017).
  * Vincent  _et al._ (2014) A. F. Vincent, N. Locatelli, J.-O. Klein, W. S. Zhao, S. Galdin-Retailleau, and D. Querlioz, IEEE Transactions on Electron Devices 62, 164 (2014).
  * Hollos (2012) S. Hollos, “Rational approximations for complete elliptic integrals,” http://www.exstrom.com/math/elliptic/ellipint.html (2012), accessed: 2019-07-31.
  * Liu  _et al._ (2014) H. Liu, D. Bedau, J. Sun, S. Mangin, E. Fullerton, J. Katine, and A. Kent, Journal of Magnetism and Magnetic Materials 358, 233 (2014).
  * Kani, Rakheja, and Naeemi (2016) N. Kani, S. Rakheja, and A. Naeemi, IEEE Transactions on Electron Devices 63, 4119 (2016).
  * inv (2016) “A file to numerically evaluate the inverse of a function in python,” https://github.com/alvarosg/pynverse/blob/master/pynverse/inverse.py (2016).
  * Sbiaa  _et al._ (2009) R. Sbiaa, R. Law, E.-L. Tan, and T. Liew, Journal of Applied Physics 105, 013910 (2009).
  * Zhao  _et al._ (2011) H. Zhao, A. Lyle, Y. Zhang, P. Amiri, G. Rowlands, Z. Zeng, J. Katine, H. Jiang, K. Galatsis, K. Wang, _et al._ , Journal of Applied Physics 109, 07C720 (2011).
  * Ament  _et al._ (2016) S. Ament, N. Rangarajan, A. Parthasarathy, and S. Rakheja, arXiv preprint arXiv:1607.04596 (2016).