# Anisotropic Hall effects in Bi2Se3/EuS interfaces

Juhi Singh  Department of Physics, Indian Institute of Technology Roorkee,
Roorkee 247667, India  Karthik V. Raman<EMAIL_ADDRESS>Tata Institute
of Fundamental Research, Hyderabad, Telengana 500046, India  Narayan Mohanta
<EMAIL_ADDRESS>Department of Physics, Indian Institute of
Technology Roorkee, Roorkee 247667, India

###### Abstract

Proximity coupling of ferromagnetic insulator EuS to the topological insulator
Bi2Se3 has been proposed to break time-reversal symmetry near the surface of
Bi2Se3, introducing an energy gap or a tilt in the surface Dirac cone. As an
inverse proximity effect, strong spin-orbit coupling available in the
topological surface states can enhance the Curie temperature of ferromagnetism
in EuS largely beyond its bulk value, and also generate a magnetic anisotropy.
This can result in a canting of the magnetic moment of Eu ions in a plane
perpendicular to the interface. Here, we investigate theoretically electronic
transport properties arising from the Bi2Se3/EuS interfaces in the planar Hall
geometry. Our analysis, based on a realistic model Hamiltonian and a semi-
classical formalism for the Boltzmann transport equation, reveals distinct
intriguing features of anisotropic planar Hall conductivity, depending on
different scenarios for the canting of the Eu moments: fixed Eu moment
canting, and freely-orientable Eu moment in response to the external in-plane
magnetic field. The anisotropy in the planar Hall conductivity arises from the
asymmetric Berry curvature of the gapped topological surface states. We also
explore topological Hall effect of the Dirac surface states, coupled to a
skyrmion crystal which can emerge in the EuS due to the interplay of
ferromagnetic Heisenberg exchange, interfacial Dzyaloshinskii-Moriya
interaction, and perpendicular alignment of the Eu moment. Our study provides
new impetus for probing complex interplay between magnetic exchange
interactions and topological surface states via anisotropic planar Hall
effects.

##  I Introduction

Three-dimensional topological insulators, when interfaced with a ferromagnet
or a superconductor, are known to generate a wide variety of intriguing
phenomena such as quantum anomalous Hall effect Chang  _et al._ (2013, 2023);
Sun  _et al._ (2019); Pan  _et al._ (2020); Moon  _et al._ (2019), topological
magnetoelectric effect Morimoto  _et al._ (2015); Dziom  _et al._ (2017) and
putative topological superconductivity Wang  _et al._ (2012); Shen  _et al._
(2017). Particularly, by interfacing Bi2Se3 with EuS, a ferromagnetic
insulator with partially filled 4$f$ orbitals and Curie temperature
$T_{c}$$\sim$$17$ K, a strong proximity effect was observed at the interface.
The $T_{c}$ in EuS film was found to enhance beyond its bulk value Katmis  _et
al._ (2016), attributed to the large enhancement in exchange coupling in EuS
near the interface mediated by the Ruderman-Kittel-Kasuya-Yosida type
interactions Kim  _et al._ (2017); Mathimalar  _et al._ (2020). As an inverse
proximity effect, the presence of topological surface states of Bi2Se3, having
strong spin-orbit coupling, modify spin ordering in EuS. As the thickness of
EuS is increased, the induced out-of-plane magnetic anisotropy turns in-plane
Katmis  _et al._ (2016); Kim  _et al._ (2017); Mathimalar  _et al._ (2020);
Assaf  _et al._ (2015); Wei  _et al._ (2013); Lee  _et al._ (2016). As a
result, the Eu moment in a thin film of EuS acquires a canting out of the
interface plane. The magnetic exchange field due to the canted Eu moment lifts
the degeneracy of the surface Dirac bands, opening up an anisotropic energy
gap in the surface Dirac bands.

Planar Hall effect (PHE), arising from a resistivity anisotropy induced by an
applied in-plane magnetic field, was found prominently on the surface of
Bi2Se3 at energies near the Dirac point Taskin  _et al._ (2017). The origin of
this resistivity anisotropy in Bi2Se3 is not very well understood. Berry
curvature of bulk conduction bands was theoretically proposed to produce a
PHE, giving rise to a standard $\cos{\theta_{in}}\sin{\theta_{in}}$ variation
of the planar Hall conductivity $\sigma_{xy}$ with the in-plane applied field
angle $\theta_{in}$ Nandy  _et al._ (2018). The in-plane magnetic field shifts
the Dirac point on the interfacing Bi2Se3 surface relative to the Dirac point
on the bottom surface. Proposals for non-topological origin of the PHE include
anisotropic back-scattering by magnetic disorders Trushin  _et al._ (2009);
Chiba  _et al._ (2017) and by the tilted Dirac cone Zheng  _et al._ (2020). In
Bi2Te3, a rise in the PHE with the thickness of the sample was found, implying
that the bulk states also contribute to the PHE Bhardwaj  _et al._ (2021).
Subsequent experiments on (Bi, Sb)2Te3/EuS interfaces reported unconventional
PHE which requires the consideration of the large proximity effect at the
interface Rakhmilevich  _et al._ (2018). Signatures of topological spin
textures such as magnetic skyrmions have also been found in PHE measurements
Suri  _et al._ .

Figure 1: (a) Planar Hall measurement geometry for the Bi2Se3/EuS interface
system, considered in this paper. Electric field ${\bf E}$ (due to applied
charge current) is oriented along the $x$ axis. The applied magnetic field
${\bf B}$ makes an angle $\theta_{in}$ with ${\bf E}$ in the $x$-$y$ plane.
The magnetization in EuS is inclined at an angle $\theta_{m}$ with the
$z$-axis and $\phi_{m}$ with the $x$ axis in the $x$-$y$ plane. (b) Band
dispersion of Bi2Se3 slab of seven quintuple layers ($N\\!=\\!7$), in the
absence of any magnetization field, showing two degenerate Dirac surface bands
in red and doubly-degenerate bulk bands in green. (c) Band dispersion of
Bi2Se3/EuS interface system with $N\\!=\\!7$ in Bi2Se3 slab,
$m\\!=\\!0.15~{}{\rm eV}$ (magnetic moment induced at the interface due to
proximity effect), $\theta_{m}=45\degree$, $B\\!=\\!2~{}$T, and
$\theta_{in}=45\degree$. (d),(e) Berry curvature and orbital magnetic moment
for the lower gapped surface band (blue band in 1(c)) with parameters
$m\\!=\\!6.9~{}\rm{\mu_{B}}$, $\theta_{m}\\!=\\!45\degree$, $B\\!=\\!2~{}$T
and $\theta_{in}=45\degree$. In (d) and (e), the colormap shows the
$z$-components, while the $x$ and $y$-components are shown by the arrows.

In this paper, motivated by these recent experimental observations, we
investigate the magnetic proximity effect at the Bi2Se3/EuS interface and the
influence of the canting of the Eu moment on the resulting PHE. We find
anisotropic PHE of different features, depending on the interplay of the Eu
moment and the applied in-plane magnetic field. We use a realistic slab
Hamiltonian for the interface system to calculate the planar Hall
conductivity, within the framework of semi-classical Boltzmann formalism,
originating from the Berry curvature of the anisotropically-gapped Dirac cone
on the interfacing Bi2Se3 surface. We find that the PHE appears predominantly
within a field range between two critical values. Furthermore, we investigate
Hall response from magnetic skyrmions which can emerge at the interface when
the Eu moment anisotropy is perpendicular to the interface plane and give rise
to a finite scalar spin chirality, another source of the PHE. In this case, we
find that with increasing the in-plane magnetic field amplitude, there is a
gradual transition from the skyrmion crystal phase to the in-plane
ferromagnetic phase and the topological Hall conductivity decreases
monotonically.

The remainder of this paper is organized as follows. In Sec. II, we describe
our theoretical model for the Bi2Se3/EuS interface and the method to calculate
the planar Hall conductivities within semi-classical Boltzmann transport
formalism. In Sec. III, we present our numerical results of the anisotropic
PHE which can originate from the momentum-space Berry curvature, for different
orientation of the Eu moment. In Sec. IV, we show the planar topological Hall
signature at the interface due to skyrmion spin texture on the EuS side of the
interface. In Sec. V, we provide further discussion on our results regarding
the connection to possible experimental observation, and summarize our
results.

##  II Model and Method

The planar Hall geometry for the Bi2Se3/EuS interface, with canted spin
moments in EuS, is shown schematically in Fig. 1(a). The interface properties
can be accessed via the Hall probes using contacts bored through EuS Suri  _et
al._ . Bi2Se3 appears in rhombohedral crystal structure, with space group
$\mathrm{D}^{5}_{3\mathrm{d}}(\mathrm{R}\bar{3}\mathrm{m})$, in which a five-
layer building block consisting of inequivalent Bi and Se atoms, known as the
‘quintuple’ layer repeats along perpendicular to these layers Zhang  _et al._
(2009). We theoretically model Bi2Se3 of a finite thickness by considering a
slab geometry, in which open boundary conditions are imposed along the
stacking direction of the quintuple layers, and periodic boundary conditions
along the other two orthogonal directions within each quintuple layer.
Therefore, each quintuple layer of Bi2Se3 can be represented by a two-
dimensional momentum-space Hamiltonian, and these Hamiltonians are coupled
along the stacking direction (considered to be the $z$ axis) in real
coordinate space. To incorporate the magnetic proximity effect due to EuS thin
film and the influence of the applied in-plane magnetic field, Zeeman exchange
coupling is introduced in the Hamiltonian for the topmost layer of ${\rm
Bi_{2}Se_{3}}$.

The total slab Hamiltonian for $N$ quintuple layers is given by Mohanta  _et
al._ (2017)

$\displaystyle H({\bf k}_{\parallel})=\begin{pmatrix}H_{0}+H_{B}&H_{1}&&&&&\\\
H^{{\dagger}}_{1}&H_{0}&H_{1}&&&&\\\ &H^{{\dagger}}_{1}&H_{0}&H_{1}&&&\\\
&&.&&&&\\\ &&&.&&&\\\ &&&&.&&\\\ &&&&H^{{\dagger}}_{1}&H_{0}&H_{1}\\\
&&&&&H^{{\dagger}}_{1}&H_{0}\\\ \end{pmatrix},$ (1)

which is written in the basis $\ket{P1^{+}_{z},\uparrow,l_{z}}$,
$\ket{P2^{-}_{z},\uparrow,l_{z}}$, $\ket{P1^{+}_{z},\downarrow,l_{z}}$,
$\ket{P2^{-}_{z},\downarrow,l_{z}}$, where $P1^{\pm}_{z}$ and $P2^{\pm}_{z}$
denote orbitals originated from the $p$ orbitals of Bi and Se Liu  _et al._
(2010), ${\bf k}_{\parallel}\\!\equiv\\!(k_{x},k_{y})$,
$(\downarrow,\uparrow)$ denotes the electron spin quantization states, and
$l_{z}$ is the layer index which varies between 0 and $N$. The Hamiltonians
$H_{0}$, $H_{1}$ and $H_{B}$ describe, respectively, a single quintuple layer,
coupling between two neighboring quintuple layers, and Zeeman exchange
coupling on the top surface layer due to proximity effect of ferromagnetic EuS
and due to the applied in-plane magnetic field. These Hamiltonians can be
written as

$\displaystyle H_{0}=\begin{pmatrix}\epsilon_{0}+M&0&0&A_{0}{\bf k_{-}}\\\
0&\epsilon_{0}-M&A_{0}{\bf k_{-}}&0\\\ 0&A_{0}{\bf k_{+}}&\epsilon_{0}+M&0\\\
A_{0}{\bf k_{+}}&0&0&\epsilon_{0}-M\\\ \end{pmatrix},$ (2)

$\displaystyle H_{1}=\begin{pmatrix}-M_{1}-C_{1}&iB_{0}/2&0&0\\\
iB_{0}/2&M_{1}-C_{1}&0&0\\\ 0&0&-M_{1}-C_{1}&-iB_{0}/2\\\
0&0&-iB_{0}/2&M_{1}-C_{1}\\\ \end{pmatrix},$ (3)

$\displaystyle H_{B}=\begin{pmatrix}h_{z}&0&h_{x}-ih_{y}&0\\\
0&h_{z}&0&h_{x}-ih_{y}\\\ h_{x}+ih_{y}&0&-h_{z}&0\\\
0&h_{x}+ih_{y}&0&-h_{z}\\\ \end{pmatrix},$ (4)

where
$\epsilon_{0}(\mathbf{k}_{\parallel})\\!=\\!C_{0}+2C_{1}+2C_{2}(2-\cos(k_{x}a)-\cos(k_{y}a))$,
$M_{0}(\mathbf{k}_{\parallel})\\!=\\!M_{0}+2M_{1}+2M_{2}(2-\cos(k_{x}a)-\cos(k_{y}a))$,
$a$ is the lattice constant, ${\bf k_{\pm}}=\sin{(k_{x}a)}\pm
i\sin{(k_{y}a)}$, $h_{x}\\!=\\!m_{x}+g\mu_{B}B_{x}/2$,
$h_{y}\\!=\\!m_{y}+g\mu_{B}B_{y}/2$, and $h_{z}\\!=\\!m_{z}$.
$m_{x}\\!=\\!m\sin\theta_{m}\cos\phi_{m}$,
$m_{y}\\!=\\!m\sin\theta_{m}\sin\phi_{m}$, $m_{z}\\!=\\!m\cos\theta_{m}$,
$B_{x}\\!=\\!B\cos\theta_{in}$, and $B_{y}\\!=\\!B\sin\theta_{in}$, ${\bf
m}\\!\equiv(m_{x},m_{y},m_{z})$ represents the $\mathrm{EuS}$ magnetic moment
which is assumed to be canted at a polar angle $\theta_{m}$ with the $z$-axis
and azimuithal angle $\phi_{m}$ with the $x$-axis on the $x$-$y$ (interface)
plane, and $B$ is the magnitude of the external in-plane magnetic field
applied at an angle $\theta_{in}$ with the $x$-axis on the $x$-$y$ plane. The
parameter values used in this paper are $A_{0}=0.8$ eV, $B_{0}=0.32$ eV,
$C_{0}=-0.0083$ eV, $C_{1}=0.024$ eV, $C_{2}=1.77$ eV, $M_{0}=-0.28$ eV,
$M_{1}=0.216$ eV, $M_{2}=2.6$ eV, $g=20$, and $a=4.14{\AA}$ Mohanta  _et al._
(2017). Fig. 1(b) and (c) show the band dispersions of a Bi2Se3 slab of seven
quintuple layers, without and with the Zeeman exchange coupling due to the Eu
moment and external applied magnetic field. The combined effect of the applied
in-plane magnetic field and the canted Eu moment, anisotropically lifts the
surface bands on the top surface of Bi2Se3, as shown in Fig. 1(c).

Figure 2: Planar Hall conductivity, $\sigma_{xy}$ normalized with respect to
its maximum value, is plotted for (a) $m=0.0~{}\mu_{B}$, (b)-(c)
$m=6.9~{}\mu_{B}$, $\phi_{m}=0\degree$, (d) $m=6.9$ to $0.4\mu_{B}$,
$\theta_{m}=\phi_{m}=0\degree$, and (e)-(f) $m=6.9~{}\mu_{B}$,
$\phi_{m}=0\degree$. The external in-plane field $B$ is fixed to $30~{}{\rm
mT}$ for plots (a)-(d) and to $B=0.5~{}{\rm T}$ for (e)-(f). The polar canting
angle $\theta_{m}$ is varied from $90\degree$ (in-plane) to $0\degree$ (out-
of-plane). The PHC is symmetric i.e., $\sigma_{xy}(B)=\sigma_{xy}(-B)$ with
$\pi$ period and follows a $\sin\theta_{in}\cos\theta_{in}$ for
$m=0.0~{}\mu_{B}$, while an angular shift is observed for
$(m_{x},my,m_{z})\equiv(0,0,m)$ as shown in (a) and (d). For cases with Eu
moment aligned in either $x$ or $y$ direction, antisymmetric behavior
$(\sigma_{xy}(B)=-\sigma_{xy}(-B))$ is observed with $2\pi$ period. For
$0\degree<\theta_{m}<90\degree$, the PHC reveals an anisotropic behavior with
$A\sin\theta_{in}\cos\theta_{in}+B\cos\theta_{in}$ dependence, where
$A,B\in\mathbb{R}$.  Figure 3: (a)-(b) The x-component of Berry curvature
$(\Omega_{x})$ at $(k_{x},k_{y})=(0,0)$ plotted as a function of in-plane
field $B$ for different values of Eu moment $m$ and azimuthal canting angle
$\phi_{m}$. For (a) $\phi_{m}=0\degree$ and the external field angle
$\theta_{in}$ and the polar canting angle $\theta_{m}$ for (a) and (b) are
both fixed at $45\degree$. (c) $\theta_{m}$ plotted as a function of the in-
plane magnetic field amplitude $B$, $\theta_{m}=\sin^{-1}(B/B_{c})$ as used in
the calculations. For $B\geq B_{c}(=20~{}{\rm mT})$, $\theta_{m}~{}{\rm
equals}~{}90\degree$ and the Eu moment turns in-plane. (d)-(f) depicts the
variation of the $x$\- component of Berry curvature $(\Omega_{x})$ at
$(k_{x},k_{y})=(0,0)$ with $B_{in}$ for different values of the Eu moment $m$,
initial canting angle $\theta_{m0}$ and azimuthal canting angle $\phi_{m}$.
The Berry curvature is maximum at a field value close to $B_{c}$ and then
drops to zero at $B_{in}\\!=\\!B_{c}$. In these plots, the polar canting angle
of the Eu moment varies with the applied field i.e.
$\sin\theta_{m}\\!=\\!B/B_{c}$. Figure 4: Planar Hall conductivity,
$\sigma_{xy}$ plotted for different magnitudes of Eu moment, (a)
$m=6.9~{}\mu_{B}$, (b) $m=3.5~{}\mu_{B}$, (c) $m=1.8~{}\mu_{B}$, and (d)
$m=0.9~{}\mu_{B}$. The polar canting angle $\theta_{m}$ follows
$\theta_{m}=\sin^{-1}(B/B_{c})$ relation and equals $90\degree$ for $B\geq
B_{c}$. The azimuthal canting angle $\phi_{m}$ is fixed to $0\degree$ and the
initial configuration is assumed to be in the perpendicular direction to the
interface i.e., $\theta_{m0}=0\degree$. The second critical field $B_{c}$ is
set to $20~{}{\rm mT}$. The first critical field, $B_{c^{\prime}}$ decreases
with decreasing $m$ and the width $\Delta B=|B_{c}-B_{c^{\prime}}|$ widens.
Other parameters are $N=7$, $T=24{\rm K}$, and $\tau=10^{-13}s$.

The planar Hall conductivity (PHC) is obtained within the semi-classical
Boltzmann formalism (see Appendix C for details) as

$\begin{split}\sigma_{xy}=&-\frac{e^{2}}{4\pi^{2}}\int d{\bf
k}~{}D\tau~{}\Big{[}\tilde{v}_{y}\tilde{v}_{x}+\frac{eB\cos\theta_{in}~{}\tilde{v}_{y}}{\hbar}(\bf{\tilde{v}_{k}}\cdot\Omega(\bf{k}))\\\
&+\frac{eB\sin\theta_{in}~{}\tilde{v}_{x}}{\hbar}({\bf\tilde{v}_{k}}\cdot{\bf\Omega(k)})\\\
&+\frac{e^{2}B^{2}\cos\theta_{in}\sin\theta_{in}}{\hbar^{2}}({\bf\tilde{v}_{k}}\cdot{\bf\Omega(k)})^{2}\Big{]}\partial_{\tilde{\epsilon}}\tilde{f}_{eq},\end{split}$
(5)

where ${\bf\tilde{v}_{k}}\equiv(\tilde{v}_{x},\tilde{v}_{y},0)$ is the
velocity of electrons at the surface of Bi2Se3 which interfaces EuS, modified
by the orbital magnetic moment (OMM),
${\bf\Omega(k)}\equiv(\Omega_{x},\Omega_{y},\Omega_{z})$ is the Berry
curvature, and $\tilde{f}_{eq}$ represents equilibrium Fermi-Dirac
distribution with modified energy(see Appendices C-D for details). The first
term on the right-hand side of Eq. (5), i.e., $\tilde{v}_{x}\tilde{v}_{y}$, is
the usual planar Hall contribution due to the presence of in-plane electric
and magnetic fields. Due to the magnetic proximity effect of EuS film on the
upper surface of Bi2Se3, a finite Berry curvature is induced in the surface
band, allowing the other three terms to contribute considerably. When Zeeman
coupling is perpendicular to the interface i.e., in the $z$-direction, the
z-component of Berry curvature, $\Omega_{z}$ is circularly symmetric centered
at $(k_{x},k_{y})=(0,0)$ in the $k$-space. With the introduction of in-plane
field components, the gapped Dirac point on the interfacing Bi2Se3 surface is
shifted in momentum compared to the gap-less Dirac point on the opposite
Bi2Se3 surface. Consequently, the maximum of $\Omega_{z}$ shifts from the
center of the Brillouin zone. Fig. 1(d) and (e) show the Berry curvature and
OMM of the anisotropically-gapped surface band, shown in Fig. 1(c). We note
that the Berry curvature and the OMM show textures analogous to anti-meron.

##  III Berry curvature driven planar Hall effect

The conventional planar Hall effect in a topological insulator is symmetric
with respect to the applied in-plane magnetic field and follows a
$\sin\theta_{in}\cos\theta_{in}$ dependence as shown in Fig. 2(a). However, at
the $\rm{Bi}_{2}\rm{Se}_{3}/\rm{EuS}$ interfaces, the planar Hall conductivity
will depend on the orientation of the Eu moment at the interface, which is
governed by the magnitude and direction of the external in-plane magnetic
field. In the following discussion, we have considered different
configurations of the Eu moment to explore different possible planar Hall
signatures at this interface.

###  III.1 Fixed orientation of the Eu moment

In this case, the Eu moment is assumed to be independent of the external in-
plane magnetic field. We consider that the azimuthal canting angle $\phi_{m}$
of the Eu moment is fixed at $0\degree$ i.e. the Eu moment is on the $x-z$
plane, while the polar canting angle is decreased from $90\degree$ (in-plane)
gradually with increasing the in-plane field to $0\degree$ (out-of-plane).
Also, the magnitude of the Eu moment is fixed at $6.9~{}\mu_{B}$, the magnetic
moment of bulk EuS as reported in Ref. Kim  _et al._ (2017). Fig. 2(b)-(c),
depicts the PHC for $\theta_{m}$ varying from $90\degree$ (Eu moment aligned
in the $x$-direction) to $3\degree$, an anti-symmetric contribution is
observed with $2\pi$ periodicity, which transitions to symmetric contribution
with $\pi$ periodicity for $\theta_{m}=0\degree$ as shown in Fig. 2(d). If the
Eu moment is considered on the $y-z$ plane i.e., $\phi_{m}=90\degree$, the
dependence shifts by $90\degree$ from $\sin\theta_{in}$ to
$\sin(90+\theta_{in})=\cos\theta_{in}$. Since $x$ and $y$ directions are
giving similar results with a $90\degree$ phase shift. Therefore, we have
considered the easy-axis alignment in the $x$-$z$ plane $(\phi_{m}=0)$ for
further results. From the numerical results, we find that the $\theta_{in}$
dependence of the PHC changes from $\sin\theta_{in}$ for
$\theta_{m}=90\degree$ to $-\cos\theta_{in}$ for $\theta_{m}=60\degree$, and
$\cos\theta_{in}$ at $\theta_{m}=30\degree$. Furthermore, for $\theta_{m}\leq
5\degree$, the PHC behaves as $\cos\theta_{in}-\sin\theta_{in}\cos\theta_{in}$
and eventually follows $-\sin\theta_{in}\cos\theta_{in}$ dependence for
$\theta_{m}=0\degree$. On generalizing the results for
$0\degree<\theta_{m}\leq 90\degree$, the PHC reveals an anisotropic behavior
with $A\sin\theta_{in}\cos\theta_{in}+B\cos\theta_{in}$ dependence, where
$A,B\in\mathbb{R}$. The transition from $\sin\theta_{in}$ to
$-\sin\theta_{in}\cos\theta_{in}$ is visible at very small $\theta_{m}$ for
low in-plane magnetic field. However, if the in-plane magnetic field is
increased to $0.5~{}\mathrm{T}$, this transition begins at higher polar
canting angles, as shown in Fig. 2(e) and (f). This happens because, at a
higher in-plane field, the effect of the magnetic moment due to EuS is
suppressed by the external field. Now, varying the magnitude of Eu moment $m$
from $0$ to $6.9~{}\mu_{B}$ for $\theta_{m}=\phi_{m}=0\degree$, the PHC shifts
from $\sin\theta_{in}\cos\theta_{in}$ to
$\sin(\theta_{in}+\theta^{\prime})\cos(\theta_{in}+\theta^{\prime})$, where
$\theta^{\prime}$ represents the angular shift. This shift in angle
$\theta^{\prime}$ increases with $m$ and equals $90\degree$ for
$m=6.9~{}\mu_{B}$ as shown in Fig. 2(d). The computed PHC discussed above is
normalized with respect to its maximum value. The anisotropy in the PHC is a
direct consequence of the anomalous velocity introduced by the anisotropic in-
plane components due to Eu moment and the external magnetic field.

To understand the variation in PHC with changing different parameters, we may
refer to the change in the magnitude of different components of the Berry
curvature ${\bf\Omega(k)}$. Fig. 3(a) shows the variation of the $x$-component
of the Berry curvature $(\Omega_{x})$ with the magnitude of Eu moment $m$. We
find that $\Omega_{x}$ increases linearly with the in-plane field and the
slope of this linear plot depends inversely on $m$. Therefore, the PHC also
varies inversely with $m$. Also, $\Omega_{x}$ increases with increasing the
azimuthal canting angle $\phi_{m}$ as shown in Fig. 3(b). However,
$\Omega_{y}$ shows an inverse relation, leading to no significant change in
the magnitude of the PHC on varying $\phi_{m}$.

The above results indicate that the anisotropy in the PHC with respect to the
applied in-plane field is due to the orientation of the Eu moment, while the
magnitude $m$ of the Eu moment is responsible for angular shift
$\theta^{\prime}$.

Figure 5: Planar Hall conductivity $\sigma_{xy}$, plotted in the plane of the
magnitude $B$ and the angle $\theta_{in}$ of the external in-plane magnetic
field, for the cases (a) $\phi_{m}\\!=\\!45\degree$, and (b)
$\phi_{m}\\!=\\!75\degree$. The initial canting angle of the Eu moment and the
critical magnetic field, considered in both cases, are
$\theta_{m0}\\!=\\!0\degree$, and $B_{c}\\!=\\!20~{}\mathrm{mT}$,
respectively. The magnetic moment is fixed to $m=6.9~{}\mu_{B}$ and the polar
canting angle follows the $\sin\theta_{m}=B/B_{c}$ relation.

###  III.2 Eu moment reorienting with in-plane field

Considering the general configuration for the magnetic moment of Eu ions ${\bf
m}\\!\equiv(m_{x},m_{y},m_{z})$ canted at a polar angle $\theta_{m}$ with the
$z$-axis and azimuthal angle $\phi_{m}$ with the $x$-axis on the $x$-$y$
plane, the energy-minimized configuration is given by (see Eq. (19) in
Appendix B) $\phi_{m}=\theta_{in}$ and

$\begin{split}\sin\theta_{m}&=\frac{StM_{s}B}{2Au^{2}m}=\frac{B}{B_{c}},~{}~{}~{}~{}B<B_{c}\\\
&=1,\hskip 71.13188ptB\geq B_{c}\end{split}$ (6)

where $t$ is the thickness of $\rm{EuS}$ film, $M_{s}$ is the saturation
magnetization, $S$ is the area of the two-dimensional system of topological
surface states, $u$ is the exchange coupling strength, and
$A=S(\sqrt{N_{c}}-\sqrt{|2N_{e}-N_{c}|})/(2\sqrt{\pi}\hbar|v_{f}|)$ with
$N_{c}$ and $N_{e}$ being the cutoff and total electron densities
respectively. The above relation indicates that the Eu moment turns in-plane
when the field amplitude $B$ exceeds the critical value $B_{c}$, which is
considered as a parameter in our discussion below.

Fig. 4(a) shows the variation in PHC with the above setting of the polar
canting angle $\theta_{m}$ and $\phi_{m}$ is fixed to $0\degree$. We observe
following three regimes of the in-plane field: (i)
$0\leq|B|<|B_{c^{\prime}}|$, (ii) $|B_{c^{\prime}}|\leq|B|\leq|B_{c}|$, and
(iii) $|B_{c}|<|B|$, where $B_{c^{\prime}}$ and $B_{c}$ are first and second
critical fields, respectively. The first critical field, $B_{c^{\prime}}$ is
the field value at which the Berry curvature becomes significant to contribute
to the Hall conductivity. The second critical field $B_{c}$, represents the
threshold field strength beyond which the Eu moment transitions to an in-plane
orientation.

The emergence of critical fields can be understood by looking at the behavior
of the components of the Berry curvature ${\bf\Omega(k)}$. For the lower
gapped surface band, the x-component of berry curvature $\Omega_{x}$ at
$(k_{x},k_{y})\equiv(0,0)$ is plotted as a function of the in-plane field in
Fig. 3(d)-(f). It is observed that $\Omega_{x}$ exhibits an increasing pattern
with $B$, attains a peak value, and subsequently diminishes to zero for $B\geq
B_{c}$. In the following discussion the second critical field $B_{c}$ is fixed
to $20~{}$mT, while other parameters are varied to observe the changes in
$B_{c^{\prime}}$, $\Delta B=|B_{c}-B_{c^{\prime}}|$, and the PHC.

Figure 6: Planar Hall conductivity $\sigma_{xy}$, plotted in the plane of the
magnitude $B$ and the angle $\theta_{in}$ of the external in-plane magnetic
field, for the cases (a) $\theta_{m0}\\!=\\!15\degree$, and (b)
$\theta_{m0}\\!=\\!45\degree$. The azimuthal canting angle and the critical
magnetic field, considered in both cases, are $\phi_{m}\\!=\\!\theta_{in}$,
and $B_{c}\\!=\\!0.3~{}\mathrm{T}$, respectively. The magnetic moment is fixed
to $m=6.9~{}\mu_{B}$ and the polar canting angle follows the
$\sin\theta_{m}=B/B_{c}$ relation.

Changing magnitude of Eu moment

The magnitude of the Eu moment is directly proportional to the thickness of
the interfaced thin film EuS, making it an important parameter for our study.
For lower values of the canted Eu moment, with a finite component along the
$z$-direction, the PHC gives a larger value as shown in Fig. 4. The increase
in the PHC with a decrease in magnetic moment is a direct consequence of a
smaller gap opening at lower field values, which leads to large Berry
curvature. Fig. 3(d) shows the variation of the $x$-component of the Berry
curvature with the in-plane magnetic field. The maximum value attained by
$\Omega_{x}$ increases with a decrease in the magnitude of the Eu moment. The
enhanced first critical field value can be understood from the increased field
range within which $\Omega_{x}$ remains large. Therefore, we can express the
dependence as

$|B_{c^{\prime}}|\propto|B_{c}|\propto{m}~{}~{}~{}{\rm and}~{}~{}~{}|\Delta
B|=|B_{c}-B_{c^{\prime}}|\propto\frac{1}{m}$

Changing azimuthal canting angle

As shown in Fig. 3(f) varying the azimuthal canting angle does not have a
significant effect on the critical fields but can affect the behavior of the
PHC significantly. To discuss this aspect, we consider two cases, as mentioned
below.

Case (a): $\phi_{m}=$ constant

As shown in Fig. 4, the PHC for $\phi_{m}=0\degree$ is dominant in the second
regime i.e., for $|B_{c^{\prime}}|\leq|B|\leq|B_{c}|$ and follows a
$\pm\cos\theta_{in}$ dependence for $\mp B$. Also, it has $2\pi$ periodicity
with antisymmetric behavior with respect to the external in-plane field. The
behavior of the PHC in this case is similar to the case where the Eu moment
does not change with the application of the in-plane field. This can be
understood by the fact that, when $B\to 0$, the polar canting angle is nearly
$0\degree$ and with increasing value of $|B|$, $\theta_{m}$ also increases. In
the second regime, $\theta_{m}$ nearly equals $60\degree$. Therefore, we get a
$-\cos\theta_{in}$ dependence. However, if $\phi_{m}$ is assumed to be
$90\degree$, $-\sin\theta_{in}$ dependence will be observed. The
$\cos\theta_{in}\tilde{v}_{y}$ and $\sin\theta_{in}\tilde{v}_{x}$ terms on the
RHS of Eq. (5) are responsible for the above behavior and it remains
consistent even at higher magnetic fields.

When $0\degree<\phi_{m}<90\degree$, both $\sin\theta_{in}\tilde{v}_{x}$ and
$\cos\theta_{in}\tilde{v}_{y}$ terms contribute and the PHC follows a linear
combination of $\sin\theta_{in}$ and $\cos\theta_{in}$. Fig. 5 depicts the
behavior of the PHC for $\phi_{m}=45\degree$ and $75\degree$. We find that the
PHC increases from zero to a finite value as the in-plane field magnitude
increases to the first critical field. However, it remains constant with the
in-plane field angle for $|B|<|B_{c^{\prime}}|$. In the
$|B_{c^{\prime}}|\leq|B|\leq|B_{c}|$ range, anisotropic PHC is observed with
$C\cos\theta_{in}+D\sin\theta_{in}$ dependence, where $C,D\in\mathbb{R}$.

Also, for $|B|>|B_{c}|$ the PHC is constant and non-zero, unlike in the
previous case; this is because $m_{y}\neq 0$ and the $x$ and the $y$
components of the Eu moment add to the Zeeman exchange due to the external in-
plane field.

Case (b): $\phi_{m}=\theta_{in}$

The azimuthal canting angle $\phi_{m}$ is now assumed to vary along with the
in-plane field angle $\theta_{in}$. The second critical field $B_{c}$ is taken
to be $0.3~{}$T for this case to understand the behavior at a higher magnetic
field. As shown in Fig. 6(a), the PHC is antisymmetric in the regime
$|B_{c^{\prime}}|\leq|B|\leq|B_{c}|$, with $\pi$ period, and follows a
$\sin(\theta_{in}+\delta)\cos(\theta_{in}+\delta)$ dependence, where $\delta$
represents the angular shift in the PHC. For smaller critical field values
$\delta\to 0\degree$, the antisymmetric behavior transitions to symmetric,
without any change in periodicity. However, to preserve the antisymmetric
behavior at smaller field values, the magnitude of Eu moment $m$ should be
lowered. For the regime $|B|>|B_{c}|$ in Fig. 6(a),
$\sin\theta_{in}\cos\theta_{in}$ dependence is observed, not affected by
further increase in the field value. This happens because $m_{z}$ becomes
zero, and $m_{x}$ and $m_{y}$ become constant for a field larger than the
second critical field. Since the magnitude of the Eu moment is much larger
than the in-plane field, we do not see any effect of increasing the in-plane
field on the PHC.

Table 1: Summary of the results discussed in Sec III. (A) Magnetic moment due
to EuS is fixed and does not change its orientation with the external in-plane
field. For different combinations of polar and azimuthal canting angles the
change in periodicity and symmetry of the PHC is mentioned. (B) Eu moment
reorients itself with the external field. The PHC is dominant in the range
$|B_{c^{\prime}}|<|B|<|B_{c}|$ of external magnetic field. The variation in
the nature of the PHC with changing the azimuthal canting angle $\phi_{m}$,
the initial canting angle $\theta_{m0}$, and the magnitude of the Eu moment
$m$ is summarized.

A. Fixed Orientation of Eu moment
---
Eu moment alignment | Hall conductivity
$m(\mu_{B})$ |  $\theta_{m}(\degree)$ |  $\phi_{m}(\degree)$ |  $\sigma_{xy}=f(\theta_{in})$ |  Period |  Symmetry with respect to negative field value
$0.0$ |  - |  - |  $\sin\theta_{in}\cos\theta_{in}$ |  $\pi$ |  Symmetric
$<6.9$ |  $0\degree$ |  - |  $\sin(\theta_{in}+\theta^{\prime})\cos(\theta_{in}+\theta^{\prime})$ |  $\pi$ |  Symmetric
$6.9$ |  $0\degree$ |  - |  $-\sin\theta_{in}\cos\theta_{in}~{}(\theta^{\prime}=90\degree)$ |  $\pi$ |  Symmetric
$90\degree$ |  $90\degree$ |  $\cos\theta_{in}$ |  $2\pi$ |  Antisymmetric
$90\degree$ |  $0\degree$ |  $\sin\theta_{in}$ |  $2\pi$ |  Antisymmetric
$0\degree<\theta_{m}<90\degree$ |  $0\degree$ |  $(A\sin\theta_{in}+B)\cos\theta_{in}$, where $A,B\in\mathbb{R}$ |  $2\pi$ |  Anisotropic

B. Polar canting angle of Eu moment is free to reorient according to
$\theta_{m}=\sin^{-1}(B/B_{c})+\theta_{m0}$
---
Varying azimuthal canting angle $\phi_{m}$ $m=6.9~{}\mu_{B}$ and $\theta_{m0}=0\degree$ |  $\phi_{m}(\degree)$ |  Hall conductivity in range $|B_{c^{\prime}}|\leq|B|\leq|B_{c}|$
$\sigma_{xy}=f(\theta_{in})$ |  Period |  Symmetry with respect to negative field value
$0\degree$ |  $-\cos\theta_{in}$ |  $2\pi$ |  Antisymmetric
$90\degree$ |  $-\sin\theta_{in}$ |  $2\pi$ |  Antisymmetric
$0\degree<\phi_{m}<90\degree$ |  $C\cos\theta_{in}+D\sin\theta_{in}$ where $C,D\in\mathbb{R}$ |  $2\pi$ |  Anisotropic
$\theta_{in}$ |  $-\sin(\theta_{in}+\delta)\cos(\theta_{in}+\delta)$ |  $\pi$ |  Antisymmetric (if $m<6.9~{}\mu_{B}$ or $B\sim 0.5~{}T$), symmetric otherwise ($\delta=0\degree$)

Varying initial canting angle $\theta_{m0}$, $m=6.9~{}\mu_{B}$ |  $\theta_{m0}(\degree)$ |  $\Delta B/B_{c}=|B_{c}-B_{c^{\prime}}|/B_{c}$ |  Varying magnitude of Eu moment $m~{}(\mu_{B})$, $\theta_{m0}=0\degree$ |  $m~{}(\mu_{B})$ |  $\Delta B/B_{c}=|B_{c}-B_{c^{\prime}}|/B_{c}$
---|---|---|---|---|---
$0\degree$ |  $0.15$ |  $6.9$ |  $0.15$
$15\degree$ |  $0.30$ |  $3.5$ |  $0.30$
$45\degree$ |  $0.57$ |  $1.8$ |  $0.50$

Results with an initial canting angle
In the above discussion, the initial orientation is chosen to be along the $z$
direction. Now, we consider the case when the Eu moment is canted initially at
a polar angle $\theta_{m0}$. In this case, the polar canting angle varies with
the in-plane field as $\theta_{m}=\sin^{-1}(B/B_{c})+\theta_{m0}$. Fig. 6
shows variation in the PHC with $\theta_{m0}$. We find that increasing
$\theta_{m0}$ (i) leads to a significant decrease in both the first and second
critical field values and (ii) results in an increase in the difference of the
critical field values. This behavior of the PHC with $\theta_{m0}$ can be
understood, again by looking at the change in Berry curvature, as shown in
Fig. 3(e). The maximum of the Berry curvature decreases and the field range
within which the Berry curvature remains large increases with an increase in
$\theta_{m0}$. In this case, we can express the dependence of the PHC on
$\theta_{m0}$ as

$|B_{c^{\prime}}|\propto|B_{c}|\propto\frac{1}{\theta_{m0}}~{}~{}~{}{\rm
and}~{}~{}~{}|\Delta B|\propto\theta_{m0}$

Therefore, with the variation in $m$, $\theta_{m}$, $\phi_{m}$, $\theta_{m0}$,
$B_{in}$, and $\theta_{in}$, the PHC may acquire symmetric, antisymmetric, or
anisotropic behavior with $\pi$ or $2\pi$ periodicity.

Summary of the results presented in Sec. III
We discussed the anomalous planar Hall effect for restricted and unrestricted
orientation of the magnetic moment due to thin film EuS in response to the
external in-plane magnetic field. The first case focuses on the behavior of
the planar Hall signal with different possible orientations of the easy axis
alignment as summarized in panel A of Table 1. Our computational results show
that the orientation and magnitude of the Eu moment play an important role in
determining the characteristics of the PHC at the interface. In the second
case, the Eu moment reorients itself with the in-plane field to minimize the
free energy of the system. Our analysis follows the energy-minimized solution
of canting in the presence of the external in-plane field. The polar canting
angle $\theta_{m}$ exhibits a monotonically-increasing behavior with respect
to changes in the in-plane magnetic field, and the azimuthal canting angle
$\phi_{m}$ is considered to be (i) fixed and (ii) changing with in-plane field
angle. The dependence of the polar canting angle on the in-plane field
introduces two critical fields. The first critical field is not only
influenced by the second critical field but also depends on factors like the
initial polar canting angle and the magnitude of the Eu moment. The second
critical field is determined by the thickness of the thin film EuS, the
magnitude of the Eu moment, and the initial polar canting angle. Planar Hall
conductivity is dominant for the in-plane field in the second regime,
$|B_{c^{\prime}}|\leq|B_{in}|\leq|B_{c}|$. The nature of the PHC in the second
regime is summarized in panel B of Table 1. The interplay between the in-plane
field and the Eu moment leads to distinct anisotropic behaviors with different
periodicity in various configurations.

##  IV Hall effect from magnetic skyrmions at the interface

Due to the presence of strong spin-orbit coupling in Bi2Se3 and broken
inversion symmetry at the interface, Dzyaloshinskii-Moriya (DM)-type
antisymmetric spin exchange interactions can appear on the EuS side of the
interface. This DM interaction is known to produce non-trivial magnetic
texture which can lead to topological Hall effect Nagaosa and Tokura (2013);
Heinze  _et al._ (2011); Yokouchi  _et al._ (2015); Yu  _et al._ (2010);
Mühlbauer  _et al._ (2009); Mohanta  _et al._ (2019). To explore the
possibility of the formation of non-trivial spin texture at the interface and
the occurrence of the topological Hall effect in our considered planar Hall
geometry, we perform Monte Carlo calculations (See Appendix F for details),
based on the following spin Hamiltonian

Figure 7: (a) Illustration of the model for skyrmion formation at the
$\rm{Bi}_{2}\rm{Se}_{3}/\rm{EuS}$ interface. The spin texture at the interface
is depicted for the cases (b) $B_{in}=0J$ and (c) $7J$ with $J=-1$, $D=1.4J$,
$m=5J$, $\theta_{m}=\phi_{m}=0\degree$, and $\theta_{in}=0\degree$. (d) Planar
topological Hall conductivity plotted for $\alpha=0.5J$ and $B_{in}$ in the
range $-2.5J$ to $2.5J$; other parameters are the same as above.

$\begin{split}H_{spin}=&-J\sum_{\langle ij\rangle}{\bf S}_{i}\cdot{\bf
S}_{j}-\sum_{\langle ij\rangle}{\bf D}_{ij}\cdot({\bf S}_{i}\times{\bf
S}_{j})\\\ &-\sum_{i}({\bf m}+{\bf B})\cdot{\bf S}_{i}\end{split}$ (7)

where $J$ is the strength of ferromagnetic Heisenberg exchange coupling, ${\bf
D}_{ij}=D(r_{i}-r_{j})\times\hat{z}/|r_{i}-r_{j}|$ is the DM vector acting
between $i$ and $j$ lattice sites in the $x$-$y$ plane with $D$ being the
strength of DM interaction, ${\bf B}$ is the external in-plane magnetic field,
and ${\bf m}$ is the magnetic moment due to EuS. Here, $i$ and $j$ represent
nearest neighbor lattice sites.

Topologically non-trivial spin textures, such as magnetic skyrmions, have been
both theoretically proposed and observed in experiments via topological Hall
measurements Hurst  _et al._ (2015); Divic  _et al._ (2022); Wu  _et al._
(2020); Zhang  _et al._ (2018); Li  _et al._ (2021). In our considered planar
Hall geometry, magnetic skyrmions can appear naturally when the Eu moment is
aligned out of the interface plane, as found in first-principle analysis of
the interface Kim  _et al._ (2017), or there is a suitable out-of-plane
contribution to the Zeeman exchange. These magnetic skyrmions can exchange-
couple to the gapped Dirac surface states on the surface of the topological
insulator and produce a planar topological Hall effect. This scenario is shown
schematically in Fig. 7(a). In the absence of the applied in-plane magnetic
field, a skyrmion crystal phase appears spontaneously in our Monte Carlo
calculations at low temperatures, within a range of values of the out-of-plane
Eu moment $m$. Fig. 7(b) shows the spin configuration on the EuS side of the
interface in the skyrmion crystal phase, obtained with parameters
$J\\!=\\!-1$, $D\\!=\\!1.4J$, $m\\!=\\!5J$, and $B_{in}\\!=\\!0$. In this
skyrmion crystal phase, when the in-plane magnetic field is turned on, the
skyrmions get deformed due to spins canted along the field direction, and
eventually beyond a critical field value an in-plane ferromagnetic phase is
established. The spin texture with a deformed skyrmion at an intermediate
field value is shown in Fig. 7(c).

To calculate the Hall conductivity arising due to the magnetic texture with a
non-zero scalar spin chirality, we consider that the spin texture is coupled
to the conduction electron on the surface of Bi2Se3. For this purpose, we
consider the below Hamiltonian,

$\begin{split}\mathcal{H}_{ex}=&-t\sum_{\langle
ij\rangle}c^{{\dagger}}_{i\sigma}c_{j\sigma}-J^{\prime}\sum_{i,\sigma,\sigma^{\prime}}({\bf
S}_{i}\cdot\bm{\sigma}_{\sigma\sigma^{\prime}})c^{{\dagger}}_{i\sigma}c_{i\sigma^{\prime}}\\\
&-\frac{i\alpha}{2a}\sum_{\langle
ij\rangle\sigma,\sigma^{\prime}}\hat{z}\cdot(\bm{\sigma}\times{\bf
d}_{ij})_{\sigma\sigma^{\prime}}c^{{\dagger}}_{i\sigma}c_{j\sigma^{\prime}},\end{split}$
(8)

where the first term corresponds to the tight-binding hopping term, the second
is the exchange coupling of the magnetic texture on the EuS side of the
interface with the itinerant electrons on the Bi2Se3 surface, and the third
term represents the Rashba hopping term with ${\bf d}_{ij}$ being the unit
vector between sites $i$ and $j$, $\alpha$ the rashba coupling strength, $a$
the lattice spacing, and
${\bm{\sigma}}\equiv(\sigma_{x},\sigma_{y},\sigma_{z})$ the Pauli matrices. By
diagonalizing the above Hamiltonian, from the eigenvalues and eigenvectors, we
calculate the transverse Hall conductivity using the Kubo formula, given by

$\sigma_{xy}=\frac{e^{2}}{h}\frac{2\pi}{N}\sum_{\epsilon_{m}\neq\epsilon_{n}}\frac{f_{m}-f_{n}}{(\epsilon_{m}-\epsilon_{n})^{2}+\eta^{2}}{\rm
Im}(\bra{m}\hat{j}_{x}\ket{n}\bra{n}\hat{j}_{y}\ket{m})$ (9)

where, $\hat{j_{x}}$ and $\hat{j_{y}}$ are the current operators in the $x$
and $y$\- directions, $f_{m/n}$ is the fermi-dirac distribution at temperature
$T$ and energy $\epsilon_{m/n}$, $\ket{m}$ is the $m$th eigenstate of
$\mathcal{H}_{ex}$, $N(L_{x}\times L_{y})$ is the lattice size, and $\eta$ is
the relaxation rate. For our calculations, we have used $\alpha=0.5$,
$T=0.001J$, $L_{x}=L_{y}=20$, and $\eta=0.1$. Fig. 7(d) shows the variation of
the topological Hall conductivity with increasing amplitude of the in-plane
magnetic field $B_{x}$, applied along the $x$ direction. The large value of
$\sigma_{xy}$ at $B_{x}\\!=\\!0$ decreases monotonically with increasing
$B_{x}$, and becomes zero at large fields.

##  V Conclusion

Planar Hall effect in topological insulators with a gapped surface band below
the Fermi level can originate from the Berry curvature-induced anomalous
velocity. On the other hand, planar topological Hall effect can arise at the
interface between a topological insulator and a magnetic insulator from a spin
texture with finite scalar spin chirality. In this paper, we explored both
anomalous and topological Hall effects in the planar Hall geometry based on
Bi2Se3/EuS interfaces. The canting of the Eu moment near the interface can
lead to an anisotropy in the Hall conductivity with respect to the applied in-
plane magnetic field. We also found that the preformed skyrmions can be
distorted in the presence of the in-plane magnetic field. When the Fermi
energy is just above the gapped Dirac surface band, the Berry curvature-driven
Hall effect dominates. However, if the Fermi energy is situated within the
bulk bands, the contribution from the non-trivial spin texture is expected to
be significant, as will be presented in another study Suri  _et al._ . The
anisotropy in the planar Hall conductivity, when the Fermi energy is within
the bulk energy gap, appears from the magnetic proximity effect of the canted
Eu moment and the in-plane applied field on the surface states of the
topological insulator. Our results can be useful to determine the nature of
the canting of the Eu moment using the anisotropy in planar Hall conductivity
in Bi2Se3/EuS interfaces, and possibly in other similar systems.

## Acknowledgement

NM acknowledges support of an initiation grant (No. IlTR/SRIC/2116/FIG) from
IIT Roorkee and SRG grant (No. SRG/2023/001188) from SERB. Calculations were
performed at the computing resources of PARAM Ganga at IIT Roorkee. JS was
supported by Ministry of Education, Government of India via a research
fellowship. KVR acknowledges support by Department of Atomic Energy,
Government of India, under Project Identification No. RTI 4007, ONRG Grant No.
N62909-23-1-2049, DST Grant No. DST/NM/TUE/QM-9/2019 (G), and SERB CRG Grant
No. CRG/2019/003810.

##  Appendix A Derivation of slab Hamiltonian

The effective low energy bulk Hamiltonian for $\mathrm{Bi}_{2}\mathrm{Se}_{3}$
in the basis $\ket{P1^{+}_{z},\uparrow}$, $\ket{P2^{-}_{z},\uparrow}$,
$\ket{P1^{+}_{z},\downarrow}$, $\ket{P2^{-}_{z},\downarrow}$ where $P1_{z}$
and $P2_{z}$ are $p_{z}$ orbitals of $\mathrm{Bi}$ and $\mathrm{Se}$, $\pm$
represents bonding and anti bonding states of orbitals is Mohanta  _et al._
(2017),

$H_{Bulk}(\textbf{k})=\begin{pmatrix}\tilde{\epsilon}_{+}(\textbf{k})&B_{1}k_{z}&0&A_{1}\tilde{k}_{-}\\\
B_{1}k_{z}&\tilde{\epsilon}_{-}(\textbf{k})&A_{1}\tilde{k}_{-}&0\\\
0&A_{1}\tilde{k}_{+}&\tilde{\epsilon}_{+}(\textbf{k})&-B_{1}k_{z}\\\
A_{1}\tilde{k}_{+}&0&-B_{1}k_{z}&\tilde{\epsilon}_{-}(\textbf{k})\\\
\end{pmatrix}$ (10)

Where,
$\tilde{\epsilon}_{\pm}(\textbf{k})=\tilde{\epsilon}_{0}(\textbf{k})\pm\tilde{M}(\textbf{k})$,
$\tilde{\epsilon}_{0}(\textbf{k})=\tilde{C}_{0}+\tilde{C}_{1}k_{z}^{2}+\tilde{C}_{2}(k_{x}^{2}+k_{y}^{2})$,
$\tilde{M}(\textbf{k})=\tilde{M}_{0}+\tilde{M}_{1}k_{z}^{2}+\tilde{M}_{2}(k_{x}^{2}+k_{y}^{2})$,
$\tilde{k}_{\pm}=k_{x}\pm ik_{y}$.

The lattice generalized Hamiltonian is obtained by substituting
$k_{i}a_{i}\simeq\sin k_{i}a_{i}$ and $(k_{i}a_{i})^{2}\simeq 2(1-\cos
k_{i}a_{i})$ in Hamiltonian given by 10.

$\displaystyle
H_{lattice}(\textbf{k})=\begin{pmatrix}\epsilon_{+}(\textbf{k})&B(k_{z})&0&A_{0}k_{-}\\\
B(k_{z})&\epsilon_{-}(\textbf{k})&A_{0}k_{-}&0\\\
0&A_{0}k_{+}&\epsilon_{+}(\textbf{k})&-B(k_{z})\\\
A_{0}k_{+}&0&-B(k_{z})&\epsilon_{-}(\textbf{k})\\\ \end{pmatrix}$ (11)

Where, $\epsilon_{\pm}(\textbf{k})=\epsilon_{0}(\textbf{k})\pm M(\textbf{k})$,
$\epsilon_{0}(\textbf{k})=C_{0}+2C_{1}(1-\cos k_{z}c)+2C_{2}(2-\cos
k_{x}a-\cos k_{y}a)$, $M(\textbf{k})=M_{0}+2M_{1}(1-\cos k_{z}c)+2M_{2}(2-\cos
k_{x}a-\cos k_{y}a)$, $k_{\pm}=\sin k_{x}a\pm i\sin k_{y}a$,
$B(k_{z})=B_{0}\sin k_{z}c$, $a$ and $c$ are lattice constants, and the
modified parameters are, $A_{0}=A_{1}/a$, $B_{0}=B_{1}/c$,
$C_{0}=\tilde{C}_{0}$, $C_{1}=\tilde{C}_{1}/c^{2}$,
$C_{2}=\tilde{C}_{2}/a^{2}$, $M_{0}=\tilde{M}_{0}$,
$M_{1}=\tilde{M}_{1}/c^{2}$, and $M_{2}=\tilde{M}_{2}/a^{2}$.

The slab Hamiltonian for $\mathrm{Bi}_{2}\mathrm{Se}_{3}$ can be obtained by
partial inverse Fourier transform of lattice generalized Hamiltonian.

$\begin{split}H_{slab}(\textbf{k})=\sum_{\textbf{k}}c_{\textbf{k}}^{{\dagger}}H_{lattice}(\textbf{k})c_{\textbf{k}}=\frac{1}{N^{2}}\sum_{\textbf{k}}\sum_{l_{z},l^{\prime}_{z}}e^{-ik_{z}l_{z}}\\\
c_{k_{x},k_{y},l_{z}}^{{\dagger}}H_{lattice}(\textbf{k})e^{ik_{z}l_{z}}c_{k_{x},k_{y},l_{z}}\end{split}$
(12)

$\begin{split}H_{slab}=\frac{1}{N}\sum_{\textbf{k}}\sum_{\alpha=1}^{4}(\sum_{l_{z}=1}^{N_{z}}c_{k\alpha
l_{z}}^{{\dagger}}H_{0}(\textbf{k})c_{k\alpha l_{z}}+\\\
\sum_{l_{z}=1}^{N_{z}-1}[c_{k\alpha l_{z}}^{{\dagger}}H_{1}c_{k\alpha
l_{z}+1}+H.c.])\end{split}$ (13)

The total slab Hamiltonian for $N$ quintuple layers in the basis
$\ket{P1^{+}_{z},\uparrow,l_{z}}$, $\ket{P2^{-}_{z},\uparrow,l_{z}}$,
$\ket{P1^{+}_{z},\downarrow,l_{z}}$, $\ket{P2^{-}_{z},\downarrow,l_{z}}$,
where $l_{z}$ is the layer index is given by:

$H(k)=\begin{pmatrix}H_{0}+H_{z}&H_{1}&&&&&\\\
H^{{\dagger}}_{1}&H_{0}&H_{1}&&&&\\\ &H^{{\dagger}}_{1}&H_{0}&H_{1}&&&\\\
&&.&&&&\\\ &&&.&&&\\\ &&&&.&&\\\ &&&&H^{{\dagger}}_{1}&H_{0}&H_{1}\\\
&&&&&H^{{\dagger}}_{1}&H_{0}\\\ \end{pmatrix}$ (14)

written in the basis $\ket{P1^{+}_{z},\uparrow,l_{z}}$,
$\ket{P2^{-}_{z},\uparrow,l_{z}}$, $\ket{P1^{+}_{z},\downarrow,l_{z}}$,
$\ket{P2^{-}_{z},\downarrow,l_{z}}$, where ${\bf
k}_{\parallel}\\!\equiv\\!(k_{x},k_{y})$, $l_{z}$ is the layer index which
varies between 0 and $N$. The Hamiltonians $H_{0}$, $H_{1}$ and $H_{z}$
describe, respectively, a single quintuple layer, coupling between two
neighboring layers, and Zeeman exchange coupling on the top surface layer due
to proximity effect of ferromagnetic EuS and due to the applied in-plane
magnetic field. These Hamiltonians are given by

$\begin{split}&H_{0}=\begin{pmatrix}\epsilon_{0}+M&0&0&A_{0}{\bf k_{-}}\\\
0&\epsilon_{0}-M&A_{0}{\bf k_{-}}&0\\\ 0&A_{0}{\bf k_{+}}&\epsilon_{0}+M&0\\\
A_{0}{\bf k_{+}}&0&0&\epsilon_{0}-M\\\ \end{pmatrix},\end{split}$ (15)

$H_{1}=\begin{pmatrix}-M_{1}-C_{1}&iB_{0}/2&0&0\\\ iB_{0}/2&M_{1}-C_{1}&0&0\\\
0&0&-M_{1}-C_{1}&-iB_{0}/2\\\ 0&0&-iB_{0}/2&M_{1}-C_{1}\\\ \end{pmatrix},$
(16)

$H_{z}=\begin{pmatrix}h_{z}&0&h_{x}-ih_{y}&0\\\ 0&h_{z}&0&h_{x}-ih_{y}\\\
h_{x}+ih_{y}&0&-h_{z}&0\\\ 0&h_{x}+ih_{y}&0&-h_{z}\\\ \end{pmatrix},$ (17)

where
$\epsilon_{0}(\mathbf{k}_{\parallel})\\!=\\!C_{0}+2C_{1}+2C_{2}(2-\cos(k_{x}a)-\cos(k_{y}a))$,
$M_{0}(\mathbf{k}_{\parallel})\\!=\\!M_{0}+2M_{1}+2M_{2}(2-\cos(k_{x}a)-\cos(k_{y}a))$,
$a$ is the lattice constant, ${\bf k_{\pm}}=\sin{(k_{x}a)}\pm
i\sin{(k_{y}a)}$, $h_{x}\\!=\\!m_{x}+g\mu_{B}B_{x}/2$,
$h_{y}\\!=\\!m_{y}+g\mu_{B}B_{y}/2$, and $h_{z}\\!=\\!m_{z}$.
$m_{x}\\!=\\!m\sin\theta_{m}\cos\phi_{m}$,
$m_{y}\\!=\\!m\sin\theta_{m}\sin\phi_{m}$, $m_{z}\\!=\\!m\cos\theta_{m}$,
$B_{x}\\!=\\!B\cos\theta_{in}$, and $B_{y}\\!=\\!B\sin\theta_{in}$, ${\bf
m}\\!\equiv(m_{x},m_{y},m_{z})$ represents the $\mathrm{EuS}$ magnetic moment
which is assumed to be canted at a polar angle $\theta_{m}$ with the $z$-axis
and azimuithal angle $\phi_{m}$ with the $x$-axis on the $x$-$y$ (interface)
plane, and $B$ is the magnitude of the external in-plane magnetic field
applied at an angle $\theta_{in}$ with the $x$-axis on the $x$-$y$ plane. The
parameter values used in this paper are
$A_{0}=0.8eV,~{}B_{0}=0.32eV,~{}C_{0}=-0.0083eV,~{}C_{1}=0.024eV,~{}C_{2}=1.77eV,~{}M_{0}=-0.28eV,~{}M_{1}=0.216eV,~{}M_{2}=2.6eV,~{}g=20,$
and $a=4.14{\AA}$ Mohanta  _et al._ (2017).

##  Appendix B Description of canting of the Eu moment

As discussed in the supplementary of Kim  _et al._ (2017), the total
Hamiltonian for $\mathrm{Bi_{2}}\mathrm{Se_{3}}/\rm{EuS}$ interfaces in the
presence of an external magnetic field can be written as

$H=\hbar v_{f}{\bm{\sigma}}\cdot({\bf
k}\times{\bf\hat{z}})+{\bm{\sigma}}\cdot(u{\bf m}-\mu_{B}{\bf B})-StM_{s}{\bf
m}\cdot{\bf B}+H_{\mathrm{EuS}}({\bf m})$ (18)

The effective total energy

$E_{tot}({\bf m})=\rm{constant}-A(u{\it m_{z}}-\mu_{B}{\it
B_{z}})^{2}-StM_{s}({\bf m\cdot B})$ (19)

Considering the general configuration for the magnetic moment due to EuS,
i.e., ${\bf m}\\!\equiv(m_{x},m_{y},m_{z})$ canted at a polar angle
$\theta_{m}$ with the $z$-axis and azimuthal angle $\phi_{m}$ with the
$x$-axis on the $x$-$y$ plane, the energy-minimized configuration with respect
to polar and azimuthal canting angles is given by $\phi_{m}=\theta_{in}$ and
$\theta_{m}=\sin^{-1}(StM_{s}B/2Au^{2}m)=B/B_{c},~{}{\rm for}~{}B<B_{c}~{}{\rm
and}~{}\theta_{m}=90\degree,~{}{\rm for}~{}B\geq B_{c}$. Here, $t$ is the
thickness of $\rm{EuS}$, $M_{s}$ is saturation magnetization and $B_{c}$ is
the critical field.

##  Appendix C Derivation of planar Hall conductivity

The semi-classical equations of motion for a Bloch electron in the presence of
electric and magnetic field are Ziman (2001); Xiao  _et al._ (2005); Dai  _et
al._ (2017); Min  _et al._ (2022); Morimoto  _et al._ (2016); Nandy  _et al._
(2018):

$\dot{{\bf r}}=\frac{1}{\hbar}\nabla_{k}\epsilon_{k}-\dot{{\bf
k}}\times{\bf\Omega({\bf k})}$ (20)

$\hbar\dot{{\bf k}}=-e{\bf E}-e\dot{{\bf r}}\times{\bf B}$ (21)

After substituting equation 21 in equation 20 and conducting algebraic
manipulation, the resulting equations are as follows:

$\dot{{\bf r}}=D({\bf B},{\bf\Omega(k))}~{}\Big{[}\tilde{{\bf
v}}_{k}+\frac{e}{\hbar}~{}{\bf
E}\times{\bf\Omega(k)}+\frac{e}{\hbar}~{}(\tilde{{\bf
v}}_{k}\cdot{\bf\Omega(k)}){\bf B}\Big{]}$ (22)

$\hbar\dot{{\bf k}}=-D({\bf B},{\bf\Omega(k)})~{}\Big{[}e{\bf E}+e(\tilde{{\bf
v}}_{k}\times{\bf B})+\frac{e^{2}}{\hbar}~{}({\bf E}\cdot{\bf
B})~{}{\bf\Omega(k)}\Big{]}$ (23)

Where, $D^{-1}({\bf B},{\bf\Omega(k)})=1+\frac{e}{\hbar}({\bf
B}\cdot{\bf\Omega(k)})$ and
$\tilde{{\bf
v}}_{k}=\frac{1}{\hbar}\nabla_{k}\tilde{\epsilon}_{k}=\frac{1}{\hbar}\nabla_{k}\epsilon_{k}-\frac{1}{\hbar}\nabla_{k}({\bf
m_{orb}}\cdot{\bf B})$ and $\tilde{\epsilon}_{k}=\epsilon_{k}-({\bf
m_{orb}}\cdot{\bf B})$ are modified velocity and energy respectively.

Here, ${\bf\Omega(k)}=(\Omega_{x},\Omega_{y},\Omega_{z})$ and ${\bf
m}_{orb}=(m_{orb}^{x},m_{orb}^{y},m_{orb}^{z})$ are the berry curvature and
orbital magnetic moment (OMM) respectively given by,

$\Omega^{n}_{xy}=-{\rm Im.}\sum_{m\neq
n}\frac{\braket{m}{\partial_{k_{x}}H}{n}\braket{n}{\partial_{k_{y}}H}{m}}{(E_{m}-E_{n})^{2}}$
(24) $m^{n}_{xy}=-\frac{e}{\hbar}{\rm Im.}\sum_{m\neq
n}\frac{\braket{m}{\partial_{k_{x}}H}{n}\braket{n}{\partial_{k_{y}}H}{m}}{E_{m}-E_{n}}$
(25)

For $n^{th}$ band, $\Omega^{n}_{i}=\epsilon_{ijk}\Omega^{n}_{jk}$ and
$m^{n}_{i}=\epsilon_{ijk}m^{n}_{jk}$.

Considering in-plane magnetic field with magnitude $B$ making an angle
$\theta_{in}$ with the electric field applied in the $x$-direction i.e., ${\bf
B}=B\cos\theta_{in}~{}\hat{{\bf x}}+B\sin\theta_{in}~{}\hat{{\bf y}}$ and
${\bf E}=E~{}\hat{{\bf x}}$. Now, substituting equations 22 and 23 in
Boltzmann transport equation,

$\dot{{\bf r}}\frac{\partial f({\bf r},{\bf k},t)}{\partial r}+\dot{{\bf
k}}\frac{\partial f({\bf r},{\bf k},t)}{\partial{\bf
k}}=-\frac{f-f_{0}}{\tau({\bf k})}=-\frac{f_{1}}{\tau({\bf k})}$ (26)

we get,

$\begin{split}-D({\bf B},{\bf\Omega(k)})~{}[\frac{e}{\hbar}{\bf
E}+\frac{e}{\hbar}(\tilde{{\bf v}}_{k}\times{\bf B})+\\\
\frac{e^{2}}{\hbar^{2}}({\bf E}.{\bf
B})~{}{\bf\Omega(k)}]~{}{\bf\nabla}_{k}\tilde{f}_{k}=-\frac{\tilde{f}_{k}-\tilde{f}_{eq}}{\tau({\bf
k})}\end{split}$ (27)

Keeping only linear terms in the applied field and using the ansatz
$\tilde{f_{k}}=\tilde{f_{eq}}+eED\tau~{}\big{[}\tilde{v}_{x}+\frac{eB\cos\theta_{in}}{\hbar}(\tilde{{\bf
v}}_{k}\cdot{\bf\Omega(k)})\big{]}~{}\partial_{\tilde{\epsilon}}\tilde{f}_{eq}$
we arrive at the following equation,

$\begin{split}De~{}\Big{[}Ev_{x}+\frac{e}{\hbar}({\bf E}\cdot{\bf
B})~{}(\tilde{{\bf
v}}_{k}\cdot{\bf\Omega(k)})~{}\partial_{\tilde{\epsilon}}\tilde{f}_{eq}+\\\
\frac{B}{\hbar}(-\tilde{v}_{z}\sin\theta_{in}\frac{\partial}{\partial
k_{x}}+\tilde{v}_{z}\cos\theta_{in}\frac{\partial}{\partial k_{y}}+\\\
(\tilde{v}_{x}\sin\theta_{in}-\tilde{v}_{y}\cos\theta_{in})\frac{\partial}{\partial
k_{z}})~{}\tilde{f}_{k}\Big{]}=-\frac{\tilde{f}_{eq}-\tilde{f}_{k}}{\tau({\bf
k})}\end{split}$ (28)

Where $\tilde{f}_{eq}$ is equilibrium Fermi-Dirac distribution with modified
energy. The charge density and current density modified to preserve
Liouville’s theorem are given by Xiao  _et al._ (2005):

$\rho=-\frac{e}{(2\pi)^{d}}\int d{\bf k}D^{-1}\tilde{f}_{k}$ (29)

$J=-\frac{e}{(2\pi)^{d}}\int d{\bf k}(D^{-1}\dot{{\bf
r}}+{\bf\nabla_{r}}\times{\bf m_{k}})\tilde{f_{k}}$ (30)

Where, $d$ is the dimension. Using the expression $J_{a}=\sigma_{ab}E_{b}$ and
using the ansatz for $\tilde{f}_{k}$, we get longitudinal and transverse PHC
expressions as,

$\begin{split}\sigma_{xx}&=-\frac{e^{2}}{(2\pi)^{d}}\int d{\bf
k}~{}D\tau\Big{[}\tilde{v}_{x}\\\
&+\frac{eB\cos\theta_{in}}{\hbar}(\bf{\tilde{v}_{k}}\cdot\Omega(\bf{k}))\Big{]}^{2}~{}\partial_{\tilde{\epsilon}}\tilde{f}_{eq},\end{split}$
(31)

$\begin{split}\sigma_{xy}=&-\frac{e^{2}}{(2\pi)^{d}}\int d{\bf
k}~{}D\tau~{}\Big{[}\tilde{v}_{y}\tilde{v}_{x}+\frac{eB\cos\theta_{in}~{}\tilde{v}_{y}}{\hbar}(\bf{\tilde{v}_{k}}\cdot\Omega(\bf{k}))\\\
&+\frac{eB\sin\theta_{in}~{}\tilde{v}_{x}}{\hbar}({\bf\tilde{v}_{k}}\cdot{\bf\Omega(k)})\\\
&+\frac{e^{2}B^{2}\cos\theta_{in}\sin\theta_{in}}{\hbar^{2}}(\bf{\tilde{v}_{k}}\cdot\bf{\Omega(k)})^{2}\Big{]}\partial_{\tilde{\epsilon}}\tilde{f}_{eq},\end{split}$
(32)

##  Appendix D Calculation of Berry curvature components $\Omega_{x}$ and
$\Omega_{y}$

The total slab Hamiltonian is $k_{z}$ independent, therefore $\Omega_{x}$ and
$\Omega_{y}$ are calculated using below approximation, treating the layer
index $l_{z}$ as the third variable:

$\frac{\partial}{\partial
k_{z}}\simeq\frac{\partial}{\partial\frac{1}{l_{z}}}=-l_{z}^{2}\frac{\partial}{\partial
l_{z}}$

$\Omega_{x}$ and $\Omega_{y}$ for perpendicular Zeeman coupling, gives Berry
curvature dipole centered at the gamma point in the $k$-space. The in-plane
components $\pm h_{y}(h_{x})$ shifts the centre of ${\bf\Omega}{\bf(k)}$ to
$\pm\delta{k_{y}}(-\delta{k_{x}})$.

##  Appendix E Monte Carlo annealing method

Low temperature spin configurations are obtained using Monte Carlo (MC)
annealing calculations on a square lattice of size 20$\times$20 with open
boundary conditions. The annealing procedure was started at a high temperature
$T\\!=\\!5J$ with an initial random spin configuration, and the system was
gradually cooled down from $T\\!=\\!5J$ to $T\\!=\\!0.001J$ in 100 steps. At
each temperature step, $10^{8}$ MC spin update steps are performed. At each
spin update step, a new orientation of the spin vector at a randomly-selected
site is chosen randomly within a small cone of angle $\Delta\theta=2^{\circ}$
around the initial spin direction. The new trial spin configuration is
accepted or rejected using the standard metropolis algorithm: accepted if the
energy difference $\Delta E$ between the two spin configurations is negative,
and accepted with the Boltzmann probability $e^{-\Delta E/(k_{B}T)}$ when
$\Delta E$ is positive; the latter is to incorporate the thermal spin
fluctuation which increases with temperature. Once the spin configuration is
stabilized (e.g. in the skyrmion crystal phase) at the lowest temperature
$T=0.001J$, a magnetic field along the in-plane $x$ direction is applied, and
$10^{10}$ MC spin updates are conducted at various values of the in-plane
field amplitude $B_{x}$ to obtain the spin configuration.

## References

  * Chang  _et al._ (2013) C.-Z. Chang, J. Zhang, X. Feng, J. Shen, Z. Zhang, M. Guo, K. Li, Y. Ou, P. Wei, L.-L. Wang, Z.-Q. Ji, Y. Feng, S. Ji, X. Chen, J. Jia, X. Dai, Z. Fang, S.-C. Zhang, K. He, Y. Wang, L. Lu, X.-C. Ma, and Q.-K. Xue, “Experimental observation of the quantum anomalous Hall effect in a magnetic topological insulator,” Science 340, 167 (2013).
  * Chang  _et al._ (2023) C.-Z. Chang, C.-X. Liu, and A. H. MacDonald, “Colloquium: Quantum anomalous Hall effect,” Rev. Mod. Phys. 95, 011002 (2023).
  * Sun  _et al._ (2019) H. Sun, B. Xia, Z. Chen, Y. Zhang, P. Liu, Q. Yao, H. Tang, Y. Zhao, H. Xu, and Q. Liu, “Rational design principles of the quantum anomalous Hall effect in superlatticelike magnetic topological insulators,” Phys. Rev. Lett. 123, 096401 (2019).
  * Pan  _et al._ (2020) L. Pan, X. Liu, Q. L. He, A. Stern, G. Yin, X. Che, Q. Shao, P. Zhang, P. Deng, C.-Y. Yang, B. Casas, E. S. Choi, J. Xia, X. Kou, and K. L. Wang, “Probing the low-temperature limit of the quantum anomalous Hall effect,” Sci. Adv. 6, eaaz3595 (2020).
  * Moon  _et al._ (2019) J. Moon, J. Kim, N. Koirala, M. Salehi, D. Vanderbilt, and S. Oh, “Ferromagnetic anomalous Hall effect in Cr-doped Bi2Se3 thin films via surface-state engineering,” Nano Lett. 19, 3409 (2019).
  * Morimoto  _et al._ (2015) T. Morimoto, A. Furusaki, and N. Nagaosa, “Topological magnetoelectric effects in thin films of topological insulators,” Phys. Rev. B 92, 085113 (2015).
  * Dziom  _et al._ (2017) V. Dziom, A. Shuvaev, A. Pimenov, G. V. Astakhov, C. Ames, K. Bendias, J. Böttcher, G. Tkachov, E. M. Hankiewicz, C. Brüne, H. Buhmann, and L. W. Molenkamp, “Observation of the universal magnetoelectric effect in a 3D topological insulator,” Nat. Commun. 8, 15197 (2017).
  * Wang  _et al._ (2012) M.-X. Wang, C. Liu, J.-P. Xu, F. Yang, L. Miao, M.-Y. Yao, C. L. Gao, C. Shen, X. Ma, X. Chen, Z.-A. Xu, Y. Liu, S.-C. Zhang, D. Qian, J.-F. Jia, and Q.-K. Xue, “The coexistence of superconductivity and topological order in the Bi2Se3 thin films,” Science 336, 52 (2012).
  * Shen  _et al._ (2017) J. Shen, W.-Y. He, N. F. Q. Yuan, Z. Huang, C.-w. Cho, S. H. Lee, Y. S. Hor, K. T. Law, and R. Lortz, “Nematic topological superconducting phase in Nb-doped Bi2Se3,” npj Quantum Mater. 2, 1 (2017).
  * Katmis  _et al._ (2016) F. Katmis, V. Lauter, F. S. Nogueira, B. A. Assaf, M. E. Jamer, P. Wei, B. Satpati, J. W. Freeland, I. Eremin, D. Heiman, P. Jarillo-Herrero, and J. S. Moodera, “A high-temperature ferromagnetic topological insulating phase by proximity coupling,” Nature 733, 513 (2016).
  * Kim  _et al._ (2017) J. Kim, K.-W. Kim, H. Wang, J. Sinova, and R. Wu, “Understanding the giant enhancement of exchange interaction in Bi2Se3/EuS heterostructures,” Phys. Rev. Lett. 119, 027201 (2017).
  * Mathimalar  _et al._ (2020) S. Mathimalar, S. Sasmal, A. Bhardwaj, S. Abhaya, R. Pothala, S. Chaudhary, B. Satpati, and K.V. Raman, “Signature of gate-controlled magnetism and localization effects at Bi2Se3-EuS interface,” npj Quantum Mater. 5, 64 (2020).
  * Assaf  _et al._ (2015) B. A. Assaf, F. Katmis, P. Wei, Cui-Zu Chang, B. Satpati, J. S. Moodera, and D. Heiman, “Inducing magnetism onto the surface of a topological crystalline insulator,” Phys. Rev. B 91, 195310 (2015).
  * Wei  _et al._ (2013) P. Wei, F. Katmis, B. A. Assaf, H. Steinberg, P. Jarillo-Herrero, D. Heiman, and J. S. Moodera, “Exchange-coupling-induced symmetry breaking in topological insulators,” Phys. Rev. Lett. 110, 186807 (2013).
  * Lee  _et al._ (2016) C. Lee, F. Katmis, P. Jarillo-Herrero, J.S. Moodera, and N. Gedik, “Direct measurement of proximity-induced magnetism at the interface between a topological insulator and a ferromagnet,” Nat. Commun. 7, 12014 (2016).
  * Taskin  _et al._ (2017) A.A. Taskin, H.F. Legg, F. Yang, S. Sasaki, Y. Kanai, K. Matsumoto, A. Rosch, and Y. Ando, “Planar Hall effect from the surface of topological insulators,” Nat. Commun. 8, 1340 (2017).
  * Nandy  _et al._ (2018) S. Nandy, A. Taraphder, and S. Tewari, “Berry phase theory of planar Hall effect in topological insulators,” Sci. Rep. 8, 14983 (2018).
  * Trushin  _et al._ (2009) M. Trushin, K. Výborný, P. Moraczewski, A. A. Kovalev, J. Schliemann, and T. Jungwirth, “Anisotropic magnetoresistance of spin-orbit coupled carriers scattered from polarized magnetic impurities,” Phys. Rev. B 80, 134405 (2009).
  * Chiba  _et al._ (2017) T. Chiba, S. Takahashi, and G. E. W. Bauer, “Magnetic-proximity-induced magnetoresistance on topological insulators,” Phys. Rev. B 95, 094428 (2017).
  * Zheng  _et al._ (2020) S.-H. Zheng, H.-J. Duan, J.-K. Wang, J.-Y. Li, M.-X. Deng, and R.-Q. Wang, “Origin of planar Hall effect on the surface of topological insulators: Tilt of Dirac cone by an in-plane magnetic field,” Phys. Rev. B 101, 041408 (2020).
  * Bhardwaj  _et al._ (2021) A. Bhardwaj, S. Prasad P., K. V. Raman, and D. Suri, “Observation of planar Hall effect in topological insulator—Bi2Te3,” Appl. Phys. Lett. 118, 241901 (2021).
  * Rakhmilevich  _et al._ (2018) D. Rakhmilevich, F. Wang, W. Zhao, M. H. W. Chan, J. S. Moodera, C. Liu, and C.-Z. Chang, “Unconventional planar Hall effect in exchange-coupled topological insulator–ferromagnetic insulator heterostructures,” Phys. Rev. B 98, 094404 (2018).
  * (23) D. Suri, S. Sasmal, A. Bhardwaj, J. Singh, S. Mundlia, A. Mishra, N. Mohanta, and K. V. Raman, “Emergence of planar topological Hall anisotropy in proximity-induced spin-canted state of Heisenberg ferromagnetic insulator EuS,” Communicated .
  * Zhang  _et al._ (2009) H. Zhang, C.-X. Liu, X.-L. Qi, X. Dai, Z. Fang, and S.-C. Zhang, “Topological insulators in $\mathrm{Bi}_{2}\mathrm{Se}_{3}$ , $\mathrm{Bi}_{2}\mathrm{Te}_{3}$, and $\mathrm{Sb}_{2}\mathrm{Te}_{3}$ with a single Dirac cone on the surface,” Nat. Phys. 5, 438 (2009).
  * Mohanta  _et al._ (2017) N. Mohanta, A. Kampf, and T. Kopp, “Emergent momentum-space skyrmion texture on the surface of topological insulators,” Sci. Rep. 7, 45664 (2017).
  * Liu  _et al._ (2010) C.-X. Liu, X.-L. Qi, H. Zhang, X. Dai, Z. Fang, and S.-C. Zhang, “Model Hamiltonian for topological insulators,” Phys. Rev. B 82, 045122 (2010).
  * Nagaosa and Tokura (2013) N. Nagaosa and Y. Tokura, “Topological properties and dynamics of magnetic skyrmions,” Nat. Nanotechnol. 8, 899 (2013).
  * Heinze  _et al._ (2011) S. Heinze, K. von Bergmann, M. Menzel, J. Brede, A. Kubetzka, R. Wiesendanger, G. Bihlmayer, and S. Blügel, “Spontaneous atomic-scale magnetic skyrmion lattice in two dimensions,” Nat. Phys. 7, 713 (2011).
  * Yokouchi  _et al._ (2015) T. Yokouchi, N. Kanazawa, A. Tsukazaki, Y. Kozuka, A. Kikkawa, Y. Taguchi, M. Kawasaki, M. Ichikawa, F. Kagawa, and Y. Tokura, “Formation of In-plane skyrmions in epitaxial MnSi thin films as revealed by planar Hall effect,” J. Phys. Soc. Jpn. 84, 104708 (2015).
  * Yu  _et al._ (2010) X. Z. Yu, Y. Onose, N. Kanazawa, J. H. Park, J. H. Han, Y. Matsui, N. Nagaosa, and Y. Tokura, “Real-space observation of a two-dimensional skyrmion crystal,” Nature 465, 901 (2010).
  * Mühlbauer  _et al._ (2009) S. Mühlbauer, B. Binz, F. Jonietz, C. Pfleiderer, A. Rosch, A. Neubauer, R. Georgii, and P. Böni, “Skyrmion lattice in a chiral magnet,” Science 323, 915 (2009).
  * Mohanta  _et al._ (2019) N. Mohanta, E. Dagotto, and S. Okamoto, “Topological Hall effect and emergent skyrmion crystal at manganite-iridate oxide interfaces,” Phys. Rev. B 100, 064429 (2019).
  * Hurst  _et al._ (2015) H. M. Hurst, D. K. Efimkin, J. Zang, and V. Galitski, “Charged skyrmions on the surface of a topological insulator,” Phys. Rev. B 91, 060401 (2015).
  * Divic  _et al._ (2022) S. Divic, H. Ling, T. Pereg-Barnea, and A. Paramekanti, “Magnetic skyrmion crystal at a topological insulator surface,” Phys. Rev. B 105, 035156 (2022).
  * Wu  _et al._ (2020) H. Wu, F. Groß, B. Dai, D. Lujan, S. A. Razavi, P. Zhang, Y. Liu, K. Sobotkiewich, J. Förster, M. Weigand, G. Schütz, X. Li, J. Gräfe, and K. L. Wang, “Ferrimagnetic skyrmions in topological insulator/ferrimagnet heterostructures,” Adv. Mater. 32 (2020).
  * Zhang  _et al._ (2018) S. Zhang, F. Kronast, G. van der Laan, and T. Hesjedal, “Real-space observation of skyrmionium in a ferromagnet-magnetic topological insulator heterostructure,” Nano Lett. 18, 1057 (2018).
  * Li  _et al._ (2021) P. Li, J. Ding, S.S.L. Zhang, J. Kally, O. G. Pillsbury, T.and Heinonen, G. Rimal, C. Bi, A. Demann, S. B. Field, W. Wang, J. Tang, J. S. Jiang, A. Hoffmann, N. Samarth, and M. Wu, “Topological Hall effect in a topological insulator interfaced with a magnetic insulator,” Nano Lett. 21, 84 (2021).
  * Ziman (2001) J. M. Ziman, “Electrons and Phonons: The Theory of Transport Phenomena in Solids,” Oxford University Press , ISBN 9780198507796 (2001).
  * Xiao  _et al._ (2005) D. Xiao, J. Shi, and Q. Niu, “Berry phase correction to electron density of states in solids,” Phys. Rev. Lett. 95, 137204 (2005).
  * Dai  _et al._ (2017) X. Dai, Z. Z. Du, and H.-Z. Lu, “Negative magnetoresistance without chiral anomaly in topological insulators,” Phys. Rev. Lett. 119, 166601 (2017).
  * Min  _et al._ (2022) W. Min, T. Daifeng, N. Yong, M. Shaopeng, G. Wenshuai, H. Yuyan, Z. Xiangde, Z. Jianhui, N. Wei, and T. Mingliang, “Novel $\pi/2$-periodic planar Hall effect due to orbital magnetic moments in $\mathrm{MnBi}_{2}\mathrm{Te}_{4}$,” Nano Lett. 22, 73 (2022).
  * Morimoto  _et al._ (2016) T. Morimoto, S. Zhong, J. Orenstein, and J. E. Moore, “Semiclassical theory of nonlinear magneto-optical responses with applications to topological Dirac/Weyl semimetals,” Phys. Rev. B 94, 245121 (2016).


* [3] D. L. Donoho, I. M. Johnstone, Ideal spatial adaptation by wavelet shrinkage, Biometrika 81 (3) (1994) 425–455.  doi:10.1093/biomet/81.3.425.
  * [4] T. Blumensath, M. E. Davies, Iterative hard thresholding for compressed sensing, Applied and Computational Harmonic Analysis 27 (3) (2009) 265–274.  doi:10.1016/j.acha.2009.04.002.
  * [5] B. Smith, Instantaneous companding of quantized signals, The Bell System Technical Journal 36 (3) (1957) 653–710.  doi:10.1002/j.1538-7305.1957.tb03858.x.
  * [6] P. Smaragdis, Dynamic range extension using interleaved gains, IEEE Audio, Speech, Language Process. 17 (5) (2009) 966–973.  doi:10.1109/tasl.2008.2012322.
URL http://dx.doi.org/10.1109/TASL.2008.2012322

  * [7] J. Abel, J. Smith, Restoring a clipped signal, in: IEEE Intl. Conf. on Acoustics, Speech and Sig. Proc. (ICASSP), 1991.  doi:10.1109/icassp.1991.150655.
URL http://dx.doi.org/10.1109/icassp.1991.150655

  * [8] J. Zhang, J. Hao, X. Zhao, S. Wang, L. Zhao, W. Wang, Z. Yao, Restoration of clipped seismic waveforms using projection onto convex sets method, Nature Sci. Rep. 6 (1).  doi:10.1038/srep39056.
URL https://doi.org/10.1038/srep39056

  * [9] A. Adler, V. Emiya, M. G. Jafari, M. Elad, R. Gribonval, M. D. Plumbley, Audio inpainting, IEEE Trans. Acoust., Speech, Signal Process. 20 (3) (2012) 922–932.  doi:10.1109/tasl.2011.2168211.
URL https://doi.org/10.1109/tasl.2011.2168211

  * [10] P. E. Debevec, J. Malik, Recovering high dynamic range radiance maps from photographs, in: ACM Siggraph, 1997.  doi:10.1145/258734.258884.
URL https://doi.org/10.1145/258734.258884

  * [11] A. Bhandari, F. Krahmer, R. Raskar, On unlimited sampling, in: Intl. Conf. on Sampling Theory and Applications (SampTA), 2017, pp. 31–35.  doi:10.1109/sampta.2017.8024471.
URL https://doi.org/10.1109/sampta.2017.8024471

  * [12] A. Bhandari, F. Krahmer, R. Raskar, On unlimited sampling and reconstruction, IEEE Trans. Sig. Proc. 69 (2020) 3827–3839.  doi:10.1109/tsp.2020.3041955.
  * [13] A. Bhandari, F. Krahmer, T. Poskitt, Unlimited sampling from theory to practice: Fourier-Prony recovery and prototype ADC, IEEE Trans. Sig. Proc. (2021) 1131–1141doi:10.1109/TSP.2021.3113497.
  * [14] D. Florescu, F. Krahmer, A. Bhandari, The surprising benefits of hysteresis in unlimited sampling: Theory, algorithms and experiments, IEEE Trans. Sig. Proc. 70 (2022) 616–630.  doi:10.1109/tsp.2022.3142507.
  * [15] A. Bhandari, F. Krahmer, R. Raskar, Unlimited sampling of sparse signals, in: IEEE Intl. Conf. on Acoustics, Speech and Signal Processing (ICASSP), 2018, pp. 4569–4573.  doi:10.1109/icassp.2018.8462231.
  * [16] A. Bhandari, F. Krahmer, R. Raskar, Unlimited sampling of sparse sinusoidal mixtures, in: IEEE Intl. Sym. on Information Theory (ISIT), 2018, pp. 336–340.  doi:10.1109/isit.2018.8437122.
  * [17] A. Bhandari, F. Krahmer, On identifiability in unlimited sampling, in: Intl. Conf. on Sampling Theory and Applications (SampTA), 2019.
  * [18] A. Bhandari, F. Krahmer, R. Raskar, Methods and apparatus for modulo sampling and recovery (May 2020).
  * [19] A. Bhandari, F. Krahmer, HDR imaging from quantization noise, in: IEEE Intl. Conf. on Image Processing (ICIP), 2020, pp. 101–105.  doi:10.1109/icip40778.2020.9190872.
  * [20] V. Shah, C. Hegde, Signal reconstruction from modulo observations, in: 2019 IEEE Global Conference on Signal and Information Processing (GlobalSIP), IEEE, 2019, pp. 1–5.
  * [21] O. Musa, P. Jung, N. Goertz, Generalized approximate message passing for unlimited sampling of sparse signals, in: IEEE Global Conf. on Signal and Information Proc., 2018, pp. 336–340.
  * [22] S. Rudresh, A. Adiga, B. A. Shenoy, C. S. Seelamantula, Wavelet-based reconstruction for unlimited sampling, in: IEEE Intl. Conf. on Acoustics, Speech and Signal Processing (ICASSP), 2018, pp. 4584–4588.
  * [23] A. Bhandari, Back in the US-SR: Unlimited sampling and sparse super-resolution with its hardware validation, IEEE Signal Process. Lett. 29 (2022) 1047 – 1051.
  * [24] D. Florescu, F. Krahmer, A. Bhandari, Event-driven modulo sampling, in: IEEE Intl. Conf. on Acoustics, Speech and Signal Processing (ICASSP), 2021, pp. 5435–5439.  doi:10.1109/icassp39728.2021.9414152.
  * [25] D. Florescu, A. Bhandari, Modulo event-driven sampling: System identification and hardware experiments, in: IEEE Intl. Conf. on Acoustics, Speech and Signal Processing (ICASSP), 2022, pp. 5747–5751.
  * [26] M. Beckmann, A. Bhandari, F. Krahmer, The modulo Radon transform: Theory, algorithms, and applications, SIAM Journal on Imaging Sciences 15 (2) (2022) 455–490.  doi:10.1137/21m1424615.
  * [27] S. Fernandez-Menduina, F. Krahmer, G. Leus, A. Bhandari, Computational array signal processing via modulo non-linearities, IEEE Trans. Sig. Proc. 70 (2022) 2168–2179.
  * [28] V. Bouis, F. Krahmer, A. Bhandari, Multidimensional unlimited sampling: A geometrical perspective, in: European Sig. Proc. Conf. (EUSIPCO), 2020, pp. 2314–2318.
  * [29] X. Zhao, J. Chen, M. Karczewicz, A. Said, V. Seregin, Joint separable and non-separable transforms for next-generation video coding, IEEE Trans. Image Proc. 27 (5) (2018) 2514–2525.
  * [30] A. Cohen, I. Daubechies, Non-separable bidimensional wavelet bases, Revista Matematica Iberoamericana 9 (1) (1993) 51–137.
  * [31] H. Tyagi, Error analysis for denoising smooth modulo signals on a graph, Applied and Computational Harmonic Analysis 57 (2022) 151–184.
  * [32] M. Fanuel, H. Tyagi, Denoising modulo samples: k-nn regression and tightness of sdp relaxation, Information and Inference: A Journal of the IMA 11 (2) (2022) 637–677.
  * [33] D. Florescu, A. Bhandari, Unlimited sampling with local averages, in: IEEE Intl. Conf. on Acoustics, Speech and Signal Processing (ICASSP), 2022, pp. 5742–5746.
  * [34] D. Florescu, A. Bhandari, Unlimited sampling via generalized thresholding, in: IEEE Intl. Symp. on Information Theory (ISIT), 2022, pp. 1606–1611.
  * [35] D. Florescu, F. Krahmer, A. Bhandari, Unlimited sampling with hysteresis, in: 55th Asilomar Conf. on Signals, Systems, and Computers, 2021, pp. 831–835.  doi:10.1109/ieeeconf53345.2021.9723306.
  * [36] E. Viscito, J. P. Allebach, The analysis and design of multidimensional FIR perfect reconstruction filter banks for arbitrary sampling lattices, IEEE Trans. Circuits Syst. 38 (1) (1991) 29–41.
  * [37] Y. M. Lu, M. N. Do, R. S. Laugesen, A computable Fourier condition generating alias-free sampling lattices, IEEE Trans. Sig. Proc. 57 (5) (2009) 1768–1782.
  * [38] M. Chiani, D. Dardari, M. K. Simon, New exponential bounds and approximations for the computation of error probability in fading channels, IEEE Trans. Wireless Commun. 2 (4) (2003) 840–845.
  * [39] J. Wozencraft, I. M. Jacobs, Principles of communication engineering., Wiley, 1965\.
  * [40] G. Gallego, T. Delbrück, G. Orchard, C. Bartolozzi, B. Taba, A. Censi, S. Leutenegger, A. J. Davison, J. Conradt, K. Daniilidis, et al., Event-based vision: A survey, IEEE Trans. Pattern Anal. Mach. Intell. 44 (1) (2020) 154–180.
  * [41] O. Graf, A. Bhandari, F. Krahmer, One-bit unlimited sampling, in: IEEE Intl. Conf. on Acoustics, Speech and Signal Processing (ICASSP), 2019, pp. 5102–5106.


# Solving combinatorial optimization problems through stochastic Landau-
Lifshitz-Gilbert dynamical systems

Dairong Chen  Center for Quantum Phenomena, Department of Physics, New York
University, New York, NY 10003 USA  Andrew D. Kent  Center for Quantum
Phenomena, Department of Physics, New York University, New York, NY 10003 USA
Dries Sels  Center for Quantum Phenomena, Department of Physics, New York
University, New York, NY 10003 USA  Center for Computational Quantum Physics,
Flatiron Institute, New York, NY 10010 USA  Flaviano Morone  Center for
Quantum Phenomena, Department of Physics, New York University, New York, NY
10003 USA

###### Abstract

We present a method to approximately solve general instances of combinatorial
optimization problems using the physical dynamics of 3d rotors obeying Landau-
Lifshitz-Gilbert dynamics. Conventional techniques to solve discrete
optimization problems that use simple continuous relaxation of the objective
function followed by gradient descent minimization are inherently unable to
avoid local optima, thus producing poor-quality solutions. Our method
considers the physical dynamics of macrospins capable of escaping from local
minima, thus facilitating the discovery of high-quality, nearly optimal
solutions, as supported by extensive numerical simulations on a prototypical
quadratic unconstrained binary optimization (QUBO) problem. Our method
produces solutions that compare favorably with those obtained using state-of-
the-art minimization algorithms (such as simulated annealing) while offering
the advantage of being physically realizable by means of arrays of stochastic
magnetic tunnel junction devices.

_Introduction._ A compelling idea linking physics and discrete mathematics is
to leverage the intrinsic dynamics of real-world physical systems to solve
mathematical optimization problem. The idea relies on correspondence between
the variables of the mathematical problem and the degrees of freedom of an
analogue physical system, where the role of the objective function to be
optimized is played by the energy function of the interacting physical system
[1, 2, 3, 4]. Once this correspondence has been successfully established, all
one has to do is to cool the system down to zero temperature and wait for it
to settle in its ground state, i.e., the state with minimal energy, which, in
turn, represents the solution to the original optimization problem [5]. In
many cases of practical interest, this is easier said than done in that during
the cooling process the physical system can get stuck in long-lived metastable
states, possibly without ever reaching the ground state on human-life
timescales. In such situations, the choice of the analogue physical system
becomes the crucial point of the entire methodology. Physical systems that are
less prone to getting trapped in the local minima of the energy function will
solve the optimization problem better than those that get easily stuck.

Many interesting optimization problems are discrete in nature, meaning that
the variables involved are represented by integers, typically taking up only
two values, say $0$ and $1$. A large body of work has thus exploited the idea
of using Ising-type physical systems to effectively map integer variables onto
binary spins and then use Monte Carlo methods or alike to simulate the
physical dynamics of these analogue Ising machines to find their ground state.
These techniques have been studied extensively in the literature (see, for
example, references [1, 3, 2, 6]).

In this article, we present a new energy-minimization framework to solve
combinatorial optimization problems based on the physical dynamics of real
macrospins described by the Landau-Lifshitz-Gilbert (LLG) equations. The idea
is to use 3-dimensional rotors as the analogue physical dynamical system,
rather than up-down Ising spins. The rationale behind our approach is the
physical gambit that trajectories away from the z-axis, in particular, those
occurring halfway between the north and south poles of the unit sphere will
help the system to escape from local solutions, thus facilitating the
discovery of the true ground-state.

We explain our method in two steps: first, we introduce the stochastic LLG
dynamical system and discuss its numerical implementation; second, we explain
how to apply it to a specific QUBO problem, called the Sherrington–Kirkpatrick
model in statistical physics [7, 8], which can be exactly solved in the
thermodynamic limit [9], and show that our method extrapolates to this exact
value while performing systematically better than discrete methods based on
Glauber dynamics of Ising spins [10].

Figure 1: Schematic description of our method. First, we encode the solution
of the optimization problem into the ground state of the spin system’s
Hamiltonian. The parameters $J_{ij}$ model the interactions between pair of
macrospins. Here we only show couplings between neighbouring spins for
illustration purposes. In general, the connections are not restricted to the
nearest neighbours. Anisotropy is set along the $z$ axis of each spin via a
magnetic anisotropy field $H_{Ai}$ favoring the north-south direction. After
initializing the macrospins in random directions on the unit sphere, the
system is cooled from high to low temperature. At high temperatures, each
macrospin experiences large thermal fluctuations with trajectories densely
covering the unit sphere. As the system cools to lower temperatures, the
fluctuations are reduced and the spins trajectories tend to confine into small
regions around the poles of the unit sphere. The dynamical trajectories of a
single macrospin during the annealing process is depicted in green on the
left-hand side of the illustration. Eventually, the system reaches the end of
the annealing schedule and spins stabilize around the north and south poles.
The spins’ final states are binarized based on their up or down direction to
get a direct readout of the solution to the original problem.

_Stochastic-LLG Dynamics._ We consider a system consisting of $n$ single-
domain ferromagnets, and apply the macrospin approximation to each domain. The
magnetization of the spins in the $i^{\mathrm{th}}$ domain are averaged to a
single giant magnetic moment $\vec{m}_{i}=m_{i}\hat{m}_{i}$, with magnitude
$m_{i}>0$ and direction $\hat{m}_{i}$, as illustrated in Fig. 1. The
micromagnetic energy $E$ of our $n$ macrospins system has three terms: the
energy from the exchange interaction between each moments $E_{\text{ex}}$, the
energy from the magnetic anisotropy determined by the crystallographic
structure of the domain, $E_{\text{anis}}$, and the Zeeman energy from an
external magnetic field applied to the system, $E_{\text{Z}}$. Summarizing:

$\displaystyle E=E_{\text{ex}}+E_{\text{anis}}+E_{\text{Z}}.$ (1)

The expression for each term is given by:

$\displaystyle E_{\text{ex}}$ $\displaystyle=-\mu_{0}\sum_{\langle
i,j\rangle}^{n}2J_{ij}H_{\text{ex}}\frac{\vec{m}_{i}\cdot\vec{m}_{j}}{m_{i}+m_{j}},$
(2) $\displaystyle E_{\text{anis}}$
$\displaystyle=-\frac{1}{2}\mu_{0}\sum_{i=1}^{n}H_{Ai}m_{i}(\vec{e}_{hi}\cdot\hat{m}_{i})^{2},$
(3) $\displaystyle E_{\text{Z}}$
$\displaystyle=-\mu_{0}\sum_{i=1}^{n}\vec{m}_{i}\cdot\vec{H}_{\text{ext},i},$
(4)

where $\mu_{0}$ is the vacuum permeability, $H_{Ai}$ characterizes the
magnetic anisotropy of moment $i$, and $\vec{e}_{hi}$ is the unit vector
characterizing the preferred crystallographic direction of moment $i$ (which
later we take to be along the z-axis, see Fig. 2). The summation over $\langle
i,j\rangle$ is over pairs of connected spins. Vector $\vec{H}_{\text{ext},i}$
is an external local magnetic field acting on the $i^{\mathrm{th}}$ magnetic
moment. $H_{\text{ex}}$ has units of field and it characterizes the general
scale of the exchange coupling strength. Coupling terms $J_{ij}$ are
dimensionless and characterize the exchange interaction between magnetic
moments.

Each magnetic moment $\vec{m}_{i}$ evolves according to the Landau-Lifshitz-
Gilbert (LLG) equation [11]:

$\displaystyle\frac{d\vec{m}_{i}}{dt}=\gamma_{i}^{\prime}\vec{m}_{i}\times\vec{H}_{\text{eff},i}-\frac{\alpha_{i}}{m_{i}}\vec{m}_{i}\times\frac{d\vec{m}_{i}}{dt},$
(5)

where $\gamma_{i}^{\prime}=\mu_{0}\gamma_{i}$, with $\gamma_{i}$ the
gyromagnetic ratio of the $i^{\mathrm{th}}$ moment, and $\alpha_{i}$ is the
damping constant. Vector $\vec{H}_{\text{eff},i}$ is the effective field
experienced by the $i^{\mathrm{th}}$ magnetic moment, given by

$\displaystyle\vec{H}_{\text{eff},i}=-\frac{1}{\mu_{0}}\frac{\partial
E}{\partial\vec{m}_{i}}.$ (6)

To incorporate thermal fluctuations into the dynamics we add a thermal field
$\vec{H}_{\text{th},i}$ to $\vec{H}_{\text{eff},i}$ by the substitution

$\displaystyle\vec{H}_{\text{eff},i}\rightarrow\vec{H}_{\text{eff},i}+\vec{H}_{\text{th},i}.$
(7)

The thermal fields are delta-correlated zero mean random variables in time,
with a correlation given by [12, 13]:

$\displaystyle\langle\vec{H}_{\text{th},i,\mu}(t)\vec{H}_{\text{th},i,\nu}(t^{\prime})\rangle=C_{i}\delta_{\mu\nu}\delta(t-t^{\prime}),$
(8) $\displaystyle
C_{i}=\frac{2\alpha_{i}k_{B}T}{\mu_{0}^{2}m_{i}\gamma_{i}(1+\alpha_{i}^{2})},$
(9)

where $\mu,\nu=x,y,z$, $k_{B}$ is the Boltzmann constant, and $T$ is the
temperature. Integrating Eq. (9) over a small time interval $\Delta t$, we
obtain the standard deviation of $\vec{H}_{\text{th},i}$ as

$\displaystyle\sigma_{i}=\frac{1}{\mu_{0}}\sqrt{\frac{2\alpha_{i}k_{B}T}{m_{i}\gamma_{i}(1+\alpha_{i}^{2})\Delta
t}}.$ (10)

In practice, for each component of $\vec{H}_{\text{th},i}$, we draw a value
from a Gaussian distribution with mean $\mu=0$ and standard deviation given by
$\sigma_{i}$ in Eq. 10. Finally, to integrate the stochastic differential
equation we apply Heun’s integration scheme [14].

In our simulation, we assume that each macrospin $i$ has the same
crystallographic structure. All parameters related to the material properties
are the same across all $n$ domains, therefore the index $i$ of the following
quantities are omitted: $H_{Ai}=H_{A}$, $\alpha_{i}=\alpha$,
$\gamma_{i}=\gamma$, $\gamma_{i}^{\prime}=\gamma^{\prime}$, and the
crystallographic direction for each spin is set along the z-axis,
$\vec{e}_{hi}=\hat{z}$. We also assume each magnetic moment has the same
magnitude, $m_{i}=m$ for all $n$ macrospins. This simplification models an
array of perfectly identical macrospins. The simulation can be further
simplified to only consider the evolution of the unit magnetic moment
$\hat{m}_{i}$. Given the assumptions above, we can think of our model as a
system of $n$ identical dynamical bits. At the end of the simulation, a
projection on the z-axis, given by $s_{i}\equiv{\rm sign}(\hat{m}_{i,z})$,
maps the macrospin onto a binary Ising-like spin $s_{i}=\pm 1$, as illustrated
in Fig. 2.

Figure 2: Illustration of Ising spin and macrospin. Ising spins can only take
two values, $s_{i}=\pm 1$. Macrospins, on the other hand, have extra degrees
of freedom when rotating on the unit sphere.

The values of the parameters used in our simulation are listed in Table 1. The
values we chose here are within those of the fluctuating, the so-called free
layer, of magnetic tunnel junction nanopillars [15].

Damping Constant | $\alpha$ | 0.01
---|---|---
Magnetic Moment | m |  $8\times 10^{-19}$J/T
Ambient Temperature | $T_{\text{amb}}$ | 300K
Exchange Field | $H_{\text{ex}}$ |  $8.24\times 10^{3}$ A/m
Anisotropy Field | $H_{A}$ |  $8.24\times 10^{3}$ A/m
Simulation Time Step | $\Delta t$ |  $0.01$ ns
Iterations at Each Temperature |  | $2\times 10^{6}$
Table 1: Values of the parameters chosen for simulation.

With these parameters, we implement an annealing schedule to mimic the
physical cooling of the system from ambient temperature down to the low-
temperature regime. Ideally, the annealing schedule should be performed slowly
from $T_{\text{amb}}=300$ to $0$ K. In practice, we choose an annealing
schedule with $n_{T}=30$ intermediate temperatures
$T_{1}>T_{2}>\cdots>T_{n_{T}}$ drawn from a distribution $P(T)\propto
1/T^{\delta}$, with $\delta=1.1$, in order to have more values in the low
temperature regime. Having explained how to simulate the general stochastic-
LLG dynamical system, we move next to discuss its application to the SK model.

_Results._ Here we show the simulation results for the SK model using the LLG
dynamics described above. We also run a Glauber dynamics for comparison [10].
The Hamiltonian of the SK model for a $n$-spin system is given by:

$\displaystyle H_{SK}(s)=-\sum_{\langle
i,j\rangle}^{n}\frac{J_{ij}}{\sqrt{n}}\ s_{i}s_{j},$ (11)

where $s_{i}=\pm 1$ are Ising spins, $J_{ij}$ is the exchange interaction
between the $i^{\mathrm{th}}$ and $j^{\mathrm{th}}$ spins. Couplings $J_{ij}$
are independent and identically distributed random variables drawn from a
Gaussian distribution with zero mean and variance $\sigma^{2}=1$. The
summation only counts a pair of $\langle i,j\rangle$ once. The ground state
energy $E_{\text{min}}\equiv\min_{s}H_{SK}(s)$ has been calculated
analytically in the limit of infinite system size $n\rightarrow\infty$ [9] and
is given by

$\lim_{n\to\infty}\frac{E_{\text{min}}}{n}=-0.76321...\ .$ (12)

We can map the Hamiltonian $H_{SK}$ onto our LLG model by directly mapping the
$J_{ij}$ terms in Eq. (11) to the $J_{ij}$ terms in Eq. (2). The
$\vec{H}_{\text{ext},i}$ term in Eq. (4) is set to zero since there is no
external field applied in the SK model. With the anisotropy term favoring each
spin to align along the z-axis, $\vec{e}_{hi}=\hat{z}$, the ground state of
the SK model is mapped onto the ground state of the LLG system once we
binarize the macrospins using their z-projections as illustrated in Fig. 2.
Thus, the analogue LLG Hamiltonian reads

$\displaystyle H_{SK}^{LLG}$
$\displaystyle=-2\mu_{0}H_{\text{ex}}\sum_{\langle
i,j\rangle}^{n}J_{ij}\frac{\vec{m}_{i}\cdot\vec{m}_{j}}{m_{i}+m_{j}}$ (13)
$\displaystyle-\frac{1}{2}\mu_{0}\sum_{i}^{n}H_{Ai}m_{i}(\hat{z}\cdot\hat{m}_{i})^{2}\
.$

We solve the SK model for several system sizes, ranging from $n=40$ to $2000$,
using both Glauber and LLG dynamics and record the value of the final energy
$E$ obtained with both methods. As shown in Fig. 3, our LLG method performs
better, i.e., lower energies are consistently found for any value of $n$.

Figure 3: Histograms of the final energies $\langle E\rangle/n$ obtained from
LLG dynamics (orange bars) and Glauber dynamics (blue bars) for $n=1200,1600$,
and $2000$. The colored vertical dashed lines represent the average values of
the corresponding distributions. To make the histograms we used 300 samples
for $n=1200$, $103$ for $n=1600$, and $112$ for $n=2000$.  Figure 4: Average
ground state energy per spin $\langle E\rangle/n$ as a function of $n^{-2/3}$
[16] for Glauber (blue line) and LLG (orange line) dynamics. The horizontal
dashed line represents the asymptotic $n\to\infty$ exact ground state energy
given by Eq. (12). Data points from the LLG dynamics plausibly extrapolate to
the true optimal value. Error bars indicate the standard error of the mean
energy per spin.

In Fig. 4 we plot the average ground state energy per spin, $\langle
E\rangle/n$, obtained from Glauber and LLG as a function of $n^{-2/3}$ [16].
We observe that the LLG dynamics results extrapolate to $\langle
E\rangle/n=-0.762(1)$ as $n\to\infty$, which is very close to the exact
optimal value Eq. (12), while the Glauber dynamics results intercept the $y$
axis at a higher sub-optimal asymptotic value $\langle E\rangle/n=-0.755(2)$.

Recently it is has been shown, assuming a widely believed conjecture about the
SK model, that there is a polynomial time algorithm for finding (with high
probability) a string whose energy is between the lowest energy and
$(1-\epsilon)$ of the lowest energy of the problem [17]. Our results suggest
we can achieve similar performance.

As a point of reference, if the algorithm would get stuck at the critical
point of the SK model (at $T=1$) the typical energy per spin would be $\langle
E\rangle/n=-0.5$. However, even greedy spin flip dynamics, in which one only
accepts flips if they lower the energy, results in a energy of around $\langle
E\rangle/n\approx-0.71$ and this algorithm terminates in a time of $O(n)$. In
contrast, spectral relaxation in which one binarizes the lowest energy
eigenvector of $J_{ij}$ only achieves an energy of $\langle
E\rangle/n\approx-0.637$ while being computationally more costly [18].
Finally, quantum methods have been proposed to solve the SK problem [19, 20].
In particular, rigorous analysis of the Quantum Approximate Optimization
Algorithm (QAOA) in Ref. [20] shows that at $p=11$ (with $p$ being the number
of layers in the quantum circuit) one can outperform the spectral relaxation
result. Although the limit of $p\rightarrow\infty$ is not established, the
results in [20] indicate convergence to the Parisi value like $1/p^{a}$ where
$a=O(1)$. However, the cost of evaluating the circuit parameters is of
$O(16^{p})$ such that one incurs an exponential cost in $1/\epsilon$ in order
to get $\epsilon$ close to the ground state. While this cost might be brought
down, it would currently be prohibitively expensive to reach the LLG results
as $p$ would have to be about $p\sim 40$. Interestingly, the LLG dynamics is
simply the semi-classical approximation of the real quantum dynamics of these
spins, suggesting that we still have much to learn about these quantum
algorithms.

As in the case of quantum annealing, our method is important because it offers
the possibility of being realized experimentally in the lab by means of arrays
of magnetic tunnel junction devices, whose physical dynamics is governed
precisely by the LLG equations. These devices operate at nanosecond timescales
and consume little power. Interesting future research directions would be to
understand the statistical limits of the robustness of our method across
different choices of the system’s parameters, as well as the inclusion of
3-spin interactions to solve 3-SAT problems, which are of central importance
in the vast context that stretches from artificial intelligence, to
cryptography, and operations research.

_Acknowledgements._ We thank Jonathan Z. Sun and Daniel L. Stein for helpful
discussions. This work was supported in part through the NYU IT High
Performance Computing resources, services, and staff expertise. This research
was supported by the Office of Naval Research (ONR) under Award No.
N00014-23-1-2771. The Flatiron Institute is a division of the Simons
Foundation. We acknowledge support from Air Force Office of Scientific
Research(AFOSR): Grant FA9550-21-1-0236.

## References

  * Mezard and Montanari [2009] M. Mezard and A. Montanari,  _Information, Physics, and Computation_ (Oxford University Press, Inc., USA, 2009).
  * Krzakala  _et al._ [2015] F. Krzakala, F. Ricci-Tersenghi, L. Zdeborova, R. Zecchina, E. W. Tramel, and L. F. Cugliandolo,  _Statistical Physics, Optimization, Inference, and Message-Passing Algorithms: Lecture Notes of the Les Houches School of Physics: Special Issue, October 2013_ (Oxford University Press, 2015).
  * Moore and Mertens [2011] C. Moore and S. Mertens,  _The Nature of Computation_ (Oxford University Press, 2011).
  * Fu and Anderson [1986] Y. Fu and P. W. Anderson, Application of statistical mechanics to np-complete problems in combinatorial optimisation, Journal of Physics A: Mathematical and General 19, 1605 (1986).
  * Kirkpatrick  _et al._ [1983] S. Kirkpatrick, C. D. Gelatt, and M. P. Vecchi, Optimization by simulated annealing, Science 220, 671 (1983).
  * Mohseni  _et al._ [2022] N. Mohseni, P. L. McMahon, and T. Byrnes, Ising machines as hardware solvers of combinatorial optimization problems, Nature Rev. Phys. 4, 363 (2022).
  * Sherrington and Kirkpatrick [1975] D. Sherrington and S. Kirkpatrick, Solvable model of a spin-glass, Phys. Rev. Lett. 35, 1792 (1975).
  * Arora  _et al._ [2005] S. Arora, E. Berger, E. Hazan, G. Kindler, and M. Safra, On non-approximability for quadratic programs, in  _Proceedings - 46th Annual IEEE Symposium on Foundations of Computer Science, FOCS 2005_, Proceedings - Annual IEEE Symposium on Foundations of Computer Science, FOCS (2005) pp. 206–215, 46th Annual IEEE Symposium on Foundations of Computer Science, FOCS 2005 ; Conference date: 23-10-2005 Through 25-10-2005.
  * Parisi [1980] G. Parisi, The order parameter for spin glasses: a function on the interval 0-1, Journal of Physics A: Mathematical and General 13, 1101 (1980).
  * Glauber [1963] R. J. Glauber, Time‐Dependent Statistics of the Ising Model, Journal of Mathematical Physics 4, 294 (1963).
  * Liftshitz and Pitaevskii [1980] E. M. Liftshitz and L. P. Pitaevskii, Statistical physics part 2, (1980).
  * Brown [1963] W. F. Brown, Thermal Fluctuations of a Single-Domain Particle, Physical Review 130, 1677 (1963).
  * García-Palacios and Lázaro [1998] J. L. García-Palacios and F. J. Lázaro, Langevin-dynamics study of the dynamical properties of small magnetic particles, Phys. Rev. B 58, 14937 (1998).
  * Ament  _et al._ [2016] S. Ament, N. Rangarajan, A. Parthasarathy, and S. Rakheja, Solving the stochastic Landau-Lifshitz-Gilbert-Slonczewski equation for monodomain nanomagnets : A survey and analysis of numerical techniques, arXiv: Computational Engineering, Finance, and Science (2016).
  * Kent and Worledge [2015] A. D. Kent and D. C. Worledge, A new spin on magnetic memories, Nature Nanotechnology 10, 187 (2015).
  * Parisi  _et al._ [1993] G. Parisi, F. Ritort, and F. Slanina, Several results on the finite-size corrections in the sherrington-kirkpatrick spin-glass model, Journal of Physics A: Mathematical and General 26, 3775 (1993).
  * Montanari [2019] A. Montanari, Optimization of the sherrington-kirkpatrick hamiltonian, in  _2019 IEEE 60th Annual Symposium on Foundations of Computer Science (FOCS)_ (2019) pp. 1417–1433.
  * Aizenman  _et al._ [1987] M. Aizenman, J. L. Lebowitz, and D. Ruelle, Some rigorous results on the sherrington-kirkpatrick spin glass model, Communications in Mathematical Physics 112, 3 (1987).
  * Harrigan  _et al._ [2021] M. P. Harrigan, K. J. Sung, M. Neeley, K. J. Satzinger, F. Arute, K. Arya, J. Atalaya, J. C. Bardin, R. Barends, S. Boixo, M. Broughton, B. B. Buckley, D. A. Buell, B. Burkett, N. Bushnell, Y. Chen, Z. Chen, B. Chiaro, R. Collins, W. Courtney, S. Demura, A. Dunsworth, D. Eppens, A. Fowler, B. Foxen, C. Gidney, M. Giustina, R. Graff, S. Habegger, A. Ho, S. Hong, T. Huang, L. B. Ioffe, S. V. Isakov, E. Jeffrey, Z. Jiang, C. Jones, D. Kafri, K. Kechedzhi, J. Kelly, S. Kim, P. V. Klimov, A. N. Korotkov, F. Kostritsa, D. Landhuis, P. Laptev, M. Lindmark, M. Leib, O. Martin, J. M. Martinis, J. R. McClean, M. McEwen, A. Megrant, X. Mi, M. Mohseni, W. Mruczkiewicz, J. Mutus, O. Naaman, C. Neill, F. Neukart, M. Y. Niu, T. E. O’Brien, B. O’Gorman, E. Ostby, A. Petukhov, H. Putterman, C. Quintana, P. Roushan, N. C. Rubin, D. Sank, A. Skolik, V. Smelyanskiy, D. Strain, M. Streif, M. Szalay, A. Vainsencher, T. White, Z. J. Yao, P. Yeh, A. Zalcman, L. Zhou, H. Neven, D. Bacon, E. Lucero, E. Farhi, and R. Babbush, Quantum approximate optimization of non-planar graph problems on a planar superconducting processor, Nature Physics 17, 332 (2021).
  * Farhi  _et al._ [2022] E. Farhi, J. Goldstone, S. Gutmann, and L. Zhou, The Quantum Approximate Optimization Algorithm and the Sherrington-Kirkpatrick Model at Infinite Size, Quantum 6, 759 (2022).