# Electron states, bound to a texture in a Néel antiferromagnet

N. Davier<EMAIL_ADDRESS>Laboratoire de Physique Théorique,
Université de Toulouse, CNRS, UPS, France  R. Ramazashvili<EMAIL_ADDRESS>tlse.fr Laboratoire de Physique Théorique, Université de Toulouse, CNRS, UPS,
France

(August 28, 2024)

###### Abstract

We study electron states, bound to topological textures such as skyrmions and
domain walls in a Néel antiferromagnet. In certain limits, we find the
dependence of bound states on the geometry of the texture, and estimate the
bound-state contribution to its energy. This contribution proves significant
compared with the purely magnetic energy, and thus substantially influences
the shape of the texture. It also considerably shifts the transition line
between the modulated and the uniform phase in favor of the latter.

##  I Introduction

Topological textures such as skyrmions and domain walls appear in various
contexts from particle physics to cold atoms and, as in the present work, to
solid-state magnetism Sky_Gobel_2021 . They also hold promise for spintronics
Jiang2017 ; Sky_Back_2020 ; Sky_Gobel_2021  as near-future information
carriers in racetrack memory devices Sky_Gobel_2021 ; Racetrack_Parkin2008 ;
Racetrack_Parkin2015 ; Fert2013 .

Both ferro- and antiferromagnets have been envisaged as materials for such
applications. However, in ferromagnets, Skyrmion Hall Effect (SHE)
Sky_Gobel_2021  makes a skyrmion deviate from driving current, likely leading
to pinning at the racetrack edges and hence information loss.
Antiferromagnetic skyrmions are a promising alternative in this regard
Sky_Gobel_2021 ; Wang ; AF_Sk_Bessarab_2019 ; AF_Sk_Zhang_2016 ;
AF_Sk_Gobel_2017 ; AF_sk_Legrand2020 ; AF_Sk_Barker_2016 . Firstly, zero net
magnetization improves robustness: being composed of two intertwined
ferromagnetic skyrmions with mutually reversed spins, antiferromagnetic
skyrmion experiences no SHE, allowing rectilinear displacement. Also,
antiferromagnetic skyrmions are expected to develop higher velocities than
their ferromagnetic counterparts. Note, however, that skyrmions in an
antiferromagnet are yet to be experimentally observed Sky_Gobel_2021 ;
AF_Sk_Bessarab_2019 ; AF_Sk_Zhang_2016 ; AF_Sk_Gobel_2017 ; AF_sk_Legrand2020
; AF_Sk_Barker_2016 ; Bogdanov_2002 .

Using antiferromagnetic skyrmions in future devices requires understanding
their interplay with other subsystems of a solid, most notably with band
electrons. Some steps in this direction have already been made: R. Cheng and
Q. Niu derived quasiclassical equations of motion for the electron coordinate,
momentum and spin near a texture in a Néel antiferromagnet Cheng_2012 . Using
a different approach, we have recently derived an effective-mass electron
Hamiltonian in the presence of a texture Davier_2023 . The latter work also
showed that, for certain locations of the electron band extrema, a texture
produces a peculiar (and anomalously strong) spin-orbit coupling. For a
skyrmion texture, this coupling generates skyrmion-electron bound states. In
the presence of dopant carriers, such bound states turn the skyrmion into a
charged particle. Reference Davier_2023  studied a tractable example of the
Belavin-Polyakov (BP) skyrmion Belavin_Polyakov_1975 ; Rajaraman , and found
the bound-state energy dependence on the skyrmion radius $R$ in the regimes
where the bound states are shallow.

Here, we extend the results of Ref. Davier_2023  in a number of ways: (i) For
the BP skyrmion, we gain quantitative insight into the range where the
skyrmion-electron bound states become profound, that is of energies comparable
with the gap $\Delta$ in the electron spectrum of the Néel state. (ii) We also
explore another profile, that becomes relevant in the presence of single-ion
anisotropy and Dzyaloshinskii-Moriya coupling: that of a domain-wall (DW)
skyrmion. Contrary to the BP skyrmion, its DW counterpart is defined by two
lengthscales rather than one: the radius $R$ and the domain wall width $w$. We
study both shallow and deep bound states of an electron and a DW skyrmion.
(iii) Locally, a DW skyrmion boundary is a domain wall. Thus, to gain insight,
we also study band electrons in the presence of a straight domain wall, and
show that it produces electron states that are bound in the transverse
direction, and itinerant along the wall. This calculation also allows us to
understand the genesis of electron bound states for a DW skyrmion. (iv) Last
but not least, at half-filling we address the influence of the skyrmion-
electron bound states on the energy balance of a domain wall.

The key results of our work are as follows: The skyrmion-electron bound states
and their energies are sensitive to the skyrmion ‘shape’. Moreover, the bound-
state contribution to the skyrmion energy proves dominant compared with its
purely ‘magnetic’ counterpart – and thus plays a key role in selecting an
‘optimal’ skyrmion profile. In particular, the BP profile with a sufficiently
small skyrmion radius induces no bound states, and hence no additional energy
cost, which favors smaller BP skyrmions. For DW skyrmions, the bound-state
contribution proves dominant, as well. We show how it enhances the stability
range of the uniform state, and shifts the transition to the modulated phase
Bogdanov_1989  to substantially higher values of the Dzyaloshinskii-Moriya
coupling.

The article has the following structure. Section II is devoted to
generalities: derivation of an electron Hamiltonian in the presence of a
texture, and of its effective-mass (low-energy) limit, with the focus on a
specific location of the band extrema. This section expands on the results of
Ref. Davier_2023  – both to make the presentation self-contained and to
provide some additional details. In Sec. III we use the low-energy Hamiltonian
to treat shallow bound states for the Belavin-Polyakov skyrmion. The deep
bound-state limit is next studied in Sec. IV, allowing us to complete the
picture of the bound state evolution as a function of the BP skyrmion radius.
In Sec. V, we address the bound-state problem for a DW skyrmion, exploring
both shallow and deep states. In Sec. VI, we estimate the contribution of
bound states to the skyrmion energy as a function of the skyrmion shape. We
also study electron states, bound to a straight domain wall, and evaluate
their contribution to the domain-wall energy. Which allows us to demonstrate
how texture-bound electron states substantially shift the transition line
between the uniform and the modulated phases Bogdanov_1989 . The article
concludes with a discussion in Sec. VII, while some of the technical details
are given in the Appendices.

##  II The Hamiltonian

###  II.1 General arguments and the uniform state

We start by writing down the effective Hamiltonian of a band electron in a
collinear Néel antiferromagnet in the second-quantization formalism:

$\mathcal{H}=\int
d^{2}r\;\Psi^{\dagger}(\mathbf{r})\left[\varepsilon\left(\hat{\mathbf{p}}\right)+J_{e}\mathbf{M}(\mathbf{r})\cdot{\bm{\sigma}}\right]\Psi(\mathbf{r})$
(1)

with $\Psi^{\dagger}(\mathbf{r})$ and $\Psi(\mathbf{r})$ the electron
creation/annihilation operators at point $\mathbf{r}$, and
$\varepsilon\left(\hat{\mathbf{p}}\right)$ the electron dispersion in the
absence of Néel order. The second term
$J_{e}\mathbf{M}(\mathbf{r})\cdot{\bm{\sigma}}$ accounts for the exchange
interaction between the spontaneous local magnetization
$\mathbf{M}(\mathbf{r})$ and the electron spin, represented by the triad of
Pauli matrices ${\bm{\sigma}}$. Without loss of generality, hereafter we
consider a square lattice with spacing $a_{0}$. Néel order can be defined via
nearly opposite sublattice magnetizations $\mathbf{M}_{1}$ and
$\mathbf{M}_{2}$, from which one can build the local Néel order parameter
$\mathbf{n}(\mathbf{r})=\frac{\mathbf{M}_{1}-\mathbf{M}_{2}}{2M_{s}}$ and the
effective local magnetization
$\mathbf{m}(\mathbf{r})=\frac{\mathbf{M}_{1}+\mathbf{M}_{2}}{2M_{s}}$, where
$M_{s}$ is the saturation magnetization. In the Néel state, the order
parameter $\mathbf{n}(\mathbf{r})$ is large compared with the effective local
magnetization: $\lVert\mathbf{n}\lVert\simeq 1\gg\lVert\mathbf{m}\lVert$;
moreover, $\mathbf{n}(\mathbf{r})$ varies smoothly in space and slowly in
time, while $\mathbf{m}(\mathbf{r})$ fluctuates wildly. Therefore, hereafter
we ignore the effective local magnetization $\mathbf{m}(\mathbf{r})$
altogether and approximate the magnetization by its dominant Fourier harmonic
as per $\mathbf{M}(\mathbf{r})\simeq
M_{s}\mathbf{n}(\mathbf{r})e^{i\mathbf{Q}\cdot\mathbf{r}}$ with
$\mathbf{Q}=\left(\pm\frac{\pi}{a_{0}},\pm\frac{\pi}{a_{0}}\right)$.

It is convenient to begin with the uniform Néel state, and choose the
magnetization to point along the $z$ axis in spin space:
$\mathbf{M}(\mathbf{r})\simeq
M_{s}\mathbf{e}_{z}e^{i\mathbf{Q}\cdot\mathbf{r}}$. Upon introducing the
notation $\Delta=J_{e}M_{s}$, the Hamiltonian reads

$\mathcal{H}=\int
d^{2}r\;\Psi^{\dagger}(\mathbf{r})\left[\varepsilon\left(\hat{\mathbf{p}}\right)+\Delta\sigma_{z}\,e^{i\mathbf{Q}\cdot\mathbf{r}}\right]\Psi(\mathbf{r}).$
(2)

Following transition to reciprocal space, the Hamiltonian takes the form

$\mathcal{H}=\sum_{\mathbf{p}\in\text{MBZ}}\begin{pmatrix}\Psi^{\dagger}_{\mathbf{p}}&\Psi^{\dagger}_{\mathbf{p+Q}}\end{pmatrix}\begin{pmatrix}\varepsilon_{\mathbf{p}}&\Delta\sigma_{z}\\\
\Delta\sigma_{z}&\varepsilon_{\mathbf{p+Q}}\end{pmatrix}\begin{pmatrix}\Psi_{\mathbf{p}}\\\
\Psi_{\mathbf{p+Q}}\end{pmatrix}.$ (3)

That is, Néel order couples electron states at any two momenta $\mathbf{p}$
and $\mathbf{p+Q}$ Kulikov1984REVIEWSOT , reducing the Brillouin zone in the
paramagnetic state to its half, as shown in Fig. 1.

Figure 1: Magnetic Brillouin zone of the Néel state (cyan square), and the
paramagnetic Brillouin zone (white square). Points $\Sigma_{i}$ with small
valleys, described by the expansion (LABEL:Developpement_Ep), are shown
symbolically in light blue. The bound sates at these valleys are spin-
polarized, with polarization depending on the skyrmion profile. The
polarization above corresponds to a Néel BP skyrmion, see main text for
details.

The Hamiltonian is block-diagonal, and can be easily diagonalized to produce
two doubly-degenerate electron bands

$E_{\pm}=\varepsilon_{+}\pm\sqrt{\Delta^{2}+\varepsilon_{-}^{2}},\hskip
19.91684pt\varepsilon_{\pm}=\frac{\varepsilon_{\mathbf{p}}\pm\varepsilon_{\mathbf{p+Q}}}{2}.$
(4)

These bands are separated by a gap of width $2\Delta$, opening at crossing
points $\varepsilon_{\mathbf{p}}=\varepsilon_{\mathbf{p+Q}}$. The Néel order
thus turns a half-filled metal into an insulator. This elementary treatment of
the uniform Néel state allows us to understand the overall electron band
structure, and puts us in a position to study a non-uniform texture
$\mathbf{n}(\mathbf{r})$.

###  II.2 A non-uniform texture

Consider the Néel order parameter $\mathbf{n}(\mathbf{r})$, varying smoothly
in space – that is, on a scale large compared with $a_{0}$. Because of this
nonuniformity, Hamiltonian (2) can no longer be easily diagonalized. However,
we can gain insight by rendering $\mathbf{n}(\mathbf{r})$ uniform via a spin
rotation $U(\mathbf{r})$ acting as per

$U^{\dagger}(\mathbf{r})\Big{[}\mathbf{n}(\mathbf{r})\cdot\bm{\sigma}\Big{]}U(\mathbf{r})=\sigma^{z}.$
(5)

Upon transformation $\Psi\rightarrow U^{\dagger}\Psi$, the Hamiltonian (2)
reads

$\mathcal{H}=\int
d^{2}r\;\Psi^{\dagger}(\mathbf{r})\left[U^{\dagger}(\mathbf{r})\varepsilon\left(\hat{\mathbf{p}}\right)U(\mathbf{r})+\Delta\sigma^{z}e^{i\mathbf{Q}\cdot\mathbf{r}}\right]\Psi(\mathbf{r}).$

The price for reducing the second term in the square brackets to that of Eq.
(2) amounts to the appearance of a Peierls substitution Volovik1987  in the
kinetic term, since

$U^{\dagger}\hat{p}_{i}U=\hat{p}_{i}+\mathbf{A}_{i}\cdot\bm{\sigma}$ (6)

where $\mathbf{A}_{i}\cdot\bm{\sigma}=A_{i}^{\alpha}\sigma^{\alpha}=-i\hbar
U^{\dagger}(\partial_{i}U)$, with $i$ the spatial indices and $\alpha$ the
spin indices (summation over repeated indices is implied). Various properties
of this gauge field will be discussed later. Switching back to real space, the
single-particle Hamiltonian, acting in the basis
$\left(\Psi_{\mathbf{p}},\Psi_{\mathbf{p+Q}}\right)$, reads

$H=\begin{pmatrix}\varepsilon_{\mathbf{p}}\left(\hat{\mathbf{p}}+\mathbf{A}\cdot\bm{\sigma}\right)&\Delta\sigma^{z}\\\
\Delta\sigma^{z}&\varepsilon_{\mathbf{p+Q}}\left(\hat{\mathbf{p}}+\mathbf{A}\cdot\bm{\sigma}\right)\end{pmatrix}$
(7)

with
$\varepsilon_{\mathbf{p}}\left(\hat{\mathbf{p}}+\mathbf{A}\cdot\bm{\sigma}\right)=U^{\dagger}\varepsilon\left(\hat{\mathbf{p}}\right)U$.
Because of the Peierls substitution, the terms on the diagonal in Eq. (7) are
no longer proportional to the unit matrix in spin space.

With the band structure of the previous section in mind, and with Hamiltonian
(7) at hand, we can now focus on the electron states near the band extrema.

###  II.3 Low-energy effective Hamiltonian

We wish to derive an effective-mass Hamiltonian Kittel1963  describing an
electron near the conduction band minimum at point $\mathbf{p}_{0}$. To this
end, we replace the electron momentum $\mathbf{p}$ by
$\hat{\mathbf{p}}=-i\hbar\bm{\nabla}$, and expand the dispersion relations
$\varepsilon(\mathbf{p})$ and $\varepsilon(\mathbf{p+Q})$, denoted as
$\varepsilon_{\mathbf{p}_{0}}(\mathbf{p})$ and
$\varepsilon_{\mathbf{p}_{0}+\mathbf{Q}}(\mathbf{p})$, around the given
conduction band minimum $\mathbf{p}_{0}$. Note that $\mathbf{p}_{0}$ must be a
crossing point
$\varepsilon(\mathbf{p}_{0})=\varepsilon(\mathbf{p}_{0}+\mathbf{Q})$ of the
uniform state. With the help of a unitary operator 111 Note that this operator
is not unique since, at point $\mathbf{p}_{0}$, Hamiltonian (3) commutes with
$\sigma^{z}$. This implies that $P$ is defined up to any two independent spin
rotations around $\mathbf{e}_{z}$, for the conduction and the valence band,
respectively. Eq. (8) corresponds to the simplest choice, with no spin
rotation around $\mathbf{e}_{z}$ in either of the two bands.

$P=\frac{1}{\sqrt{2}}\begin{pmatrix}1\\!\\!1&1\\!\\!1\\\
\sigma^{z}&-\sigma^{z}\end{pmatrix}$ (8)

we now switch to the basis where Hamiltonian (7) becomes diagonal at
$\mathbf{p}_{0}$:

$H\to P^{\dagger}HP=\begin{pmatrix}\hat{\bm{\gamma}}_{+}+\Delta
1\\!\\!1&\hat{\bm{\gamma}}_{-}\\\
\hat{\bm{\gamma}}_{-}&\hat{\bm{\gamma}}_{+}-\Delta 1\\!\\!1\end{pmatrix},$ (9)

with

$\hat{\bm{\gamma}}_{\pm}(\hat{\mathbf{p}})=\frac{1}{2}\left[\varepsilon_{\mathbf{p}_{0}}\left(\hat{\mathbf{p}}+\mathbf{A}\cdot\bm{\sigma}\right)\pm\sigma^{z}\varepsilon_{\mathbf{p_{0}+Q}}\left(\hat{\mathbf{p}}+\mathbf{A}\cdot\bm{\sigma}\right)\sigma^{z}\right].$
(10)

Figure 2: Deriving Hamiltonian $H_{eff}$ in Eq. (13) from the full $4\times 4$
Hamiltonian of Eq. (9). Hamiltonian (13) provides an effective description
near the band extremum.

From Hamiltonian (9) acting on a bispinor, we now wish to construct an
effective Hamiltonian, that would describe only states near the bottom of
conduction band (see Fig. 2), and thus would act only on a single spinor. This
is akin to taking the Dirac Hamiltonian to its non-relativistic Pauli-
Schrödinger limit Landau_vol4 ; Ryder1996Quantum . To this end, we apply
Hamiltonian (9) to the bispinor $\left(\Phi,\chi\right)$. We then consider an
electron at energy $E=\Delta+\varepsilon$ near the conduction band minimum
$+\Delta$, with $|\varepsilon|\ll\Delta$. This yields two equations

$\begin{split}\hat{\bm{\gamma}}_{-}\chi&=\left(\varepsilon-\hat{\bm{\gamma}}_{+}\right)\Phi,\\\
\hat{\bm{\gamma}}_{-}\Phi&=\left(2\Delta+\varepsilon-\hat{\bm{\gamma}}_{+}\right)\chi.\end{split}$
(11)

Treating $\hat{\bm{\gamma}}_{\pm}(\hat{\mathbf{p}})$ as small against
$\Delta$, we limit ourselves to leading terms of their momentum expansion near
$\mathbf{p}_{0}$. To first order in $\frac{\hat{\bm{\gamma}}_{-}}{\Delta}\ll
1$, one finds $\chi\simeq\frac{\hat{\bm{\gamma}}_{-}}{2\Delta}\Phi$, which
means $\Phi\gg\chi$. This allows us to write down the effective Schrödinger
equation acting only on the spinor $\Phi$:

$\left(\hat{\bm{\gamma}}_{+}+\frac{\hat{\bm{\gamma}}_{-}^{2}}{2\Delta}\right)\Phi=\varepsilon\Phi,$
(12)

which yields the sought effective Hamiltonian to first order in
$\frac{\varepsilon}{\Delta}\ll 1$:

$H_{eff}=\hat{\bm{\gamma}}_{+}+\frac{\hat{\bm{\gamma}}_{-}^{2}}{2\Delta}.$
(13)

The leading terms in the momentum expansion of $\hat{\bm{\gamma}}_{\pm}$
depend on the location of the minimum $\mathbf{p}_{0}$ in the Brillouin zone.
Here we focus on points $\Sigma_{i}$: the face centers of the magnetic
Brillouin zone in Fig. 1. The peculiarity of these points is the presence of a
linear term in the momentum expansion of
$\varepsilon_{\mathbf{p_{0}}}\left(\mathbf{p}\right)$ and
$\varepsilon_{\mathbf{p_{0}+Q}}\left(\mathbf{p}\right)$:

$\begin{split}&\varepsilon_{\mathbf{p_{0}}}\left(\mathbf{p}\right)=vp_{y}+\frac{p_{i}^{2}}{2m_{i}}+\mathcal{O}\left(\left\|\mathbf{p}^{3}\right\|\right),\\\
&\varepsilon_{\mathbf{p_{0}+Q}}\left(\mathbf{p}\right)=-vp_{y}+\frac{p_{i}^{2}}{2m_{i}}+\mathcal{O}\left(\left\|\mathbf{p}^{3}\right\|\right).\end{split}$
(14)

Here $v$ is the Fermi velocity of the parent paramagnetic state along the
$\hat{y}$ direction, and $m_{i}$ are the effective masses along $\hat{x}$ and
$\hat{y}$, as shown in Fig. 1.

Equations (LABEL:Developpement_Ep) yield the expansion of
$\hat{\bm{\gamma}}_{\pm}(\hat{\mathbf{p}})$ in Eq. (10) to second order in
momentum:

$\begin{split}&\hat{\bm{\gamma}}_{+}(\hat{\mathbf{p}})=v\mathbf{A}_{y}^{\parallel}\cdot\bm{\sigma}+\frac{1}{2m_{i}}\left[\left(\hat{p}_{i}+A_{i}^{z}\sigma^{z}\right)^{2}+\left(A_{i}^{\parallel}\right)^{2}\right],\\\
&\hat{\bm{\gamma}}_{-}(\hat{\mathbf{p}})=v(\hat{p}_{y}+A_{y}^{z}\sigma^{z})+\frac{1}{2m_{i}}\\{\hat{p}_{i},\mathbf{A}_{i}^{\parallel}\cdot\bm{\sigma}\\},\end{split}$
(15)

where $\mathbf{A}_{i}^{\parallel}=(A_{i}^{x},A_{i}^{y},0)^{t}$ are the gauge-
field components in the $xy$-plane in spin space. Which, in turn, allows us to
rewrite Hamiltonian (13) to second order in momentum:

$H=\frac{1}{2m_{i}^{*}}\left(\hat{p}_{i}+A_{i}^{z}\sigma^{z}\right)^{2}+v\mathbf{A}_{y}^{\parallel}\cdot\bm{\sigma}+\frac{\left(A_{i}^{\parallel}\right)^{2}}{2m_{i}},$
(16)

where the second term in Eq. (13) renormalizes the effective mass $m_{y}^{*}$
relative to its band value in the expansion (LABEL:Developpement_Ep) as per

$m_{y}^{-1}\;\longrightarrow\;(m_{y}^{*})^{-1}=m_{y}^{-1}+\frac{v^{2}}{\Delta},$
(17)

with $\frac{m_{y}^{*}}{m_{y}}\sim\frac{\Delta}{\varepsilon_{F}}\ll 1$. As we
will see below, it is thanks to this renormalization that the kinetic term
$\frac{1}{2m_{y}^{*}}\left(\hat{p}_{y}+A_{y}^{z}\sigma^{z}\right)^{2}$ becomes
comparable with the dominant potential term
$v\mathbf{A}_{y}^{\parallel}\cdot\bm{\sigma}$ at smaller texture length scales
$a\ll L\ll\xi$.

Hamiltonian (16) admits further simplification: within continuum description,
the Néel order parameter $\mathbf{n}$ of an antiferromagnetic texture varies
on a length scale $L$, that is large compared with the lattice spacing
$a_{0}$. The latter is of the order of the ‘effective lattice spacing’
$a\equiv\frac{\hbar}{m_{y}v}$. As the reader will see below, the two key
length scales in the problem are $a$ and the ‘coherence length’
$\xi\equiv\frac{\hbar v}{\Delta}$. Their ratio is
$\frac{a}{\xi}=\frac{\Delta}{\varepsilon_{F}}\ll 1$, with
$\varepsilon_{F}\equiv m_{y}v^{2}$ the characteristic electron bandwidth. For
an electron wave function varying on the characteristic length scale
$L_{\psi}\gg a$, we find

$vA_{i}^{\alpha}\sim\Delta\frac{\xi}{L}\gg\frac{\left(A_{i}^{\alpha}\right)^{2}}{2m_{i}},\frac{\\{\hat{p}_{i},A_{i}^{\alpha}\\}}{2m_{i}}\sim\Delta\left\\{\frac{\xi
a}{L^{2}},\frac{\xi a}{LL_{\psi}}\right\\}.$ (18)

As a result, the last term in Eq. (16) can be safely omitted at all relevant
texture length scales below, and the leading terms of Hamiltonian (13) to
second order in momentum thus read

$H=\frac{1}{2m_{i}^{*}}\left(\hat{p}_{i}+A_{i}^{z}\sigma^{z}\right)^{2}+v\mathbf{A}_{y}^{\parallel}\cdot\bm{\sigma}.$
(19)

Notice that the effective mass $m_{i}^{*}$ above is substantially anisotropic:
$m_{x}^{*}=m_{x}$ is of the order of the electron band mass or greater (on a
square lattice with only the nearest-neighbour hopping, $m_{x}$ is infinite),
and thus $m_{x}^{*}\gg m_{y}^{*}$. While mass anisotropy at points $\Sigma$ is
generic as dictated by their position at the magnetic Brillouin zone boundary,
strong mass anisotropy in Hamiltonian (19) can be traced back to separation of
scales: antiferromagnetic gap $\Delta$ in the electron spectrum being small
relative to the band width $\varepsilon_{F}$.

As we show below, the second term
$v\mathbf{A}_{y}^{\parallel}\cdot\bm{\sigma}$ of Hamiltonian (19) tends to
produce electron states, bound to topological textures such as skyrmions and
domain walls. On the one hand, in a Néel antiferromagnet, such a term owes its
existence to lower symmetry of band extrema $\Sigma_{i}$ in the Brillouin
zone: at higher-symmetry extrema (e.g. points $\Gamma$ and $X$ in Fig. 1),
such a term is not allowed. On the other hand, in a ferromagnet, such a term
is not allowed at all, which underlines a key difference between ferro- and
antiferromagnets. We briefly compare the two problems in Appendix A.

Hamiltonian (19) may appear odd to a reader used to working with
electromagnetic field. Thus we now take a closer look at transformation
properties of the gauge field $\mathbf{A}$.

###  II.4 Gauge transformations

Unitary operator $U(\mathbf{r})$ of Eq. (5) is defined up to a spin rotation
$V_{z}(\mathbf{r})=e^{i\chi(\mathbf{r})\sigma^{z}}$ around the $z$ axis, with
a single-valued field $\chi(\mathbf{r})$. Indeed, such a rotation leaves Eq.
(5) intact:

$V_{z}^{\dagger}U^{\dagger}\Big{[}\mathbf{n}\cdot\bm{\sigma}\Big{]}UV_{z}=V_{z}^{\dagger}\sigma^{z}V_{z}=\sigma^{z},$
(20)

which means that $V_{z}(\mathbf{r})$ amounts to a gauge transformation. Under
such a transformation, $\mathbf{A}_{i}\cdot\bm{\sigma}$ of Eq. (6) varies as
per

$\begin{split}\tilde{\mathbf{A}_{i}}\cdot\bm{\sigma}&=-i\hbar
V_{z}^{\dagger}U^{\dagger}\partial_{i}\left(UV_{z}\right)\\\
&=V_{z}^{\dagger}\mathbf{A}_{i}\cdot\bm{\sigma}V_{z}-i\hbar
V_{z}^{\dagger}\partial_{i}V_{z}\\\
&=V_{z}^{\dagger}\mathbf{A}_{i}^{\parallel}\cdot\bm{\sigma}V_{z}+A_{i}^{z}\sigma^{z}+\hbar\partial_{i}\chi\sigma^{z}.\end{split}$
(21)

That is, the $A_{i}^{z}$ component transforms as electromagnetic vector
potential, while the transverse component
$\mathbf{A}_{i}^{\parallel}=(A_{i}^{x},A_{i}^{y},0)^{t}$ simply rotates by
$2\chi(\mathbf{r})$ around the $z$ axis in spin space. While the different
spin components $\mathbf{A}_{i}^{\alpha}\sigma^{\alpha}$ do not commute, they
split into $A_{i}^{z}\sigma^{z}$ and
$\mathbf{A}_{i}^{\parallel}\cdot\bm{\sigma}$, each obeying its distinct
abelian transformation law.

Another transformation that leaves Eq. (5) invariant is a spin rotation
$W(\mathbf{r})=e^{i\phi(\mathbf{r})\mathbf{n}\cdot\bm{\sigma}}$ around the
local $\mathbf{n}(\mathbf{r})$:

$U^{\dagger}W^{\dagger}\Big{[}\mathbf{n}\cdot\bm{\sigma}\Big{]}WU=U^{\dagger}\mathbf{n}\cdot\bm{\sigma}U=\sigma^{z},$
(22)

which implies gauge transformation

$\tilde{\mathbf{A}_{i}}\cdot\bm{\sigma}=\mathbf{A}_{i}\cdot\bm{\sigma}-i\hbar
U^{\dagger}\left[W^{\dagger}\partial_{i}W\right]U.$ (23)

Just as for transformation $V_{z}$ above, a simple calculations shows that,
under $W$, $\mathbf{A}_{i}^{\|}$ undergoes an in-plane rotation by $2\phi$,
while $A_{i}^{z}$ transforms as electromagnetic vector potential (see Appendix
B).

Note that, for infinitesimal $\chi(\mathbf{r})$ and $\phi(\mathbf{r})$, the
relations (LABEL:transformation_de_jauge) and (23) read, respectively

$\delta\mathbf{A}_{i}\cdot\bm{\sigma}=i\hbar\chi[\mathbf{A}_{i}\cdot\bm{\sigma},\sigma^{z}]+\hbar\partial_{i}\chi\sigma^{z}$
(24)

and 222Here we use the identity
$-i\hbar\partial_{i}\left(U^{\dagger}\mathbf{n}\cdot\bm{\sigma}U\right)=[\sigma^{z},\mathbf{A}_{i}\cdot\bm{\sigma}]-i\hbar
U^{\dagger}\partial_{i}\mathbf{n}\cdot\bm{\sigma}U=0.$

$\begin{split}\delta\mathbf{A}_{i}\cdot\bm{\sigma}&=\hbar\phi\,U^{\dagger}[\partial_{i}\mathbf{n}\cdot\bm{\sigma}]U+\hbar\partial_{i}\phi\sigma^{z}\\\
&=i\hbar\phi\,[\mathbf{A}_{i}\cdot\bm{\sigma},\sigma^{z}]+\hbar\partial_{i}\phi\sigma^{z}\end{split}$
(25)

and thus become identical. Invariant under these two transformations,
Hamiltonians (16) and (19) enjoy the peculiar gauge symmetry above.

To perform calculations for a given texture, we will have to choose a concrete
$U(\mathbf{r})$ – that is, to fix a gauge. We will define $\mathbf{n}$ via its
polar angle $\theta$ and azimuthal angle $\phi$ as per

$\mathbf{n}=\left(\sin\theta\cos\phi,\sin\theta\sin\phi,\cos\theta\right)^{t},$
(26)

and choose Tatara_2008

$U(\mathbf{r})=\mathbf{m}\cdot\bm{\sigma}$ (27)

with $\mathbf{m}$ being the bisector between $\mathbf{n}$ and
$\mathbf{e}_{z}$. That is, $\mathbf{m}\cdot\bm{\sigma}$ amounts to a
$\pi$-rotation around $\mathbf{m}$. Such a choice makes calculations
straightforward with the help of the textbook Pauli-matrix property
$(\mathbf{n}\cdot\bm{\sigma})\,(\mathbf{m}\cdot\bm{\sigma})=\mathbf{n}\cdot\mathbf{m}+i\,(\mathbf{n}\times\mathbf{m})\cdot\bm{\sigma}$,
yielding

$\mathbf{A}_{i}\cdot\bm{\sigma}=-i\hbar\,(\mathbf{m}\cdot\bm{\sigma})\,(\partial_{i}\mathbf{m}\cdot\bm{\sigma})=\hbar\,(\mathbf{m}\times\partial_{i}\mathbf{m})\cdot\bm{\sigma}.$
(28)

A simple calculation leads to the following expression:

$\hbar^{-1}\mathbf{A}_{i}=\frac{\partial_{i}\theta}{2}\begin{pmatrix}-\sin\phi\\\
\cos\phi\\\
0\end{pmatrix}+\frac{\partial_{i}\phi}{2}\begin{pmatrix}-\sin\theta\cos\phi\\\
-\sin\theta\sin\phi\\\ 1-\cos\theta\end{pmatrix}.$ (29)

As with electromagnetic field, certain combinations of the $\mathbf{A}_{i}$
components and their gradients prove gauge-invariant. In appendix C, we show
that these can be expressed via the gradient energy density
$({\bm{\nabla}}{\bf n})^{2}$ and the the skyrmion number density $({\bf
n}\cdot\left[\partial_{x}{\bf n}\times\partial_{y}{\bf n}\right])$.

Now that we constructed an effective-mass Hamiltonian for a generic
antiferromagnetic texture, and elucidated its gauge-transformation properties,
we are in a position to study an electron in the presence of a skyrmion.

##  III Belavin-Polyakov profile, shallow bound states

###  III.1 Profile properties

Note that a skyrmion has its own localized eigenmodes that, generally, must be
considered on an equal footing with the electron degrees of freedom. There is,
however, a realistic limit, where proper frequencies of the skyrmion prove to
be small relative to those of the electron motion. Then, to first
approximation, the texture can be treated as static input to the electron
problem. Once the electron eigenstates are found, a comparison of
characteristic electron frequencies with those of the proper modes of the
skyrmion allows one to define the validity range of such a calculation, as it
was done in our prior work for a skyrmion close to the Belavin-Polyakov (BP)
profile Davier_2023 .

The BP skyrmion naturally emerges in a perfectly isotropic magnet, described
by the continuum-limit energy density
$J\left(\bm{\nabla}\cdot\mathbf{n}\right)^{2}$ with stiffness $J$. In a
topological sector labelled by the winding number $Q=\pm 1,\pm 2,\pm 3,...$,
the lowest-energy solution is the BP skyrmion of radius $R$, with
$R$-independent energy $4\pi J|Q|$ Belavin_Polyakov_1975 ; Rajaraman .

The profile of a skyrmion is specified by the dependence of polar and
azimuthal angles $\theta$ and $\phi$ in Eq. (26) on the coordinates in the
plane. For the latter, we use polar coordinates $r=\sqrt{x^{2}+y^{2}}$ and
$\alpha=\arctan\frac{y}{x}$. For high-symmetry skyrmion textures such as the
BP skyrmion, the polar angle $\theta$ depends only on the radial coordinate
$r$ in the plane, while the azimuthal angle $\phi$ depends on angle $\alpha$
in the plane as per $\phi=Q\alpha+\gamma$. The offset $\gamma$ is called
helicity, while the coefficient $Q$ is the skyrmion winding number
Sky_Gobel_2021 . A skyrmion with $\gamma=0$ is often called Néel skyrmion,
while its counterpart with $\gamma=\frac{\pi}{2}$ is referred to as Bloch
skyrmion.

Here we limit ourselves to $Q=1$, and consider the BP profile defined as per
Belavin_Polyakov_1975 ; Rajaraman

$\sin\theta=\frac{2z}{1+z^{2}},\hskip 28.45274ptz=\frac{r}{R},$ (30)

depending on a single length scale, the skyrmion radius $R$. At first sight,
the definition (30) does not imply a choice of boundary conditions for the
orbital angle $\theta(r)$. Furthermore, for an antiferromagnet, the choice of
boundary conditions should be physically irrelevant since it relies on an
arbitrary choice of definition for the Néel order parameter
$\mathbf{n}=\pm\frac{\mathbf{M}_{1}-\mathbf{M}_{2}}{2M_{s}}$. However, when
looking for example at the $z$ components of the gauge field

$A_{i}^{z}=\hbar\frac{\partial_{i}\phi}{2}\left(1-\cos\theta\right),$ (31)

we see that the choice $\theta(0)=\pi$ makes the origin singular. The only
possible choice of boundary conditions is thus $\theta(0)=0$ and
$\lim_{r\to\infty}\theta(r)=\pi$, which allows to obtain
$\cos\theta=2\frac{1-z^{2}}{1+z^{2}}$, leading to

$A^{z}_{x}=\hbar\frac{-y}{R^{2}+r^{2}}\,\,,\hskip
28.45274ptA^{z}_{y}=\hbar\frac{x}{R^{2}+r^{2}}.$ (32)

Encircling the skyrmion along a large closed contour makes the electron
accumulate a geometric flux of $2\pi$: $\sigma^{z}\oint
A^{z}_{i}dl_{i}=2\pi\hbar\sigma^{z}$, which gives rise to the skyrmion Hall
Effect TSHE_Yin_2015 . This is a consequence of a more general relation
(LABEL:invariant_de_jauge): an electron encircling a skyrmion of topological
charge $Q$ accumulates a geometric flux $2\pi Q$, independent of the skyrmion
profile.

The BP profile has a remarkable property that greatly simplifies the problem.
To illustrate it, consider the expression for the $A_{y}^{x}$, produced by a
high-symmetry $Q=1$ skyrmion with $\phi=\alpha+\gamma$:

$\begin{split}A_{y}^{x}&=-\frac{\hbar}{2}\cos\alpha\sin\alpha\sin\gamma\left(\frac{\sin\theta}{r}-\frac{d\theta}{dr}\right)\\\
&-\frac{\hbar}{2}\cos\gamma\left(\cos^{2}\alpha\frac{d\theta}{dr}+\sin^{2}\alpha\frac{\sin\theta}{r}\right).\end{split}$
(33)

According to Eq. (33), variables $r$ and $\alpha$ separate for a profile,
satisfying

$\frac{d\theta}{dr}=\frac{\sin\theta}{r}.$ (34)

Remarkably, this very relation is obeyed by the $Q=1$ BP profile John_1995 ,
implying

$\mathbf{A}_{y}^{\parallel}=-\hbar\frac{\sin\theta}{2r}\mathbf{e}_{\sigma},$
(35)

where unit vector $\mathbf{e}_{\sigma}=(\cos\gamma,\sin\gamma,0)^{t}$ defines
the orientation of $\mathbf{A}_{y}^{\parallel}$ in spin space. Note that
$\mathbf{e}_{\sigma}$ only depends on the helicity, and remains uniform in
real space. Without loss of generality, this allows us to limit ourselves to a
Néel skyrmion ($\gamma=0$), for which $\mathbf{e}_{\sigma}=\mathbf{e}_{x}$ 333
For a $Q\neq 1$ skyrmion with $\phi=Q\alpha+\gamma$, the BP profile can be
defined via  $\frac{d\theta}{dr}=Q\frac{\sin\theta}{r}.$ The dominant
potential term in Eq. 19 then reads
$v\mathbf{A}_{y}^{\|}\cdot\bm{\sigma}=\frac{\xi}{R}\frac{Q\Delta}{1+z^{2}}\,\mathbf{e}_{\sigma}^{Q}\cdot\bm{\sigma},$
where $\mathbf{e}_{\sigma}^{Q}=(\cos\phi,\sin\phi,0)^{t}$, and
$\phi\equiv\gamma+(Q-1)\alpha$. Thus, for $Q\neq 1$, orientation of
$\mathbf{A}_{y}^{\|}$ in spin space becomes non-uniform. .

We will now study shallow electron states, bound to the BP skyrmion as
described by Hamiltonian (19).

###  III.2 Shallow bound states

The first thing to note is that strong anisotropy $m_{x}^{*}\gg m_{y}^{*}$
makes $y$ a ‘fast’ variable relative to $x$, which calls for the Born-
Oppenheimer approximation Tully2001 . Recalling that
$m_{y}^{*}\sim\frac{\Delta}{v^{2}}$, Hamiltonian (19) can be recast as

$\begin{split}\frac{H}{\Delta}=\frac{\xi^{2}}{2R^{2}}\left(-\frac{\partial^{2}}{\partial\tilde{y}^{2}}+\frac{\tilde{x}^{2}}{(1+z^{2})^{2}}\right)-\frac{\xi}{R}\frac{1}{1+z^{2}}\sigma^{x},\end{split}$
(36)

with $\tilde{x}_{i}=\frac{x_{i}}{R}$. Equation above allows us to demonstrate
the appearance of skyrmion-electron bound states, polarized along
$\mathbf{e}_{x}$.

Above, anticommutator $\left\\{\hat{p}_{y},A_{y}^{z}\sigma^{z}\right\\}$ has
been omitted as proportional to $\sigma^{z}$, and thus only acting in second
order of perturbation theory, the dominant term being proportional to
$\sigma^{x}$. This can be seen directly from the physics picture in Fig. 6:
the term $\left\\{\hat{p}_{y},A_{y}^{z}\sigma^{z}\right\\}$ has matrix
elements only across the gap, hence its contribution is suppressed by large
perturbation-theory denominator $\Delta$.

Considering the $|+\rangle_{x}$ component of the wave function, the
corresponding low-lying states can be evaluated by expanding the potential
around the bottom of the well. This yields the effective harmonic potential
for the ‘slow’ coordinate $\tilde{x}=\frac{x}{R}$:

$\frac{\varepsilon^{n}_{x}}{\Delta}=\frac{\xi}{R}\left(-1+\left(n+\frac{1}{2}\right)\sqrt{\frac{2\xi}{R}}+\left[1+\frac{\xi}{2R}\right]\tilde{x}^{2}\right).$
(37)

Reintroducing the kinetic energy term $\frac{\hat{p}_{x}^{2}}{2m_{x}^{*}}$
then yields the low-lying states

$\begin{split}\frac{\varepsilon^{n}_{m}}{\Delta}&=\frac{\xi}{R}\left[-1+\left(n+\frac{1}{2}\right)\sqrt{\frac{2\xi}{R}}\right]\\\
&\hskip
28.45274pt+\left(m+\frac{1}{2}\right)\sqrt{2\epsilon}\left(\frac{\xi}{R}\right)^{3/2}\sqrt{1+\frac{\xi}{2R}}\end{split}$
(38)

with $\epsilon=\frac{m_{y}}{m_{x}}\frac{a}{\xi}\lll 1$. Parameter $\epsilon$
governs the ratio $\frac{\omega_{x}}{\omega_{y}}$ of the characteristic
frequencies of the ‘slow’ and ‘fast’ degrees of freedom. The wave functions of
these states are those of a two dimensional anisotropic harmonic oscillator
with different frequencies along directions $x$ and $y$. These wave functions
are thus localized at the skyrmion core, as shown in Fig. 3 for the ground
state (energy $\varepsilon_{0}^{0}$ in Eq. (38)).

Figure 3:  Density plot of the ground-state wave function (see Eq. (38)) of an
electron, bound to a Belavin-Polyakov skyrmion of radius $R$, for a large
($R\gg\xi$) and intermediate-size ($R\sim\xi$) skyrmion. For $R\gg\xi$, the
wave function is localized well within the skyrmion core. It starts spreading
beyond the skyrmion core as the skyrmion radius $R$ decreases to become of the
order of the coherence length $\xi$, before spilling over $R$ for
$R\lesssim\xi$ (not shown here). For a large ($R\gtrsim\xi$) Néel skyrmion,
the low-lying bound states are those of an anisotropic harmonic oscillator.
They are spin-polarized essentially uniformly along $-\mathbf{e}_{x}$, as
shown in Fig. 1, as sketched by red arrows here. Along $y$, the characteristic
length scale of the ground-state wave function is $\bar{y}\sim
R\sqrt[4]{\frac{\xi}{R}}$, while the scale along $x$ is $\bar{x}\sim
R\sqrt[4]{\epsilon\frac{\xi}{R}}\ll\bar{y}$: it remains small compared with
$R$ even for $R\sim\xi$.

The bound states remain shallow, and the harmonic expansion continues to hold
for radia $R\gg\xi$. With $R$ approaching $\xi$, the frequency
$\hbar\omega_{y}$ becomes comparable to $\Delta$, and the low-energy
Hamiltonian (19) breaks down. That is, for radia $R\lesssim\xi$, the complete
Hamiltonian (9) must be considered.

At the same time, with $R$ decreasing further below $\xi$, the repulsive term
$\frac{\left(A_{y}^{z}\right)^{2}}{2m_{y}^{*}}\sim\Delta\frac{\xi^{2}}{R^{2}}$
in the low-energy Hamiltonian (19) grows relative to the attraction
$vA_{y}\sim\Delta\frac{\xi}{R}$. Using approximation
$A_{y}^{z}\sim\hbar\frac{x}{R^{2}}$ and estimating $\langle\hat{x}^{2}\rangle$
within harmonic approximation shows that the repulsion compensates the
attraction $vA_{y}^{x}$ for $R\sim\sqrt{\epsilon}\xi$. (Notice that
$\sqrt{\epsilon}\xi\gg a$, hence at this scale the continuum description of
the texture still applies.) Therefore, bound states become shallow again and
disappear at an $R=\bar{R}\sim\sqrt{\epsilon}\xi\ll\xi$. This estimate can be
obtained more rigorously Davier_2023 , and does not rely on strong mass
anisotropy. For a perfectly isotropic toy model, our calculation of $\bar{R}$
in Appendix D yields the same result.

We now turn to a study of deep bound states with the help of the full
Hamiltonian (9).

##  IV Belavin-Polyakov profile: Deep bound states

Above and in Ref. Davier_2023 , we have analysed shallow bound states at both
small $R\gtrsim\bar{R}$ and large BP skyrmion radia $R\gg\xi$. Now we turn to
intermediate radia $\bar{R}\lesssim R\lesssim\xi$, where the bound-state
energy becomes comparable with the gap $\Delta$.

###  IV.1 Particle-hole symmetry

Before turning to deep bound states, we would like to highlight the electron-
hole symmetry of Hamiltonian (9). If we go through steps (11) to (13), now
considering energies near the valence band maximum and setting
$E=-\Delta+\varepsilon$ with $|\varepsilon|\ll\Delta$, this time we obtain

$\Phi\simeq-\frac{\hat{\bm{\gamma}}_{-}}{2\Delta}\chi.$ (39)

The full 4-spinor wave function $(\phi,\chi)$ is now dominated by $\chi$. To
first order in $\frac{\hat{\bm{\gamma}}_{-}}{\Delta}\ll 1$, the effective
Hamiltonian for $\chi$ reads

$H_{-}=\hat{\bm{\gamma}}_{+}-\frac{\hat{\bm{\gamma}}_{-}^{2}}{2\Delta}.$ (40)

Comparing Eqs. (13) and (40), we see that shallow bound states near the top of
the valence band have polarization $|-\rangle_{\mathbf{e}_{\sigma}}$ and
energy $\varepsilon=E+\Delta>0$, while their conduction-band counterparts have
polarization $|+\rangle_{\mathbf{e}_{\sigma}}$ and energy
$\varepsilon=E-\Delta<0$. Both for large radia $R\gg\xi$ and for
$R\gtrsim\bar{R}$, shallow bound states near the bottom of the conduction band
and and their counterparts near the top of the valence bands obey the same
Schrödinger equation, and their energies are perfectly symmetric relative to
the gap center $E=0$. In other words, shallow bound states near the edge of
the conduction band are related to their valence-band partners by particle-
hole symmetry

$|+\rangle_{\mathbf{e}_{\sigma}}\to|-\rangle_{\mathbf{e}_{\sigma}},\hskip
28.45274pt\varepsilon\to-\varepsilon.$ (41)

In the intermediate range $R\sim\xi$, the full $4\times 4$ Hamiltonian (9) has
to be considered. Particle-hole symmetry (41) corresponds to wave function
transformation
$(\phi|+\rangle_{\mathbf{e}_{\sigma}},\chi|+\rangle_{\mathbf{e}_{\sigma}})^{t}\rightarrow(-\chi|-\rangle_{\mathbf{e}_{\sigma}},\phi|-\rangle_{\mathbf{e}_{\sigma}})^{t}$,
represented by operator

$M=\begin{pmatrix}0&-\sigma^{z}\\\ \sigma^{z}&0\end{pmatrix}.$ (42)

Applying this transformation to Hamiltonian (9) yields

$\tilde{H}=M^{\dagger}HM=\begin{pmatrix}\sigma^{z}(\hat{\bm{\gamma}}_{+}-\Delta)\sigma^{z}&-\sigma^{z}\hat{\bm{\gamma}}_{-}\sigma^{z}\\\
-\sigma^{z}\hat{\bm{\gamma}}_{-}\sigma^{z}&\sigma^{z}\left(\hat{\bm{\gamma}}_{+}+\Delta\right)\sigma^{z}\end{pmatrix}.$
(43)

Comparing Hamiltonians (43) and (9), we see that if

$\\{\hat{\bm{\gamma}}_{+},\sigma^{z}\\}=0,\hskip 14.22636pt\text{and}\hskip
14.22636pt\left[\hat{\bm{\gamma}}_{-},\sigma^{z}\right]=0$ (44)

(that is, if only linear terms are kept in the momentum expansion of
$\hat{\bm{\gamma}}_{\pm}$ in Eqs. (15)), then

$\tilde{H}=-H,$ (45)

which means exact particle-hole symmetry. Quadratic terms in the momentum
expansion break this symmetry, producing an asymmetry of the order of
$\Delta/\epsilon_{F}\ll 1$ relative to the otherwise particle-hole symmetric
background.

###  IV.2 Solving the full Hamiltonian

Now we turn to deep bound states in the intermediate range $R\sim\xi$ . In
this preliminary section, we focus on a $Q=1$ skyrmion, and keep the
presentation as general as possible. We consider Hamiltonian (9) and use the
expansion (15), only keeping its dominant terms, which is valid for
$R\sim\xi\gg a$. In this approximation, only differential operators acting on
$y$ are kept, since the ‘slow’ coordinate $x$ will only slightly affect the
energy; thus we fix it at its mean value which is, by symmetry, $x=0$.
Technically, this implies setting $A_{y}^{z}=0$ in $\hat{\bm{\gamma}}_{-}$,
which is further supported by the fact that the dominant potential term is
proportional to $\mathbf{e}_{\sigma}\cdot\bm{\sigma}$. With
$\mathbf{e}_{\sigma}$ assumed to be orthogonal to $\mathbf{e}_{z}$,
$A_{y}^{z}\sigma^{z}$ will only act in second order of perturbation theory.
Stated differently, by analyzing the electron motion only along the ‘fast’
coordinate $y$, we study the ‘gross’ structure of the bound-state energy
levels, temporarily ignoring the ‘fine’ structure: the sublevels, appearing
due to quantization of electron motion along the ‘slow’ coordinate $x$. The
approximations above imply

$\begin{split}\hat{\bm{\gamma}}_{+}&=vA_{y}^{\parallel}\mathbf{e}_{\sigma}\cdot\bm{\sigma}\equiv-\gamma_{+}\mathbf{e}_{\sigma}\cdot\bm{\sigma},\\\
\hat{\bm{\gamma}}_{-}&=v\hat{p}_{y}=\hat{\gamma}_{-}.\end{split}$ (46)

The Hamiltonian (9) can now be applied to bispinor
$(\phi_{1},\psi_{1},\psi_{2},\phi_{2})^{t}$ where the spinor components
$\phi_{i}$ and $\psi_{i}$ are given in the spin eigen basis
$|\pm\rangle_{\mathbf{e}_{\sigma}}$ along $\mathbf{e}_{\sigma}$.

The corresponding Schrödinger equation thus reads

$\begin{split}&(1)\hskip
22.76228pt(\Delta-\gamma_{+}-E)\phi_{1}+\hat{\gamma}_{-}\psi_{2}=0,\\\
&(2)\hskip
22.76228pt(\Delta+\gamma_{+}-E)\psi_{1}+\hat{\gamma}_{-}\phi_{2}=0,\\\
&(3)\hskip
22.76228pt\hat{\gamma}_{-}\phi_{1}-(\Delta+\gamma_{+}+E)\psi_{2}=0,\\\
&(4)\hskip
22.76228pt\hat{\gamma}_{-}\psi_{1}-(\Delta-\gamma_{+}+E)\phi_{2}=0.\end{split}$
(47)

The equations (2) and (3) allow us to express $\psi_{1}$ and $\psi_{2}$ as

$\begin{split}\psi_{1}&=-\frac{1}{\Delta+\gamma_{+}-E}\hat{\gamma}_{-}\phi_{2},\\\
\psi_{2}&=\frac{1}{\Delta+\gamma_{+}+E}\hat{\gamma}_{-}\phi_{1}.\end{split}$
(48)

Paying attention to non-commutativity of $\hat{\gamma}_{-}$ and $\gamma_{+}$,
we find

$\left[\hat{\gamma}_{-}^{2}-\frac{\hat{\gamma}_{-}(\gamma_{+})}{\Delta+\gamma_{+}\pm
E}\gamma_{-}+\Delta^{2}-\left(\gamma_{+}\pm E\right)^{2}\right]\phi=0.$ (49)

The particle-hole symmetry (41) is explicit here, with polarization inversion
$|+\rangle_{\mathbf{e}_{\sigma}}\to|-\rangle_{\mathbf{e}_{\sigma}}$
corresponding to replacing $\phi_{1}$ by $\phi_{2}$, with sign $+$ in equation
(49) associated with $\phi_{1}$.

This one-dimensional second-order differential equation can be easily solved
numerically. Searching for energies $E$, for which the wave function is
physical, leads to a discrete energy spectrum $E_{n}(R)$. These energies are
depicted as functions of the BP skyrmion radius $R$ in Figs. 4 and 5. Figure 4
shows numerical solutions for the ground-state energy and the energy of the
first excited state (with respect to the ‘fast’ degree of freedom $y$) in
comparison with the analytical result (38) obtained in the large-radius limit
$R\gg\xi$. The agreement is excellent for $R\gtrsim 8\xi$, as expected. Note
that the ground state is the only state whose energy crosses zero.

Figure 4: Ground-state energy (with polarization
$|+\rangle_{\mathbf{e}_{\sigma}}$) of Hamiltonian (9) as a function of
$\frac{R}{\xi}$, shown in units of $\Delta$. The numerical solution in blue is
compared with the analytical one (38) obtained with Harmonic approximation,
plotted in cyan and valid for $\frac{R}{\xi}\gg 1$. The energy of the first
excited state relative to the ‘fast’ degree of freedom is represented in dark
green, and again compared with its analytical counterpart (38) in green. The
abscissa of the crossing point between ground states energies $E_{0}^{c}$ and
$E_{0}^{v}$ is $R^{*}\simeq 0.37\xi$.

Fig. 5 shows numerical results for the first excited states. Their structure
is non-trivial for $R$ comparable to $\xi$, where the frequencies
$E_{c}^{i+1}-E_{c}^{i}$ are far from constant. The harmonic description holds
only for low-lying excited states, and for $R$ large against the coherence
length $\xi$.

Figure 5: Excited-state energies $E_{c}^{i}$ (in $\Delta$ units) as functions
of $\frac{R}{\xi}$, obtained by numerical solution of equation (49). Indices
$i$ represent the number of the wave-function zeros along $y$. The shown
states are arising from the conduction band, and have symmetric counterparts,
generated by the valence band.

Note that this solution describes only the $y$-dependence of the wave
function. As we already pointed out, it encapsulates the ‘gross’ structure of
bound-state spectrum. Each state of energy $E_{n}(R)$ is accompanied by a set
of excited states separated by frequencies of the order of $\hbar\omega_{x}$.
This is due to the fact that for the ‘slow’ degree of freedom $x$, the kinetic
term remains small compared with main potential term even in the range
$R\lesssim\xi$. In the picture above, the energy $E_{c}^{0}$ of the lowest
conduction-band bound state crosses zero at an $R\sim\xi$. Combined with the
particle-hole symmetry, this implies level crossing with bound states arising
from the valence band. We will now see if this level crossing is avoided.

###  IV.3 Perturbation-theory analysis of crossing points

Now we will treat the previously neglected terms in $\hat{\bm{\gamma}}_{\pm}$
of Eq. (15) as perturbation of the dominant terms (46). Figure 5 suggests that
the only states allowed to intersect are those whose wave function has no
zeros along the $y$ direction. Focusing on these states, we denote them as

$\Psi_{i}^{c}=\begin{pmatrix}\phi_{i}\\\ 0\\\ \psi_{i}\\\
0\end{pmatrix},\hskip 28.45274pt\Psi_{i}^{v}=\begin{pmatrix}0\\\ -\psi_{i}\\\
0\\\ \phi_{i}\end{pmatrix},$ (50)

where indices $i$ stand for the number of wave-function zeros along $x$. In
first approximation, these states can be presented as
$\phi_{i}(x,y)=\phi(y)\langle x|i\rangle$ and $\psi_{i}(x,y)=\psi(y)\langle
x|i\rangle$, with $|i\rangle$ the one-dimensional harmonic oscillator
eigenstates, and $\phi(y)$ the ground-state eigenfunction of Eq. (49). Here
$\psi=\frac{\hat{\gamma}^{0}_{-}}{\Delta+\gamma^{0}_{+}}\phi$, with
$\hat{\gamma}_{\pm}^{0}$ representing the zero-order expansion (46) of
$\left\|\hat{\bm{\gamma}}_{\pm}\right\|$.

An avoided crossing can only appear due to off-diagonal terms in the spin
basis defined by $\mathbf{e}_{\sigma}$. In $\hat{\bm{\gamma}}_{-}$, among such
terms the dominant one is $vA_{y}^{z}\sigma^{y}$ 444We make a change of basis
that amounts to a permutation of spin space indices $(x,y,z)$, replacing
$(x,y,z)$ by $(z,x,y)$.. This will split the crossing energies

$E_{i}^{0}=\langle\Psi_{i}^{c}|H_{0}|\Psi_{i}^{c}\rangle=\langle\Psi_{i}^{v}|H_{0}|\Psi_{i}^{v}\rangle$
(51)

with $H_{0}$ made of $\hat{\gamma}_{\pm}^{0}$, by the amount Landau1981Quantum

$\delta_{ij}=\langle\Psi_{i}^{c}|V|\Psi_{j}^{v}\rangle,\hskip
19.91684ptV=vA_{y}^{z}\sigma^{y}\otimes\sigma_{x},$ (52)

with $\sigma_{x}$ coupling spinors. To estimate the splitting, we approximate
$A_{y}^{z}\sim\Delta\frac{\xi}{R}\frac{\tilde{x}}{1+\tilde{y}^{2}}$, and
observe that the integral

$\int_{-\infty}^{+\infty}d\tilde{y}\,\frac{|\phi|^{2}+|\psi|^{2}}{1+\tilde{y}^{2}}$
(53)

appears to be numerically near unity for $R\sim\xi$. Using standard properties
of harmonic oscillator eigenstates, we find

$\delta_{ij}\sim\Delta\frac{\xi}{R}\langle
i|\hat{\tilde{x}}|j\rangle\sim\Delta\frac{\xi}{R}\left(\epsilon\frac{\xi}{R}\right)^{1/4}\left(\delta_{i,j+1}+\delta_{i+1,j}\right).$
(54)

In reality, the wave functions of crossing states are more complicated, and
$A_{y}^{z}$ is not a separable function. However, since $A_{y}^{z}$ is an odd
function of $x$, it couples states of different parities, meaning
$\delta_{ij}$ is expected to be non-zero if $|i-j|$ is an odd integer, and to
realize its maximum for $|i-j|=1$.

The other off-diagonal term is proportional to $A_{x}^{y}$, an odd function of
$y$, and is irrelevant since there is no crossing between eigenstates of
different parity along $y$, see Fig. 5.

These arguments suggest that, to second order in momentum expansion of Eq.
(15), there is no avoided crossing between levels of the same spatial
symmetry. Looking at the ratio

$\frac{\delta_{01}}{\hbar\omega_{x}}\simeq\frac{1}{2}\left[\epsilon\frac{\xi}{R}\left(1+\frac{\xi}{2R}\right)^{3}\right]^{-1/4},$
(55)

we see that it is larger than one as long as the skyrmion radius is larger
than

$R_{c}\sim\epsilon^{1/4}\xi\ll\xi$ (56)

with $\epsilon=\frac{m_{y}}{m_{x}}\frac{a}{\xi}\lll 1$. This appears to be the
radius for which the characteristic frequency of the ‘slow’ degree of freedom

$\hbar
w_{x}=\Delta\left(\frac{\xi}{R}\right)^{3/2}\sqrt{2\epsilon}\sqrt{1+\frac{\xi}{2R}}$
(57)

becomes of the order of $\Delta$. Thus, as long the present discussion makes
sense, the gaps opening at avoided crossing points are much larger that the
characteristic spacing between the energy levels. This observation implies
that the off-diagonal coupling term $A_{y}^{z}\sigma^{z}$ will avoid all
possible crossings between levels of the same spatial symmetry along the $x$
axis. Since the number of crossing levels is large, the precise energy
structure is difficult to identify. However, the qualitative picture of the
bound-state evolution as a function of the skyrmion radius $R$ is afforded by
the textbook argument for electron terms in a diatomic molecule
Landau1981Quantum  and is shown in Fig. 6. The argument states that, for a
single-parameter Hamiltonian, two electron terms of the same symmetry cannot
cross. In the present problem, the single parameter is the skyrmion radius
$R$, and the avoided crossing is ensured by the term
$V=vA_{y}^{z}\sigma^{y}\otimes\sigma_{x}$ 555This term only couples states of
opposite polarization and different spatial symmetry along $x$. Applying it as
a perturbation does however produce states of superposed polarizations,
meaning they will interact with states of any spatial symmetry due to other
terms in the Hamiltonian that have been previously ignored.

###  IV.4 Physics picture

The analysis above allows us to draw the following physics picture. For large
BP skyrmion radia $R\gg\xi$, the bound states remain shallow: that is, they
lie near the edge of the band they originated from: $|\varepsilon|\ll\Delta$.
For both $x$ and $y$ coordinates, the low-lying states are close to those of a
harmonic oscillator, see Eq. (38). As $R$ decreases to approach the coherence
length $\xi$, these low-lying states become more and more profound, and
approach the gap center at $R\sim R^{*}$. At lower radia $R\lesssim R^{*}$,
the energy gap is filled by a ladder of bound states of nearly harmonic
structure along $x$, with a characteristic level spacing $\hbar\omega_{x}$. As
$R$ decreases further, the harmonic structure along $x$ is lost, and levels
split to reach an energy difference of about $\Delta$ at $R\sim R_{c}$, before
the last bound state disappears at $R=\bar{R}$, with
$a\ll\bar{R}\sim\sqrt{\frac{m_{y}}{m_{x}}}\sqrt{a\xi}\ll\xi$. For $R<\bar{R}$,
there are no bound states any more, a result that holds as long as the
continuum approximation is valid, that is for $R\gg a$. Evolution of bound
states with $R$ is sketched in Fig. 6.

Figure 6:  Evolution of the bound-state levels as a function of the BP
skyrmion radius $R$, key features. $R^{*}$ represents the largest radius
allowing for level crossing before perturbation, see Fig. 4. And $R_{c}$ is
the radius for which the ‘slow’ frequency $\omega_{x}$ becomes comparable with
$\Delta$, see Eq. (56).

Note that, generally, the eigenstates of Hamiltonian (9) are spin-polarized
non-uniformly. In sections III and IV above, the Hamiltonian and the resulting
eigenstates were expressed in the basis obtained via transformations $P$ and
$U(\mathbf{r})$. The physical wave function must then be reconstructed from
the found solution $\Psi$ via

$\Psi_{\text{physical}}(\mathbf{r})=U(\mathbf{r})P\Psi(\mathbf{r}).$ (58)

Spin rotation $U(\mathbf{r})$ is non-uniform. Thus, generally, so is the spin
polarization of the wave function. However, for a large BP skyrmion
($R\gg\xi$), the wave functions of low-lying states are localized well within
the skyrmion core ($r\ll R$), where the Néel order parameter points along
$\hat{z}$, and thus $U(\mathbf{r})$ is close to $\pi$-rotation around
$\hat{z}$. Therefore, for $R\gg\xi$, the physical wave function is spin-
polarized opposite to $\hat{x}$, as shown in Figs. 3 and 8. However, for a
small skyrmion ($R\lesssim\xi$), spin polarization of the wave function
becomes substantially non-uniform even for low-lying bound states Davier_2023
.

Being parameterized by a single length scale and having a simple and elegant
analytic profile, Belavin-Polyakov skyrmion is a tractable example par
excellence, allowing one to analyze skyrmion-electron bound states. Moreover,
weak single-ion anisotropy affects the skyrmion shape only far away from the
center ($r\gg R$) Komineas_2020 ; Voronov_1983 , and thus does not influence
the low-lying bound states whose wave functions are localized well within the
skyrmion core. However, the BP profile is relevant in a parameter range where
the skyrmion lifetime appears to be short AF_Sk_Bessarab_2019 . Skyrmion-
electron bound states only increase the skyrmion energy, and thus reduce its
stability even further. We now turn to a different profile, which is described
by two length scales rather than one: that of a ‘domain wall’ (DW) skyrmion.
The DW profile becomes relevant near the transition to the modulated phase
Bogdanov_1989 , where skyrmions tend to have longer lifetimes
AF_Sk_Bessarab_2019 .

##  V Domain wall profile

In this section, we study the so-called Domain Wall (DW) skyrmion profile
defined by

$\theta(r)=\pi-2\arctan\left(\frac{\sinh(R/w)}{\sinh(r/w)}\right),$ (59)

with the skyrmion radius $R$ and the domain wall width $w$. The name ‘domain
wall skyrmion’ points to a nearly constant $\theta(r)\ll 1$ within the
skyrmion radius $R$, undergoing a smooth variation within the distance $w\ll
R$ around $r=R$, to become nearly constant ($\theta\approx\pi$) again at
$r>R$. That is, such a skyrmion can be seen as a circular domain wall of
radius $R$ and width $w$. We consider a DW skyrmion with winding number $Q=1$
and an arbitrary helicity $\gamma$. In the relevant limit $R\gg w$
komineas2020Large_radius , we can make an approximation

$\theta^{\prime}(r)\simeq\frac{1}{w\cosh\left(\frac{R-r}{w}\right)}\simeq\frac{\sin\theta}{w},$
(60)

which implies, near $r=R$, that

$\frac{\sin\theta}{r}\sim\frac{\sin\theta}{R}\ll\frac{\sin\theta}{w}\simeq\theta^{\prime}.$
(61)

The inequality above allows us to simplify the expression for the gauge field
in Eq. (33). For the dominant term, this leads to 666This is equivalent to the
inequality $\partial_{i}\theta\gg\partial_{i}\phi$, which allows us to keep
only the first term in the expression (29) for $\mathbf{A}_{i}$.

$v\mathbf{A}_{y}^{\parallel}\cdot\bm{\sigma}\simeq-\Delta\frac{\xi}{2w}\sin\theta\sin\alpha\,\mathbf{e}_{\sigma}^{DW}\cdot\bm{\sigma}$
(62)

with
$\mathbf{e}_{\sigma}^{DW}=(\sin(\alpha+\gamma),-\cos(\alpha+\gamma),0)^{t}$.
Helicity then plays the same role as for the BP profile: it only defines the
polarisation of the bound states. This potential has a non-trivial shape,
depicted in Fig. 7. It has two minima of depth $\Delta\frac{\xi}{2w}$ located
at points $\pm\mathbf{r}_{0}=\pm R\mathbf{e}_{y}$ with opposite polarization
relative to $\mathbf{e}_{\sigma}^{DW}$.

Figure 7: A sketch of the dominant potential term
$v\mathbf{A}_{y}^{\parallel}\cdot\mathbf{e}_{\sigma}^{DW}$, for a domain wall
profile of radius $R$, large compared with the wall width $w$.

###  V.1 Shallow bound states

We start by treating the limit $w\gg\xi$, where the potential is small against
$\Delta$. Direct inspection shows that, near the minima $\pm\mathbf{r}_{0}$,
the dominant repulsive term can be estimated as

$\begin{split}\frac{\left(A_{y}^{z}\right)^{2}}{2m_{y}^{*}}\sim\Delta\frac{\xi^{2}}{8R^{2}}x^{2}\ll
vA_{y}^{\|}\sim\Delta\frac{\xi}{2w}.\end{split}$ (63)

That is, repulsion remains small compared with the dominant attractive
potential $v\mathbf{A}_{y}^{\parallel}\cdot\bm{\sigma}$, and can be neglected
altogether. Cross-terms $\\{\hat{p}_{i},A_{i}^{z}\\}\sigma^{z}$ are also small
and off-diagonal in the basis of $\mathbf{e}_{\sigma}^{DW}$. Thus, in the
first approximation, they can also be ignored.

The potential can be diagonalized with the help of an operator

$P(\alpha)=\frac{1}{\sqrt{2}}\left(1-i\mathbf{e}_{r}\cdot\bm{\sigma}\right)$
(64)

with
$\mathbf{e}_{r}=\left(\cos(\gamma+\alpha),\sin(\gamma+\alpha),0\right)^{t}$,
which produces additional terms when sandwiching the operator
$\hat{p}_{i}^{2}\rightarrow P^{\dagger}(\alpha)\hat{p}_{i}^{2}P(\alpha)$
777This modifies the Peierls substitution in Eq. (6). . Nevertheless,
depending only on the polar angle $\alpha$, operator $P$ varies on the
characteristic scale $R$, and thus produces only terms of the order
$\Delta\frac{\xi^{2}}{R^{2}}$ (and even smaller for $\hat{p}_{x}$), negligible
compared with $vA_{y}^{\parallel}$.

The resulting Schrödinger equation for bound states near the point
$+\mathbf{r}_{0}$ is then

$\begin{split}\left[-\frac{\xi^{2}}{2w^{2}}\frac{\partial^{2}}{\partial\tilde{y}^{2}}-\frac{\epsilon}{2}\frac{\xi^{2}}{R^{2}}\frac{\partial^{2}}{\partial\tilde{x}^{2}}-\frac{\xi}{2w}\sin\theta\sin\alpha\right]\psi^{\uparrow}=\frac{\varepsilon}{\Delta}\psi^{\uparrow}\end{split}$
(65)

with $\tilde{x}=\frac{x}{R}$ and $\tilde{y}=\frac{y-R}{w}$, and the spin
projection $\uparrow$ along $\mathbf{e}_{\sigma}^{DW}$. Note that here the
anisotropy is much stronger than in the BP case: on top of the effective-mass
anisotropy, the potential itself is now strongly anisotropic. Such a strong
anisotropy again invites the Born-Oppenheimer approximation. The low-energy
limit $w\gg\xi$ is also appropriate for a harmonic expansion, since the
kinetic term is small compared with the depth of the potential well. The
potential can be expanded locally as

$-\sin(\alpha)\,\text{sech}\frac{r-R}{w}\simeq-1+\frac{\tilde{x}^{2}+\tilde{y}^{2}}{2}$
(66)

yielding the energy levels

$\frac{\varepsilon_{m}^{n}}{\Delta}=-\frac{\xi}{2w}\left[1-\left(n+\frac{1}{2}\right)\sqrt{\frac{2\xi}{w}}\right]+\left(m+\frac{1}{2}\right)\frac{\xi}{R}\sqrt{\frac{\epsilon}{2}\frac{\xi}{w}}.$
(67)

The spectrum above is similar to that of shallow bound states for the BP
profile in Eq. (38), except for an overall factor $\frac{1}{2}$: as a function
of $w$, the potential well is a half of what it was for the BP profile as a
function of $R$. The frequency of the ‘slow’ degree of freedom is however
smaller, due to the additional anisotropy of the potential. The wave functions
of these states are those of a two-dimensional anisotropic harmonic
oscillator. The key difference relative to the BP skyrmion is that, here, the
electron is predominantly localized at the two opposites edges of the DW
skyrmion rather than near its center, compare Figs. 3 and 8.

Figure 8: Density plot of the ground-state wave function (see Eq. (67)) of an
electron, bound to a domain-wall skyrmion in limits $w\gg\xi$ and $w\sim\xi$.
For a domain-wall skyrmion with a large wall width $w\gg\xi$, lowest electron
bound states lie deep in the effective potential well formed by the domain
wall. Their wave functions are thus close to those of a substantially
anisotropic two-dimensional harmonic oscillator, which implies strong
localization of the trapped electron inside the domain wall, as depicted in
the left panel of the figure. With the wall width $w$ decreasing to approach
$w\sim\xi$, the effective potential well becomes more and more shallow, which
makes the wave function of the low-lying electron states spill outside the
domain wall, as depicted in the right panel above. In both cases, the states
localized near the bottom well at $(0,\pm R)$ are polarized along
$\pm\mathbf{e}_{\sigma}^{DW}$, pointing along $-\mathbf{e}_{x}$ at the bottom
well. Away from $(0,\pm R)$, local spin polarization of the bound state
remains tangential to the domain wall, as shown in the right panel. The
characteristic spread of the ground state wave function along $y$ is
$\bar{y}\sim w\left(\frac{\xi}{w}\right)^{1/4}$, while along $x$ it is
$\bar{x}\sim R\left(\epsilon\frac{w\xi}{R^{2}}\right)^{1/4}$. Thus, by
contrast with $\bar{y}$, the bound-state range $\bar{x}$ along the domain wall
is sensitive not only to the domain wall width $w$, but also to the skyrmion
radius $R$.

Just as for the BP skyrmion, this low-energy description breaks down with $w$
decreasing to become comparable with $\xi$: at this point, the full
Hamiltonian (9) has to be solved.

###  V.2 Deep bound states

Equation (49) that we derived to treat the BP profile becomes even more
accurate for the DW configuration, since the ratio $|A_{y}^{z}|/|A_{y}^{\|}|$
is now smaller than for BP skyrmion. Again, we solve Eq. (49) numerically, and
find the ‘gross’ spectrum $E_{n}^{0}$ with $n$ the number of wave-function
zeros along $y$. The results are shown in Fig. 9 for the ground state, and
Fig. 10 for the excited states.

Figure 9: Ground state energy (in units of $\Delta$) given as a function of
the domain wall width $w$ in units of $\xi$. The numerical solution (in blue)
is compared with the analytical result (67), in red, obtained for $w\gg\xi$.
The agreement is excellent for $w\gtrsim 4\xi$, as indeed it should be.

Relative to the BP profile, the key difference can be observed already for the
ground state: it never crosses zero energy, even in the limit $w\to a$,
meaning there is no crossing between states generated by the valence and
conduction bands. Thus, the in-gap structure will be much simpler here than
for a BP skyrmion. The difference arises from the factor $\frac{1}{2}$ in
front of the square brackets in Eq. (67) as compared with Eq. (38): for a DW
skyrmion with a given $w$, the attractive potential is twice weaker than for
an equal-radius ($R=w$) BP skyrmion. The DW excited-state structure is very
similar to the BP case, even if the bound states are about twice as shallow.

Figure 10: Numerical results obtained from the resolution of equation (49)
with DW profile. The energies $E_{n}^{0}(R)$ (in units of $\Delta$) are given
as a function of wall with $w$ (in units of $\xi$). Subscript $n$ stands for
the number of wave-function zeros along $y$. Superscript 0 indicates that
corresponding states have no zeros along the $x$ coordinate.

The harmonic structure for sub-levels associated with the ‘slow’ degree of
freedom $x$ holds longer for the DW profile. Indeed, this time

$\frac{\hbar\omega_{x}}{\Delta}=\frac{\xi}{R}\sqrt{\frac{\epsilon}{2}\frac{\xi}{w}}$
(68)

becomes comparable with unity for

$w_{c}=\frac{m_{y}}{m_{x}}\frac{\xi^{2}}{2R^{2}}a.$ (69)

Thus, for $R\gtrsim\sqrt{\frac{m_{y}}{m_{x}}}\xi$, the energy
$\hbar\omega_{x}$ remains always small compared with the gap $\Delta$. Which
means that, in the continuum approximation, harmonic structure persists for
any radius $R$.

As in the BP case, the dominant repulsive term
$\frac{\left(A_{y}^{z}\right)^{2}}{2m_{y}^{*}}$ could be responsible for the
loss of bound states at a lower wall width $\bar{w}$. However, since harmonic
description holds for the slow degree of freedom, this term can be estimated
as

$\left\langle\frac{\left(A_{y}^{z}\right)^{2}}{2m_{y}^{*}}\right\rangle\sim\Delta\frac{\xi^{2}}{R^{2}}\langle\tilde{x}^{2}\rangle\sim\Delta\frac{\xi^{3}}{R^{3}}\sqrt{\epsilon\frac{w}{\xi}}.$
(70)

Therefore, this repulsion can never compensate the attraction
$vA_{y}^{\parallel}\sim\Delta\frac{\xi}{w}$. Which means that, for a DW
profile with $w\ll R$, bound states are always present, at least in the
continuum description.

For the DW profile, evolution of the bound-state spectrum as a function of $w$
looks similar to what one sees for the BP skyrmion as a function of $R$:
compare Figs. 11 and 6. At the same time, we observe two important
differences: (i) there is no critical DW width $\bar{w}$ below which bound
states disappear – at least, as long as $w\ll R$. And, since there is no
avoided crossing, (ii) polarisation of bound states remains pure, that is
$|+\rangle_{\mathbf{e}_{\sigma}}$ or $|-\rangle_{\mathbf{e}_{\sigma}}$. Note
that, in the leading approximation, the skyrmion radius $R$ has no impact on
the ‘gross’ structure of bound states: see the first term in Eq. (67).
However, $R$ does affect the ‘fine’ structure of bound sates, see
$\hbar\omega_{x}$ in Eq. (68).

Figure 11: The in-gap structure for a DW profile. There is no crossing between
states arising from valence- and conduction bands. The numerical solutions of
equation (49) hold until $w$ becomes comparable with critical width $w_{c}$,
where the fine-structure level spacing becomes comparable with $\Delta$, see
Eq. (69).

Having gained insight into the bound-state spectrum for both the BP and the DW
profiles, we now turn to the contribution of bound states to the skyrmion
energy at half-filling.

##  VI Energetics at half-filling

We consider a material with one electron per site, at zero temperature. In the
Néel phase, all the negative-energy states are occupied, while all the
positive-energy states are empty, see Figs. 6 and 11. Filled bound states,
generated by the valence band, rise above $E=-\Delta$ (the top of the valence
band in the absence of the skyrmion), thus making a positive contribution to
the skyrmion energy.

The electron contribution to the skyrmion energy may be defined as the
difference between the sums of occupied-state energies with and without the
skyrmion

$E_{T}=\sum_{i}\underbrace{(E^{v}_{i})^{\text{Bloch}}+(E_{i}^{v})^{\text{bound}}}_{\text{with\;skyrmion}}-\sum_{i}\underbrace{(E^{v}_{i})^{\text{Bloch}}}_{\text{without\;skyrmion}}.$

Calculation of $E_{T}$ requires understanding the effect of the skyrmion on
itinerant electron states. Let $N_{\text{bs}}(R)$ be the microscopic number of
valence-band bound states with polarization $|-\rangle_{\mathbf{e}_{\sigma}}$,
for a single valley $\Sigma_{i}$ (see fig 1). On the other hand, in Sec. II we
saw that in the uniform Néel state all the electron eigenstates are doubly
degenerate. Thus there are $N_{\text{bs}}$ states with the same unperturbed
energy, but polarization $|+\rangle_{\mathbf{e}_{\sigma}}$, that are repelled
from the skyrmion, resulting in a smaller density near the skyrmion center,
and a smaller energy. The intra-band energy shift for one of these states can
be estimated as

$\delta\varepsilon_{i}=v\int d^{2}rA_{y}^{\parallel}|u_{i}|^{2}\sim v\int
d^{2}rA_{y}^{\parallel}|u_{0}|^{2}$ (71)

with $u_{i}(\mathbf{r})$ the $i^{\text{th}}$ Bloch function. Since the
relevant states are close to the valence band maximum, we make the
approximation $u_{i}\sim u_{0}\sim\frac{1}{L}$ (with $L$ the characteristic
system size) in the spirit of the rigid band theorem Kittel1963 . This yields
the energy shift

$\delta\varepsilon\simeq\Delta\frac{\xi}{R}\frac{R^{2}}{L^{2}}\int_{0}^{L/R}\frac{z}{1+z^{2}}dz\simeq\Delta\frac{\xi}{R}\frac{R^{2}}{L^{2}}\ln\frac{L}{R},$
(72)

which is a factor of $\frac{R^{2}}{L^{2}}\ln\frac{L}{R}$ smaller than the
characteristic energy $\Delta\frac{\xi}{R}$ of the bound states. In the
thermodynamic limit $L\ggg R$, the characteristic energy of these
$N_{\text{bs}}$ $|+\rangle_{\mathbf{e}_{\sigma}}$-polarized states is thus
negligible relative to the characteristic energy $\Delta\frac{\xi}{R}$ of the
bound states. The impact of the skyrmion on the energy of Bloch states is
therefore tiny compared with the effect of bound state formation, since there
is an equal number $N_{\text{bs}}$ of states of each polarization. The
$N_{\text{bs}}$ $|+\rangle_{\mathbf{e}_{\sigma}}$-polarized states can then be
ignored when calculating the electron contribution $E_{T}$ to the skyrmion
energy.

Finally, consider the other states of the same valley, that are described by
the effective-mass Hamiltonian, but lie sufficiently deep in the band, so that
the corresponding $|-\rangle_{\mathbf{e}_{\sigma}}$-polarized states do not
become bound. These states will experience an energy shift similar to (72)
with a sign, dependent on their polarization. Since states of opposite
polarization are initially degenerate, the energy shifts of a pair of states
will tend to compensate. The resulting statement is that the total energy
$E_{T}$ can be approximated by the sum over occupied bound-state energies
only.

Before estimating $E_{T}$, a remark concerning electron-electron interaction.
While its careful treatment is beyond the scope of this work, for a
sufficiently large skyrmion electric neutrality will enforce uniform electron
density: in each valley, one spin polarization is attracted to the skyrmion
center to form bound states, while the opposite polarization is repelled and
forms itinerant states whose density is reduced near the skyrmion center. The
estimate of $E_{T}$ in the following subsection neglects electron-electron
interaction altogether. We now evaluate $E_{T}$ for the two skyrmion profiles
discussed above.

###  VI.1 Belavin-Polyakov skyrmion

####  VI.1.1 Harmonic states

We start by computing the contribution of harmonic states in the limit of
large skyrmion radius $R\gg\xi$. Treating the valence band maximum as energy
reference, we rewrite Eq. (38) as

$\varepsilon_{m}^{n}=\varepsilon_{0}^{0}-\left(\frac{1}{2}+n\right)\hbar\omega_{y}-\left(\frac{1}{2}+m\right)\hbar\omega_{x},$
(73)

with $\varepsilon_{0}^{0}=\Delta\frac{\xi}{R}$. Which allows us to write

$E_{T}=\sum_{m,n}\varepsilon_{m}^{n}=\sum_{n=0}^{N_{y}}\sum_{m=0}^{M_{n}}\left(\varepsilon_{0}^{n}-\left(\frac{1}{2}+m\right)\omega_{x}\right),\\\
$ (74)

with $N_{y}=\frac{\varepsilon_{0}^{0}}{\hbar\omega_{y}}$ and
$M_{n}=\frac{\varepsilon_{0}^{n}}{\hbar\omega_{x}}$. Since $N_{y}\gg 1$, this
yields

$E_{T}\simeq\frac{\left(\varepsilon_{0}^{0}\right)^{3}}{6\hbar^{2}\omega_{x}\omega_{y}}=\frac{\Delta}{6}\frac{(\xi/R)^{3}}{2\sqrt{\epsilon}(\xi/R)^{3}}\frac{1}{\sqrt{1+\xi/2R}},$
(75)

which implies asymptotic behaviour

$\lim_{R\to\infty}E_{T}(R)\simeq
E_{0},\;\;\text{with}\;\;E_{0}\equiv\frac{\Delta}{12\sqrt{\epsilon}}.$ (76)

Notice that, at large $R$, the increasing number of bound states and their
decreasing characteristic energy compensate each other. Also notice that
$E_{0}\gg\Delta$ emerges as a characteristic scale, defining the contribution
of bound states to the skyrmion energy. Note that the estimate above treats
bound states as harmonic (equidistant) and underestimates their number, thus
giving a lower bound of the bound-state contribution $E_{T}$ to the skyrmion
energy.

We now evaluate $E_{T}$ for intermediate radia $R_{c}\ll R\lesssim\xi$.

####  VI.1.2 Contribution of deep bound states

In the following section, we focus on the density of in-gap states since, in
the first approximation, the sum over energies does not depend on the fine
structure of energy levels. This amounts to ignoring the perturbation of
harmonic structure, discussed in subsection IV.3. In this context, the physics
picture is roughly the one presented in Fig. 6. Following the steps of Chap.
IV and assuming that, in the intermediate regime $R\gtrsim R^{*}$, the
harmonic description holds only for the $x$ direction, this yields

$E_{T}\simeq\frac{1}{2\hbar\omega_{x}}\sum_{n}(\varepsilon_{n}^{0})^{2}$ (77)

which can be evaluated numerically, as shown in Fig. 12.

Figure 12: Numerical estimate of the bound-state contribution to the skyrmion
energy $E_{T}$, shown in units of $E_{0}=\frac{\Delta}{12\sqrt{\epsilon}}$.
Orange line corresponds to Eq. (77), where the sum is calculated numerically
for the first seven excited states, the maximum of the 7${}^{\text{th}}$ level
energy being already smaller than $4.10^{-3}\Delta$ for this radius range. The
apparent asymptotic value, close to $1.34E_{0}$, is of the same order of
magnitude but slightly larger than the one obtained with harmonic
approximation (76), as expected. The blue part is calculated for radius
smaller than the crossing radia $R<R^{*}$ using Eq. (79).

As the skyrmion radius decreases below $R^{*}$, the bound states produced by
valence and conduction bands begin to hybridize, complicating the problem.
However, for $R\gg R_{c}$, the ‘fine-structure’ levels produced by
quantization of motion along $x$ remain roughly harmonic (equidistant). That
is, each state of the ‘gross’ structure, labelled relative to the $y$
coordinate, can be seen as having a ‘fine’-structure band with a density of
states

$\nu_{0}(R)=\frac{1}{\hbar\omega_{x}}=\frac{E_{0}}{\sqrt{2}\Delta^{2}}\left(\frac{R}{\xi}\right)^{3/2}\left(1+\frac{\xi}{2R}\right)^{-1/2}$
(78)

attached to it. In this picture, for $R\geq\bar{R}_{1}$ (with $\bar{R}_{1}$
being the radius where the first excited bound state vanishes) the band
$\varepsilon_{1}^{v}$ associated with first excited state can be considered
completely filled. The bands $\varepsilon_{0}^{v}$ and $\varepsilon_{0}^{c}$
are filled up to energy $\Delta$, see Fig. 13. The resulting total energy cost
can thus be estimated as

$\begin{split}E_{T}&\sim\frac{\nu}{2}\left(\varepsilon_{1}^{0}\right)^{2}+\nu\int_{0}^{\Delta}\varepsilon
d\varepsilon+\nu\int_{2\Delta-\varepsilon_{0}^{0}}^{\Delta}\varepsilon
d\varepsilon\\\
&\sim\frac{\nu}{2}\left(\varepsilon_{1}^{0}\right)^{2}+\nu\Delta^{2}\left(-1+2\frac{\varepsilon_{0}^{0}}{\Delta}-\frac{1}{2}\left(\frac{\varepsilon_{0}^{0}}{\Delta}\right)^{2}\right).\end{split}$
(79)

For $R<R_{1}$ this result holds by simply taking $\varepsilon_{1}^{0}=0$, the
corresponding plot is given in Fig. 12.

Figure 13: Schematic representation of occupied states at half filling.
Occupied “bands” are coloured. The top of the valence band is taken as zero
energy, as for calculations made in this chapter. The bound state spectrum
remaining approximately equidistant for $R\sim R^{*}$ allows us to approximate
the density of states as being constant as a function of energy inside the
gap, and equal to $1/\hbar\omega_{x}$ in region II. In the region I, where
valence-band states $\varepsilon_{v}^{0}$ mix with conduction-band states
$\varepsilon_{c}^{0}$, the density of states is expected to be about
$2/\hbar\omega_{x}$. Note that the length $R_{c}$ is not shown here, since it
is vanishingly small compared with the coherence length $\xi$, as is the loss
radius $\bar{R}$.

For $R<R_{c}$, the numerical results in Fig. 12 can no longer be trusted.
However, at $R<\bar{R}$, there are no bound states any more, and thus the
bound-state contribution to the skyrmion energy goes to zero at $R\to\bar{R}$.
Therefore, to minimize the bound-state contribution to its energy, a BP
skyrmion shall shrink to a radius comparable to $\bar{R}$. Note that these
results are only qualitative, but the overall picture of the energy cost
$E_{T}(R)$ should be captured by the above derivation. We now turn to
estimating the bound-state contribution to the energy of a DW skyrmion.

###  VI.2 Domain-wall skyrmion

Figure 14: The bound-state contribution to the DW skyrmion energy $E_{T}$,
shown in units of $E_{0}\frac{R}{\xi}$ as a function of the wall width $w$.
The asymptotic behaviour is $E_{T}\simeq 1.26E_{0}\frac{R}{w}$, with a
coefficient slightly greater that the factor $1$ predicted within harmonic
approximation, as expected.

The shallow bound-state limit $w\gg\xi$ can be treated as for the BP profile,
with the help of Eqs. (67) and (75), yielding

$E_{T}\underset{w\gg\xi}{\sim}2\frac{\varepsilon_{0}^{3}}{6\hbar^{2}\omega_{x}\omega_{y}}=E_{0}\frac{R}{w},$
(80)

where the factor of 2 above accounts for the two potential wells generated by
the skyrmion. The bound-state contribution $E_{T}$ to the skyrmion energy thus
decreases with $w$ increasing relative to the DW skyrmion radius $R$, but only
for $w\ll R$. Meaning that (for a given $R$) $E_{T}$ has a lower bound of the
order of $E_{0}$. For the DW profile, there is no crossing of bound states
that emerge from the conduction- and valence bands. Therefore, the energy cost
for $w\lesssim\xi$ can be estimated as we did it for the BP profile for
$R>R^{*}$ (see Eq. (77)), leading here to

$E_{T}\simeq\frac{2}{2\hbar\omega_{x}}\sum_{n}(\varepsilon_{n}^{0})^{2}=12\sqrt{2}E_{0}\frac{R}{\xi}\sqrt{\frac{w}{\xi}}\sum_{n}\frac{(\varepsilon_{n}^{0})^{2}}{\Delta^{2}}.$
(81)

Estimating the above sum numerically yields Fig. 14. The energy cost $E_{T}$
tends to zero as $w\rightarrow 0$ and, naturally, grows linearly with the DW
skyrmion radius $R$.

In a way, the above behaviour of $E_{T}(R,w)$ favors the skyrmion shrinking to
reduce both its radius $R$ and the wall width $w$. In other words, the bound-
state contribution favours compact profiles similar to the BP one, producing a
single potential well.

###  VI.3 Straight domain wall

Figure 15: Straight domain wall at an angle $\alpha$ to the $x$-axis in real
space, the axis defined for a given valley as shown in Fig. 1. Unit normal
$\mathbf{u}$ to the wall points at the same angle $\alpha$ with respect to the
$x$ axis. Polar angle $\theta$ varies along $\mathbf{u}$ as a function of $s$
in Eq. (82). Angle $\alpha$ governs the bound-state formation: the attractive
potential $v(\mathbf{A}_{y}\cdot\bm{\sigma})$ of Eq. (84) is the strongest for
$\alpha=\pi/2$, that is for $\theta$ dependent only on $y$ while the wall
stretches along $x$. By contrast, at $\alpha=0$ attraction ceases to exist:
$v(\mathbf{A}_{y}\cdot\bm{\sigma})\equiv 0$ for $\theta$ dependent only on $x$
while the wall stretches along $y$.

Now we consider an electron in the presence of a straight domain wall. The
problem is interesting for several reasons: firstly, in its own right.
Secondly – because it is quasi-exactly solvable. And, last but not least –
since it allows one to gain insight into the origin of the peculiar form of
the effective potential in Eq. (62) (see Fig. 7) as well as into the effect of
electron bound states on the transition to the modulated phase Bogdanov_1989 .

Orientation of a straight domain wall in real space is defined by an in-plane
unit vector $\mathbf{u}=(\cos\alpha,\sin\alpha,0)^{t}$, orthogonal to the
wall, see Fig. 15, while its profile is described by

$\theta(s)=2\arctan\left(e^{-s}\right),\hskip
19.91684pts\equiv\frac{\cos(\alpha)\,x+\sin(\alpha)\,y}{w}.$ (82)

According to (29), it produces a simple gauge field

$\mathbf{A}_{i}\cdot\bm{\sigma}=\frac{1}{2}(\mathbf{e}_{\sigma}\cdot\bm{\sigma})\partial_{i}\theta\,$
(83)

where $\mathbf{e}_{\sigma}=(-\sin\phi,\cos\phi,0)^{t}$ with
$\phi=\alpha+\gamma$. Parameter $\gamma$ is similar to the helicity of a
skyrmion, and defines bound-state polarization. The only potential term reads

$v(\mathbf{A}_{y}\cdot\bm{\sigma})=(\mathbf{e}_{\sigma}\cdot\bm{\sigma})\frac{\Delta}{2}\frac{\xi}{w}\sin\alpha\,\text{sech}\,s\,,$
(84)

it vanishes for $\sin\alpha=0$. This is hardly surprising: the presence of
$\mathbf{A}_{y}=-i\hbar U^{\dagger}(\partial_{y}U)$ in the l.h.s. above means
that, to produce attraction and bound states for electrons in the given
valley, the Néel order parameter $\mathbf{n}$ must depend on $y$. Viewing a
large-radius DW skyrmion as a circle, and approximating it as a regular
inscribed polygon with straight domain-wall segments, we see from Eq. (84)
that the attraction is the strongest on segments that are nearly normal to the
$y$ axis – indeed, as we already saw in Fig. 7.

The limit $w\gg\xi$ can again be treated with effective Hamiltonian (19) in
the harmonic approximation, yielding the low-lying spectrum

$\begin{split}\varepsilon_{k}^{n}=&\;-\Delta\frac{\xi}{w}\sin\alpha+\left(\frac{1}{2}+n\right)\frac{\hbar}{w}\sqrt{\frac{\Delta\xi\sin\alpha}{m_{s}w}}\\\
&+\left(\frac{1}{2m_{q}}-\frac{m_{s}}{8M^{2}}\right)\frac{\hbar^{2}k^{2}}{w^{2}}.\end{split}$
(85)

We introduced here effective masses

$\begin{split}\frac{1}{m_{s}}=\frac{\cos^{2}\alpha}{m_{x}}+\frac{\sin^{2}\alpha}{m_{y}^{*}},\hskip
14.22636pt\frac{1}{m_{q}}=\frac{\sin^{2}\alpha}{m_{x}}+\frac{\cos^{2}\alpha}{m_{y}^{*}},\\\
\frac{1}{M}=\sin 2\alpha\left(\frac{1}{m_{y}^{*}}-\frac{1}{m_{x}}\right)\hskip
42.67912pt\end{split}$ (86)

along pertinent directions, and the wave vector $k$ along the domain wall. The
energy cost of a domain wall in this limit can now be evaluated as

$E_{T}=-\sum_{n}\int_{\varepsilon_{0}^{n}}^{0}\nu(\varepsilon)\varepsilon
d\varepsilon\simeq
E_{DW}\sum_{n}\left(\frac{|\varepsilon_{0}^{n}|}{\Delta}\right)^{3/2}$ (87)

with $\nu$ the DOS associated with delocalized states along the wall, and
$E_{DW}=\frac{16\sqrt{2}}{\pi}E_{0}\frac{L}{\xi}\sin\alpha$ the characteristic
energy of a DW. A simple calculation shows that for $w\gg\xi$ this expression
tends to

$E_{T}\simeq\frac{32\sqrt{2}}{5\pi}\sin^{2}\alpha E_{0}\frac{L}{w},$ (88)

showing that the energy cost grows linearly with the domain wall length $L$,
indeed as expected.

For $w\lesssim\xi$, Eq. (49) can be applied directly: note that this time
there is no need to neglect the term $vA_{y}^{z}\sigma^{z}$ in
$\hat{\bm{\gamma}}_{-}$ since now it is identically equal to zero. Eq. (49) is
equivalent to the one for a DW skyrmion, up to a redefinition of the coherence
length $\xi\to\sin\alpha\,\xi$. The r.h.s. of Eq. (87) can then be evaluated
using the energy derived numerically for a circular skyrmion, the result is
given in units of $E_{DW}$ in Fig. 16. The limit $w\geq\xi$ shows a
$\left(\frac{w}{\sin(\alpha)\xi}\right)^{-1}$ behaviour as expected. Note that
the total energy cost of a DW always varies linearly with its length $L$, as
does $E_{DW}$. In the region $w\ll\xi$, there is a single non-zero term in the
sum (87) that tends toward 1, corresponding to a unique occupied bound-state
band. The bound-state contribution to the energy of the domain wall is then
expected to converge toward $E_{DW}$ as $w$ goes to zero.

Figure 16: Numerical estimation of the energy $E_{T}$ of a straight DW,
performing the sum (87). The result is given in units of $E_{DW}$, the natural
energy unit for a straight DW. The asymptotic behaviour for $w\geq\xi$ is in
$\left(\frac{w}{\sin(\alpha)\xi}\right)^{-1}$, see dashed red curve, as
expected from the harmonic approximation (88). The limit of $E_{T}$ for
$w\ll\xi$ is finite and precisely equal to $E_{D}$, see main text. The dashed
green curve give the estimation of the sum (87), by fitting linearly
$\varepsilon_{0}^{0}(w)$.

Not surprisingly, for the wall orientation $\alpha=0$, the $E_{T}$ in Eq. (88)
vanishes, simply because the Néel order parameter does not vary along $y$,
hence no attractive potential and no bound states. However, this is true only
for the valley $\Sigma_{i}$ that we considered to define $\alpha$: for the
next valley, $\Sigma_{i+1}$, one has $\alpha_{i+1}=\alpha_{i}+\frac{\pi}{2}$.
That is, taking into account contributions from all the four $\Sigma$-valleys,
and considering the $w\gg\xi$ limit, $E_{T}$ would contain two terms with
angle dependence $\sin^{2}\alpha$ and two terms with
$\sin^{2}\left(\alpha+\frac{\pi}{2}\right)=\cos^{2}\alpha$, which yields

$E_{T}\simeq\frac{64\sqrt{2}}{5\pi}E_{0}\frac{L}{w},$ (89)

meaning that the bound-state contribution to the domain-wall energy from all
the four $\Sigma$ valleys does not depend on the orientation of the wall.

Equation (89) and its counterpart for a circular DW skyrmion, Eq. (80), imply
that skyrmion-electron bound states strongly influence the transition to the
modulated phase Bogdanov_1989 . To see this, recall that the bound-state
contribution to the skyrmion energy complements the purely magnetic part
$E[\mathbf{n}]$ of the energy of a texture Bogdanov_1989 :

$E[\mathbf{n}]=\int
d^{2}r\left[J\left(\bm{\nabla}\mathbf{n}\right)^{2}+\frac{K}{a_{0}^{2}}(1-n_{z}^{2})+\frac{D}{a_{0}}\epsilon_{DM}\right],$
(90)

where $J$ is stiffness, $K$ the easy-axis anisotropy, $D$ the
Dzyaloshinskii–Moriya coupling, and the Dzyaloshinskii–Moriya energy density
$\epsilon_{DM}$ depends linearly on $\nabla\mathbf{n}$. For a domain wall of
length $L$ and width $w$, one has $|\bm{\nabla}\mathbf{n}|\sim 1/w$ within the
range of about $w$ from the center of the wall and zero otherwise, while
$(1-n_{z}^{2})\sim 1$ within the same range of about $w$, and zero elsewhere.
With this estimate in mind, naïve integration of the different terms in Eq.
(90) along the wall and transversely to it yields the following characteristic
dependence of $E[\mathbf{n}]$ on $L$ and $w$:

$E[\mathbf{n}]\sim
L\left[\frac{J}{w}+\frac{K}{a_{0}^{2}}w-\frac{D}{a_{0}}\right].$ (91)

Eq. (91) is written up to numerical coefficients of the order of unity,
dependent on the details of the domain wall profile, just as the coefficients
in Eqs. (89) and (80). These coefficients do not affect the conclusion below.
The first two terms in Eq. (91) define the optimal domain wall width
$w_{0}=a_{0}\sqrt{\frac{J}{K}}$ that minimizes $E[\mathbf{n}]$, thus producing
the characteristic dependence of the domain wall energy on $J$, $K$ and $D$:

$E(w_{0})\sim\frac{L}{a_{0}}\left[\sqrt{JK}-D\right].$ (92)

A more accurate calculation yields the energy $E(w_{0})$ becoming negative at

$D_{c}=\frac{4}{\pi}\sqrt{JK},$ (93)

marking the transition to a modulated state Bogdanov_1989 .

Now, comparison of Eq. (91) with Eqs. (89) and (80) shows that the bound-state
contribution $E_{T}$ to the domain-wall energy renormalizes and, in fact,
makes a dominant contribution to stiffness $J$ in Eqs. (91) and (93) as per

$J\rightarrow\tilde{J}=J+\frac{\Delta}{\sqrt{\epsilon}}=J+\sqrt{\frac{m_{x}}{m_{y}}}\sqrt{\varepsilon_{F}\Delta}.$
(94)

The significant upward renormalization of $J$ in Eq. (94) increases the domain
wall width. Even more importantly, it substantially increases $D_{c}[J,K]$ in
Eq. (93) while keeping its $\sqrt{K}$ scaling, and thus extends the stability
range of the uniform state, as shown in Fig. 17. Note that this picture holds
for $w\gtrsim\xi$, that is for sufficiently small values of
$K\lesssim\tilde{J}\left[\frac{a_{0}}{\xi}\right]^{2}$.

Figure 17: The effect of skyrmion-electron bound states on the phase diagram:
the critical value $D_{c}[J,K]$, marking the transition between the uniform
and the modulated phases, increases from $D_{c}=\frac{4}{\pi}\sqrt{JK}$ to
$\tilde{D}_{c}=\frac{4}{\pi}\sqrt{\tilde{J}K}$ as per Eq. (94), while keeping
its $\sqrt{K}$ dependence. Notice that this picture holds for $w\gtrsim\xi$,
which limits it to small values of
$K\lesssim\tilde{J}\left[\frac{a_{0}}{\xi}\right]^{2}$.

Numerically, how large a renormalization may arise from Eq. (94)? As a
concrete example, consider optimally-doped high-temperature superconductor
NCCO (Nd1.85Ce0.15CuO4). The presence of Néel order in this compound remains
controversial Luke_1990 ; Mang_2004 ; Motoyama_2007 ; Saadaoui_2015 , yet a
number of experiments Yamada_2003 ; Kang_2005 ; Yu_2007 ; Dorantes_2018 ;
RR_Zeeman_SOC_2021  report static or quasi-static Néel order. The electron
subsystem of NCCO is nearly two-dimensional, with the low-temperature ratio of
the $c$-axis and $ab$-plane resistivities of about $10^{3}$ to $10^{4}$
Armitage_2010 . The two-dimensional band structure of NCCO involves small
carrier pockets at points $\Sigma$, with anisotropic dispersion as in Fig. 1.
Assuming quadratic dispersion
$\frac{p_{x}^{2}}{2m_{x}}+\frac{p_{y}^{2}}{2m_{y}}$, the mass anisotropy can
be read off ARPES maps in Refs Armitage_2002 ; He_2019  as $m_{x}/m_{y}\sim
25$. The $\varepsilon_{F}/\Delta\sim\xi/a$ ratio can be evaluated, knowing the
NCCO lattice constant $a\approx 3.95$Å Armitage_2010 , and estimating $\xi\sim
9$ nm – for instance, from magnetic quantum oscillation data
RR_Zeeman_SOC_2021 . As a result, for $\tilde{D}_{c}/D_{c}$ in Fig. 17, with
numerical parameters of NCCO we find

$\frac{\tilde{D}_{c}}{D_{c}}\sim\sqrt{\frac{\tilde{J}}{J}}\sim\sqrt[4]{\frac{m_{x}}{m_{y}}\frac{\varepsilon_{F}}{\Delta}}\sim
5.$

The bound-state contribution to the energy of a domain-wall texture in Eqs.
(89) and (80) proves large compared with the ‘magnetic’ contribution of Eq.
(90). This is in contrast with what one finds for electron states, bound to a
vortex in a type-II superconductor DeGennes_Book ; tinkham2004 : in the latter
problem, the gradient contribution to the vortex energy has the form
$\left(\frac{\phi_{0}}{4\pi\lambda}\right)^{2}\ln\frac{\lambda}{\xi}$, where
$\xi$ is the superconducting coherence length, and $\lambda\gg\xi$ the
penetration depth. The electron states, bound to the vortex core, only correct
the large argument of the logarithm by a coefficient of the order of unity
DeGennes_Book ; tinkham2004 , thus making a parametrically small contribution
to the energy of the vortex.

##  VII Discussion

###  VII.1 Validity of static-skyrmion approximation

Above, all the textures were implicitly treated as static input to the
electron problem. This assumption must be taken with caution, since an
isolated antiferromagnetic skyrmion possesses its own excitation modes
Kravchuk_2019  that couple to the bound electron. To justify viewing the
skyrmion as static, these eigenmodes must have frequencies far below those of
electron bound states. The proper skyrmion frequencies are closely bound from
above Kravchuk_2019  by the gap $\hbar\Omega_{0}$ in the magnon spectrum,
which is about $\hbar\Omega_{0}\sim\sqrt{2JK}$ Rezende2019 ; Kittel_1951  for
systems where single-ion anisotropy $K$ is much smaller than stiffness $J$.
Thus, it is $\Omega_{0}$ that must be compared with characteristic electron
frequencies, whose scales are different for the ‘gross’ and ‘fine’ structures
of bound-electron levels.

For the BP profile, the ‘gross’ structure is defined by the frequency
$\hbar\omega_{y}\sim\Delta\left(\frac{\xi}{R}\right)^{3/2}$ in the first term
of Eq. (38). Inequality $\omega_{y}\gg\Omega_{0}$ limits the validity of the
static-texture approximation to not-too-large skyrmions

$R\ll R_{y}^{BP}=\xi\left(\frac{\Delta^{2}}{JK}\right)^{1/3}.$

At the same time, consistency with the effective-Hamiltonian approach requires
that the bound states remain shallow. According to Eq. (38), this implies
$R\gg\xi$. Therefore, the upper bound $R_{y}^{BP}$ above must be large
compared with $\xi$, which thus demands Davier_2023

$\frac{K}{J}\ll\left(\frac{\Delta}{J}\right)^{2}.$ (95)

The r.h.s. of the inequality (95) may be of the order of unity or smaller:
recall that, for a lattice model, $\Delta$ scales linearly with the lattice
spin $S$ (ordered moment), whereas $J\sim S^{2}$. Nevertheless, for a not-too-
large $S$, inequality (95) simply requires that $K/J$ be small against unity.

For the ‘fine’ structure of the bound-state spectrum, the condition to meet
for the BP profile is more stringent: now it is the ‘low’ frequency
$\omega_{x}$ in the second term of Eq. (38) that must be large relative to
$\Omega_{0}$:

$\hbar\omega_{x}\sim\Delta\sqrt{\epsilon}\left(\frac{\xi}{R}\right)^{3/2}\gg\sqrt{2JK}.$
(96)

This limits the validity of the static-texture approach to

$R\ll R_{x}^{BP}=\epsilon^{1/3}\xi\left(\frac{\Delta^{2}}{JK}\right)^{1/3}.$

Again, for consistency with the effective-mass approximation, $R_{x}^{BP}$
must be large against $\xi$, which thus demands

$\frac{K}{J}\ll\epsilon\left(\frac{\Delta}{J}\right)^{2}.$ (97)

Inequality (97) is similar to (95), but has an extra factor of
$\epsilon=\frac{m_{y}}{m_{x}}\frac{a}{\xi}\ll 1$ in the right-hand side.

Numerically, how stringent are inequalities (95) and (97)? As in the
immediately preceding section, we turn to NCCO as a concrete example. The
values used in Subsection VI.3 yield $\epsilon\sim 2\cdot 10^{-3}$. As a
result, the ‘gross’ structure of electron states, bound to a not-too-large BP
skyrmion, can be treated in the static-skyrmion approximation for the $K/J$
ratio sufficiently small compared with unity. By contrast, treating the fine
structure of levels in the same approximation requires a much smaller
anisotropy $K/J\ll 10^{-3}$.

For a domain-wall (DW) profile, the wall width $w$ is defined by anistropy as
per $w\sim a_{0}\sqrt{\tilde{J}/K}$, see Sec. VI.3. Just as for the BP
skyrmion, one constraint on the validity of our results stems from Eq. (67)
and from the requirement that the bound states remain shallow, which implies
$w\gg\xi$. The latter inequality can be recast as

$\frac{K}{J}\ll\frac{\tilde{J}}{J}\left(\frac{a_{0}}{\xi}\right)^{2}.$

Using Eq. (94) and treating the second term in its r.h.s. as dominant, one
finds

$\frac{K}{J}\ll\frac{\Delta}{J}\sqrt{\frac{m_{x}}{m_{y}}}\left(\frac{a_{0}}{\xi}\right)^{3/2}.$
(98)

At the same time, the validity of the static-skyrmion approximation is limited
by the condition that the characteristic electron frequency be large compared
with that of skyrmion eigenmodes.

For the ‘gross’ structure of electron levels this means that the scale
$\Delta(\xi/w)^{3/2}$ of the ‘high’ oscillator frequency in the first term of
Eq. (67) must be large compared with the gap $\hbar\Omega_{0}\sim\sqrt{2JK}$
in the magnon spectrum. Taking into account the characteristic domain wall
width $w\sim a_{0}\sqrt{\tilde{J}/K}$, this limits the single-ion anisotropy
$K$ from below rather than from above:

$\frac{K}{J}\gg\frac{J}{\Delta}\left(\frac{a_{0}}{\xi}\right)^{9/2}\left(\frac{m_{x}}{m_{y}}\right)^{3/2}.$
(99)

Together, the inequalities (98) and (99) read

$\frac{J}{\Delta}\left(\frac{a_{0}}{\xi}\right)^{9/2}\left(\frac{m_{x}}{m_{y}}\right)^{3/2}\ll\frac{K}{J}\ll\frac{\Delta}{J}\sqrt{\frac{m_{x}}{m_{y}}}\left(\frac{a_{0}}{\xi}\right)^{3/2}.$
(100)

What numerical bounds can this produce? Again, for a rough estimate, consider
NCCO as a concrete (even if not necessarily universal) example. The material
parameters above yield

$10^{-4}\ll\frac{K}{J}\ll 5\cdot 10^{-2}.$

Notice that the upper-to-lower bound ratio in inequality (100) is

$\left(\frac{\Delta}{J}\right)^{2}\frac{m_{y}}{m_{x}}\left(\frac{\xi}{a_{0}}\right)^{3}\gg
1.$

It contains (i) the $\left(\frac{\Delta}{J}\right)^{2}$ factor that is likely
of the order of unity, (ii) the $\frac{m_{y}}{m_{x}}$ ratio that tends to be
small albeit not parametrically so (about 0.04 for NCCO) and, finally, (iii)
the factor $\left(\frac{\xi}{a_{0}}\right)^{3}$ that is indeed parametrically
large (over $10^{4}$ for NCCO) and overwhelms all the rest.

To treat the ‘fine’ structure of electron levels in the static-skyrmion
approximation, one has to demand that the ‘low’ frequency in the second term
of Eq. (67) be large compared with the gap $\hbar\Omega_{0}\sim\sqrt{2JK}$ in
the magnon spectrum:

$\Delta\sqrt{\epsilon}\frac{\xi^{3/2}}{w^{1/2}R}\gg\sqrt{JK}.$ (101)

With $w\sim a_{0}\sqrt{\tilde{J}/K}$, this inequality translates into another
upper bound, but far more stringent than (98):

$\frac{K}{J}\ll\left(\frac{m_{y}}{m_{x}}\right)^{5/2}\left(\frac{a_{0}}{\xi}\right)^{1/2}\left(\frac{\xi}{R}\right)^{4}\left(\frac{\Delta}{J}\right)^{3}.$
(102)

Again assuming $\Delta/J\sim 1$, for NCCO the product of the first two factors
above falls below $10^{-4}$. Now, recall that $\xi\ll w\ll R$. Even for
$R/\xi\sim 10$, the factor $\left(\xi/R\right)^{4}$ above would then yield
$K/J\ll 10^{-8}$. That is, realistically, for a DW skyrmion, fine structure of
bound electron states can be treated in the static-skyrmion approximation only
for perfectly isotropic magnets.

While the constraints of this section limit the validity of the static-
skyrmion approximation as applied to shallow skyrmion-electron bound states,
note that the appearance of texture-bound electron states is a much broader
phenomenon, extending to regimes where electron degrees of freedom and those
of a texture must be treated on an equal footing.

###  VII.2 Experimental signatures

Combined with the appearance of skyrmion-electron bound states, doping a half-
filled antiferromagnetic insulator turns a skyrmion into a charged particle.
The latter can be detected and manipulated with the help of electric field.

Even exactly at half-filling, resonance transitions between different bound-
state levels would have a distinct frequency spectrum. The latter may involve
transitions at $\hbar\omega<2\Delta$ (see Figs. 6 and 11), which could serve
as a fingerprint of skyrmion-electron bound states. Away from half-filling,
resonance transitions between the bound-state levels at frequencies
$\hbar\omega\ll\Delta$ could serve as another spectroscopic hallmark. A study
of such signatures will be the subject of forthcoming work.

###  VII.3 Conclusion

We examined a smooth texture in a Néel antiferromagnet with a specific
location of the electron band extrema, and showed that it produces a peculiar
spin-orbit coupling that has no analogue in a ferromagnet. The coupling has
the scale $\hbar v/L$, where $v$ is the Fermi velocity of the underlying
paramagnetic state, and $L$ is the relevant characteristic length scale of the
texture. For topological textures such as skyrmions and domain walls, this
coupling generates bound electron states in the gap of the electron spectrum.

Motivated by fundamental and technological interest alike, we focused on an
interesting limit, where the proper frequencies of the skyrmion are small
against those of the bound electron. In this limit, bound-state energies are
sensitive to the skyrmion ‘shape’; the latter serves as a static input to the
electron problem. We studied evolution of bound states with variation of the
skyrmion profile for the Belavin-Polyakov and domain-wall skyrmions. Bound
states prove to be energetically costly, and this (a) favors smaller BP
skyrmions that produce no bound states, and (b) shifts the transition line
between the uniform and the modulated states to substantially larger values of
the Dzyaloshinskiy-Moriya coupling.

With a major research effort focusing on magnetic topological textures, it is
vital to understand their interplay with other degrees of freedom of the
crystal, most notably with band electrons. We hope that our results contribute
to a better understanding of magnetic textures and to putting them to use.

##  VIII Acknowledgement

It is our pleasure to thank G. E. Volovik and V. Geshkenbein for enlightening
discussions.

##  Appendix A Comparison with the problem of a texture in a ferromagnet

It is instructive to compare the problem we studied above with its
ferromagnetic counterpart: band electron in the presence of a texture in a
ferromagnet. Smoothly varying as a function of coordinate $\bf{r}$, local
magnetization $\bf{M_{r}}$ couples to the band electron spin $\bm{\sigma}$ via
exchange term $J(\bf{M_{r}}\cdot\bm{\sigma})$ (cf. Subsection II.1). For the
ease of comparison, we introduce notation
$\Delta(\bf{n_{r}}\cdot\bm{\sigma})\equiv J(\bf{M_{r}}\cdot\bm{\sigma})$, with
unit vector $\bf{n_{r}}$ pointing in the direction of local magnetization
$\bf{M_{r}}$.

Figure 18: Band dispersion in the uniform ferromagnetic state, where a spin-
degenerate band of the parent paramagnetic phase is Stoner-split into two non-
degenerate subbands. Panel (a) shows the limit $\Delta\ll\varepsilon_{F}$,
where the minima of both subbands lie deep at the bottom of the Fermi sea, and
thus are inaccessible to low-energy and low-temperature experimental probes.
Panel (b) shows the limit $2\Delta\gtrsim\varepsilon_{F}$ at half-filling,
where the lower subband is completely filled while the upper subband is empty.
Lightly doping such a state renders the physics of electron states near the
upper band minimum accessible to low-energy measurements. Such states can be
described by the effective Hamiltonian $\mathcal{H}_{eff}$ of Eq. 104.

In the uniform state, we choose ${\bf n_{r}}={\bf n}_{0}$ to point along the
$z$ axis. The band electron spectrum is thus simply

$\mathcal{H}=\varepsilon({\bf p})+\Delta\sigma_{z},$

where $\varepsilon({\bf p})$ is the notional band dispersion in the absence of
ferromagnetism, and the second term above describes a uniform Stoner splitting
White2007  of $2\Delta$ between the spin-up and spin-down subbands with
energies $\varepsilon_{\pm}({\bf p})=\varepsilon({\bf p})\pm\Delta$, see Fig.
18.

In the presence of a smooth non-uniform texture $\bf{n_{r}}$, it is convenient
to employ the unitary transformation of Eq. (5) and recast the Hamiltonian as

$\mathcal{H}=\varepsilon\left(\hat{\bf
p}+\mathbf{A}\cdot\bm{\sigma}\right)+\Delta\sigma_{z},$ (103)

see Eqs. (6) and (7). As a result, the Stoner splitting term is rendered
uniform, and now the texture influences the electron motion via the gauge
potential $\mathbf{A}\cdot\bm{\sigma}$.

A smooth texture induces a perturbation on an energy scale that is small
against both $\Delta$ and the Fermi energy $\varepsilon_{F}$. Since we are
interested in possible appearance of bound states, let us notice that such
states become visible to low-temperature measurements only for
$\Delta\gtrsim\varepsilon_{F}$ – see Fig. 18. It is thus fitting to compare
our results for shallow bound states in an antiferromagnet with a toy model of
a ferromagnet with an entirely filled lower subband of width
$\varepsilon_{F}\lesssim 2\Delta$, as shown in Fig. 18(b): at low
temperatures, a single dopant electron will find itself at low-lying levels of
the otherwise empty upper subband, which will correspond to the setting that
we studied above for an antiferromagnet.

As in Section II.2, we are now in a position to write down an effective-mass
Hamiltonian that will describe low-lying electron states near the bottom of
the upper subband in Fig. 18(b). To this end, notice that spin-$x$ and -$y$
components of the vector potential $\mathbf{A}\cdot\bm{\sigma}$ produce matrix
elements only between the spin-up and spin-down subbands, but not within
either of the two. As a result, they can contribute to the effective-mass
Hamiltonian only via perturbation theory, producing terms of the order of
$\frac{\hbar^{2}}{mL^{2}\Delta}\ll 1$ relative to those of the first term in
Eq. (103). Therefore, to zeroth order in $\hbar^{2}/(mL^{2}\Delta)$, the
upper-subband ($\sigma_{z}=+1$) effective-mass Hamiltonian reads

$\mathcal{H}_{\text{eff}}=\frac{1}{2m}\left(\hat{p}_{i}+A_{i}^{z}\right)^{2},$
(104)

where, without loss of generality, we assume isotropic effective mass and set
zero energy at the bottom of the upper subband.

Hamiltonian (104) is precisely what one would obtain in an antiferromagnet
upon omitting all but the first term in Hamiltonians (16) and (19), and
keeping to the spin-up problem. Now we will show that, for a $Q=1$ Belavin-
Polyakov skyrmion, this Hamiltonian does not produce skyrmion-electron bound
states, contrary to Hamiltonians (16) and (19). To this end, consider the
problem for a Belavin-Polyakov skyrmion in the gauge of Eq. (27), that gives
rise to $A_{i}^{z}$ of Eq. (32). Given isotropic mass above, it is convenient
to switch to polar coordinates $r=\sqrt{x^{2}+y^{2}}$ and
$\alpha=\arctan\frac{y}{x}$ as in Subsection III.1. We also switch to polar
components of $A_{i}^{z}$, that is from $(A_{x}^{z},A_{y}^{z})$ to
$(A_{r}^{z},A_{\alpha}^{z})$, where
$A_{r}^{z}=(\textbf{e}_{r}\cdot\textbf{A}^{z})$ and
$A_{\alpha}^{z}=(\textbf{e}_{\alpha}\cdot\textbf{A}^{z})$, while
$\textbf{e}_{r}=(\cos\alpha,\sin\alpha)$ and
$\textbf{e}_{\alpha}=(-\sin\alpha,\cos\alpha)$. Using Eq. (32), we find

$A_{r}^{z}=0,\,\,A_{\alpha}^{z}=\frac{\hbar r}{R^{2}+r^{2}}$ (105)

Hamiltonian (104) now reads

$\displaystyle\mathcal{H}_{\text{eff}}$ $\displaystyle=$
$\displaystyle\frac{\hbar^{2}}{2m}\left[-\frac{1}{r}\frac{\partial}{\partial
r}r\frac{\partial}{\partial
r}+\left(\frac{l}{r}+A_{\alpha}^{z}\right)^{2}\right]=$ (106) $\displaystyle=$
$\displaystyle\frac{\hbar^{2}}{2m}\left[-\frac{1}{r}\frac{\partial}{\partial
r}r\frac{\partial}{\partial
r}+\left(\frac{l}{r}+\frac{r}{R^{2}+r^{2}}\right)^{2}\right]$ (107)

Here, we took into account $\mathcal{H}_{\text{eff}}$ being independent of
$\alpha$, and replaced the orbital momentum operator
$\hat{l}=-i\frac{\partial}{\partial\alpha}$ by its integer eigenvalue $l=0,\pm
1,\pm 2,..$. We see that $A_{\alpha}^{z}$ only modifies the centrifugal energy
term while keeping it positive-definite. The first term is the kinetic energy
of a free particle in an $l=0$ state. Hence no bound states, certainly not for
the $Q=1$ skyrmion above: it is indeed the spin-orbit coupling term $v{\bf
A}_{y}^{\|}\cdot\bm{\sigma}$ in Hamiltonians (16) and (19) that generates
skyrmion-electron bound states in a Néel antiferromagnet.

##  Appendix B Gauge transformation

Operator $W$ can be written as

$W=e^{i\phi\mathbf{n}\cdot\bm{\sigma}}=\cos\phi+i\sin\phi\,\mathbf{n}\cdot\bm{\sigma},$
(108)

leading to

$W^{\dagger}\partial_{i}W=i\partial_{i}\phi\,\mathbf{n}\cdot\bm{\sigma}+i\sin\phi\,W^{\dagger}\partial_{i}\mathbf{n}\cdot\bm{\sigma}$
(109)

and then to

$\begin{split}&\tilde{\mathbf{A}_{i}}=-i\hbar
U^{\dagger}W^{\dagger}\partial_{i}\left(WU\right)=\mathbf{A}_{i}+\hbar\partial_{i}\phi\sigma^{z}\\\
&+\hbar U^{\dagger}\left[\frac{\sin
2\phi}{2}\partial_{i}\mathbf{n}\cdot\bm{\sigma}+\sin^{2}\phi(\mathbf{n}\times\partial_{i}\mathbf{n})\cdot\bm{\sigma}\right]U.\end{split}$
(110)

The last term lies in the plane orthogonal to $\mathbf{e}_{z}$, and thus only
contributes to the transformation of $\mathbf{A}_{i}^{\|}$.

Acting with $\partial_{i}$ on Eq. (5) yields

$\partial_{i}U^{\dagger}\mathbf{n}\cdot\bm{\sigma}U+U^{\dagger}\partial_{i}\mathbf{n}\cdot\bm{\sigma}U+U^{\dagger}\mathbf{n}\cdot\bm{\sigma}\partial_{i}U=0.$
(111)

A right multiplication by $U^{\dagger}\mathbf{n}\cdot\bm{\sigma}U=\sigma^{z}$
then leads to

$\partial_{i}U^{\dagger}U-iU^{\dagger}(\mathbf{n}\times\partial_{i}\mathbf{n})\cdot\bm{\sigma}U+\sigma^{z}U^{\dagger}\partial_{i}U\sigma^{z}=0.$
(112)

Finally, the identity $\partial_{i}\left(U^{\dagger}U\right)=0$ allows to
write

$U^{\dagger}\partial_{i}U-\sigma^{z}U^{\dagger}\partial_{i}U\sigma^{z}=-iU^{\dagger}(\mathbf{n}\times\partial_{i}\mathbf{n})\cdot\bm{\sigma}U.$
(113)

Which, by definition of gauge field $A_{i}^{\alpha}$ in Eq. (6), means that
its in-plane components reads

$\mathbf{A}_{i}^{\parallel}=-\frac{\hbar}{2}U^{\dagger}\left[(\mathbf{n}\times\partial_{i}\mathbf{n})\cdot\bm{\sigma}\right]U.$
(114)

Together with (110), this implies

$\begin{split}\tilde{\mathbf{A}}_{i}^{\|}=\hbar U^{\dagger}\left[\frac{\sin
2\phi}{2}\partial_{i}\mathbf{n}\cdot\bm{\sigma}\right.\hskip 85.35826pt\\\
+\left.\left(\sin^{2}\phi-\frac{1}{2}\right)(\mathbf{n}\times\partial_{i}\mathbf{n})\cdot\bm{\sigma}\right]U.\end{split}$
(115)

Using the identity $\sin^{2}\phi-\frac{1}{2}=-\frac{1}{2}\cos 2\phi$, we then
find

$\left(\tilde{\mathbf{A}}_{i}^{\|}\right)^{2}=\frac{\hbar^{2}}{4}||\partial_{i}\mathbf{n}||^{2}=\left(\mathbf{A}_{i}^{\|}\right)^{2}.$
(116)

The norm $\left(\mathbf{A}_{i}^{\|}\right)^{2}$ being conserved, the in plane
components $\mathbf{A}_{i}^{\|}$ thus only undergo an in-plane rotation, as
they do under gauge transformation by operator $V_{z}$.

From the gauge properties above, we can deduce the allowed terms

$\displaystyle\left(\mathbf{A}_{i}^{\parallel}\right)^{n}$ $\displaystyle\to$
$\displaystyle\left(\tilde{\mathbf{A}}_{i}^{\parallel}\right)^{n}$
$\displaystyle\left(\hat{p}_{i}+A_{i}^{z}\sigma^{z}\right)^{n}$
$\displaystyle\to$
$\displaystyle\left(\hat{p}_{i}+\tilde{A}_{i}^{z}\sigma^{z}\right)^{n}$ (117)
$\displaystyle\\{\hat{p}_{i},\mathbf{A}_{i}^{\parallel}\cdot\bm{\sigma}\\}$
$\displaystyle\to$ $\displaystyle\\{\hat{p}_{i}+\delta
A_{i}^{z}\sigma^{z},\tilde{\mathbf{A}}_{i}^{\parallel}\cdot\bm{\sigma}\\}=\\{\hat{p}_{i},\tilde{\mathbf{A}}_{i}^{\parallel}\cdot\bm{\sigma}\\}$

where in the last line we used the identity $\\{\sigma^{i},\sigma^{j}\\}=0$
for $i\neq j$. We observe that all the terms present in
$\hat{\bm{\gamma}}_{\pm}$ in Eqs. (15) are thus allowed, meaning that
Hamiltonian (9) and its low-energy limit (19) are indeed gauge-invariant.

##  Appendix C Gauge invariants

In this Appendix, we set $\hbar=1$. In Subsection II.4, we saw that gauge
transformations (LABEL:transformation_de_jauge) and (23) amount only to local
rotation of in-plane components $\mathbf{A}_{i}^{\parallel}$, thus conserving
$\left(\mathbf{A}_{i}^{\parallel}\right)^{2}$ as well as
$\left[\mathbf{A}_{x}^{\parallel}\times\mathbf{A}_{y}^{\parallel}\right]$ and
$\left(\mathbf{A}_{x}^{\parallel}\cdot\mathbf{A}_{y}^{\parallel}\right)$. The
absolute value
$\left(\mathbf{A}^{\parallel}\right)^{2}=\sum_{i=x,y}\left(\mathbf{A}_{i}^{\parallel}\right)^{2}$
can be expressed as (with implicit sum over repeated indices)

$4\left(\mathbf{A}^{\parallel}\right)^{2}=\left(\partial_{i}\theta\partial^{i}\theta+\sin^{2}\theta\,\partial_{i}\phi\partial^{i}\phi\right)=\left(\bm{\nabla}\cdot\mathbf{n}\right)^{2},$
(118)

that is the gradient energy density of the isotropic antiferromagnet. Vector
product
$\left[\mathbf{A}_{x}^{\parallel}\times\mathbf{A}_{y}^{\parallel}\right]$ can
be written as

$\begin{split}\mathbf{A}_{x}^{\parallel}\times\mathbf{A}_{y}^{\parallel}&=\frac{1}{4}\sin\theta\left(\partial_{x}\theta\partial_{y}\phi-\partial_{y}\theta\partial_{x}\phi\right)\\\
&=\frac{1}{4}\left(\mathbf{n}\cdot\left[\partial_{x}\mathbf{n}\times\partial_{y}\mathbf{n}\right]\right)=\pi\mathcal{D}_{s}(\mathbf{r}),\end{split}$
(119)

where $\mathcal{D}_{s}$ is the skyrmion density Tatara_2008 ; Jiang2017 ,
defined so that

$\int d^{2}r\mathcal{D}_{s}(\mathbf{r})=\int
d^{2}r\frac{1}{4\pi}\left(\partial_{x}\mathbf{n}\times\partial_{y}\mathbf{n}\right)\cdot\mathbf{n}=\mathcal{N}_{s}$
(120)

where $\mathcal{N}_{s}\in\mathbb{Z}$ is the total topological charge of
skyrmions present in the plane. Another invariant, suggested by the
transformation law of $A_{i}^{z}$ in Eq. (LABEL:transformation_de_jauge) is
the fully antisymmetric tensor

$F^{z}_{ij}=\partial_{i}A_{j}^{z}-\partial_{j}A_{i}^{z}.$ (121)

The latter invariant can be recast as

$\begin{split}F_{ij}^{z}&=\partial_{i}\left(\partial_{j}\phi\frac{1-\cos\theta}{2}\right)-\partial_{j}\left(\partial_{i}\phi\frac{1-\cos\theta}{2}\right)\\\
&=\frac{1}{2}\sin\theta\left(\partial_{i}\theta\partial_{j}\phi-\partial_{j}\theta\partial_{i}\phi\right)=2\pi\epsilon_{ij}\mathcal{D}_{s}(\mathbf{r}),\\\
\end{split}$ (122)

with $\epsilon_{ij}$ the fully antisymmetric tensor. This means that there are
only two distinct invariants associated with gauge transformation
(LABEL:transformation_de_jauge). Note that it also implies that the textbook
non-abelian gauge field tensor is identically equal to zero:

$\mathbf{F}_{ij}\equiv\partial_{i}\mathbf{A}_{j}-\partial_{j}\mathbf{A}_{i}+i\left[\mathbf{A}_{i},\mathbf{A}_{j}\right]=0.$
(123)

This is due to the ”electromagnetic” part exactly cancelling the anti-
commutator, both being proportional to skyrmion density.

##  Appendix D Isotropic toy model

In this appendix, we show that the loss of bound states at small BP skyrmion
radia is not related to the mass anisotropy in Eq. (17). A perfectly isotropic
toy model allows us to capture the scale $\bar{R}=\sqrt{\epsilon}\xi$, below
which the bound states are lost. To this end, consider Hamiltonian (19) with
isotropic mass

$H_{i}=\sum_{i=x,y}\frac{1}{2m^{*}}\left(\hat{p}_{i}+A_{i}^{z}\sigma^{z}\right)^{2}+vA_{y}^{x}\sigma_{x}.$
(124)

The repulsive potential

$\frac{1}{2m^{*}}\sum_{i}\left(A_{i}^{z}\right)^{2}=\frac{\hbar^{2}}{2m^{*}}\frac{(1-\cos\theta)^{2}}{4r^{2}}$
(125)

does not depend on polar angle $\alpha$, while the cross term is proportional
to

$\sum_{i}A_{i}^{z}\partial_{i}=\frac{(1-\cos\theta)}{2r}\mathbf{e}_{\alpha}\cdot\bm{\nabla}=\frac{(1-\cos\theta)}{2r^{2}}\frac{\partial}{\partial\alpha}$
(126)

and thus only acts on $\alpha$. That is, Hamiltonian (124) allows separation
of variables, and the angular momentum $l$ is a good quantum number. For
$l=0$, Hamiltonian (124) reads

$\begin{split}H_{i}^{l=0}=&\;\frac{\hbar^{2}}{2m^{*}R^{2}}\left(-\frac{1}{z}\frac{\partial}{\partial
z}\left(z\frac{\partial}{\partial
z}\right)+\frac{z^{2}}{(1+z^{2})^{2}}\right)\\\ &-\frac{\hbar
v}{R}\frac{1}{1+z^{2}}\sigma^{x}.\end{split}$ (127)

To keep the same density of states as for Hamiltonian (19), we choose the
effective mass $m^{*}=\sqrt{m_{x}^{*}m_{y}^{*}}$. Now treat Hamiltonian (127),
considering the potential as a perturbation relative to the kinetic energy
Landau1981Quantum . To this end, we find the wave function in two separate
regions. First, denote the sum of the potential terms in Eq. (127) as $U(z)$,
and consider the skyrmion core, where $|\varepsilon|\ll|U(z)|$: the
Schrödinger equation reads

$\frac{\hbar^{2}}{2m^{*}R^{2}}\frac{1}{z}\frac{d}{dz}\left(z\frac{d\psi}{dz}\right)=U(z)\psi(z).$
(128)

The above equation can now be integrated along $z$ to the upper limit $z_{1}$
of the attractive region:

$\left.\frac{d\psi}{dz}\right|_{z_{1}}\simeq\frac{1}{z_{1}}\frac{2m^{*}R^{2}}{\hbar^{2}}\psi(z_{1})\int_{0}^{z_{1}}U(z)zdz.$
(129)

Here we took advantage of the observation that the sought wave function is
spread out, and does not vary much within the potential well. By contrast, far
from the attractive region, i.e. for $z\gg 1$, the Schrödinger equation
describes a free particle, whose wave function is proportional to the Hankel
function $H_{0}^{(1)}(i\kappa z)$, with
$\kappa=\sqrt{\frac{2m^{*}R^{2}|\varepsilon|}{\hbar^{2}}}\ll 1$. For $\kappa
z\ll 1$, this function behaves asymptotically as $\ln(\kappa z)$. The value of
$z_{1}$ being of the order of unity and thus $\kappa z_{1}$ being small allows
us to match the logarithmic derivative $\frac{1}{\psi}\frac{d\psi}{dz}$ of
this asymptotic form to that in Eq. (129) at $z=z_{1}$, which yields

$\begin{split}\varepsilon&\simeq-\frac{\hbar^{2}}{2m^{*}R^{2}z_{1}^{2}}\exp\left[\frac{\hbar^{2}}{m^{*}R^{2}}\left(\int_{0}^{z_{1}}U(z)zdz\right)^{-1}\right]\\\
&\simeq-\Delta\frac{R_{0}}{R}\frac{\xi}{R}\frac{e^{2/\ln
z_{1}}}{z_{1}^{2}}\exp\left[\frac{1}{1-\frac{R}{R_{0}}}\right]\end{split}$
(130)

with $R_{0}\equiv\frac{\sqrt{\epsilon}}{2}\xi$. The expression above
demonstrates the existence of a skyrmion radius $\bar{R}\sim R_{0}$, where the
lowest bound state merges into the continuum and disappears. The same result
can be obtained by solving Hamiltonian (127) numerically.

## References

  * (1) B. Göbel, I. Mertig, and O. A. Tretiakov, “Beyond skyrmions: Review and perspectives of alternative magnetic quasiparticles,” _Physics Reports_ , vol. 895, p. 1–28, Feb 2021. [Online]. Available: http://dx.doi.org/10.1016/j.physrep.2020.10.001
  * (2) W. Jiang, G. Chen, K. Liu, J. Zang, S. G. te Velthuis, and A. Hoffmann, “Skyrmions in magnetic multilayers,” _Physics Reports_ , vol. 704, pp. 1–49, 2017, skyrmions in Magnetic Multilayers. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0370157317302934
  * (3) C. Back, V. Cros, H. Ebert, K. Everschor-Sitte, A. Fert, M. Garst, T. Ma, S. Mankovsky, T. L. Monchesky, M. Mostovoy, N. Nagaosa, S. S. P. Parkin, C. Pfleiderer, N. Reyren, A. Rosch, Y. Taguchi, Y. Tokura, K. von Bergmann, and J. Zang, “The 2020 skyrmionics roadmap,” _Journal of Physics D: Applied Physics_ , vol. 53, no. 36, p. 363001, jun 2020. [Online]. Available: https://doi.org/10.1088/1361-6463/ab8418
  * (4) S. S. P. Parkin, M. Hayashi, and L. Thomas, “Magnetic domain-wall racetrack memory,” _Science_ , vol. 320, no. 5873, pp. 190–194, 2008. [Online]. Available: https://www.science.org/doi/abs/10.1126/science.1145799
  * (5) S. Parkin and S.-H. Yang, “Memory on the racetrack,” _Nature Nanotechnology_ , vol. 10, no. 3, pp. 195–198, 2015. [Online]. Available: https://doi.org/10.1038/nnano.2015.41
  * (6) A. Fert, V. Cros, and J. Sampaio, “Skyrmions on the track,” _Nature Nanotechnology_ , vol. 8, no. 3, pp. 152–156, 2013. [Online]. Available: https://doi.org/10.1038/nnano.2013.29
  * (7) X. Wang, H. Yuan, and X. Wang, “A theory on skyrmion size,” _Communications Physics_ , vol. 1, 07 2018.
  * (8) P. F. Bessarab, D. Yudin, D. R. Gulevich, P. Wadley, M. Titov, and O. A. Tretiakov, “Stability and lifetime of antiferromagnetic skyrmions,” _Physical Review B_ , vol. 99, no. 14, p. 140411(R), Apr 2019. [Online]. Available: http://dx.doi.org/10.1103/PhysRevB.99.140411
  * (9) X. Zhang, Y. Zhou, and M. Ezawa, “Antiferromagnetic skyrmion: Stability, creation and manipulation,” _Scientific Reports_ , vol. 6, no. 1, p. 24795, Apr 2016. [Online]. Available: http://dx.doi.org/10.1038/srep24795
  * (10) B. Göbel, A. Mook, J. Henk, and I. Mertig, “Antiferromagnetic skyrmion crystals: Generation, topological hall, and topological spin hall effect,” _Phys. Rev. B_ , vol. 96, p. 060406(R), Aug 2017. [Online]. Available: https://link.aps.org/doi/10.1103/PhysRevB.96.060406
  * (11) W. Legrand, D. Maccariello, F. Ajejas, S. Collin, A. Vecchiola, K. Bouzehouane, N. Reyren, V. Cros, and A. Fert, “Room-temperature stabilization of antiferromagnetic skyrmions in synthetic antiferromagnets,” _Nature Materials_ , vol. 19, no. 1, pp. 34–42, 2020. [Online]. Available: https://doi.org/10.1038/s41563-019-0468-3
  * (12) J. Barker and O. A. Tretiakov, “Static and dynamical properties of antiferromagnetic skyrmions in the presence of applied current and temperature,” _Phys. Rev. Lett._ , vol. 116, p. 147203, Apr 2016. [Online]. Available: https://link.aps.org/doi/10.1103/PhysRevLett.116.147203
  * (13) A. N. Bogdanov, U. K. Rößler, M. Wolf, and K.-H. Müller, “Magnetic structures and reorientation transitions in noncentrosymmetric uniaxial antiferromagnets,” _Phys. Rev. B_ , vol. 66, p. 214410, Dec 2002. [Online]. Available: https://link.aps.org/doi/10.1103/PhysRevB.66.214410
  * (14) R. Cheng and Q. Niu, “Electron dynamics in slowly varying antiferromagnetic texture,” _Phys. Rev. B_ , vol. 86, p. 245118, Dec 2012. [Online]. Available: https://link.aps.org/doi/10.1103/PhysRevB.86.245118
  * (15) N. Davier and R. Ramazashvili, “Texture-induced spin-orbit coupling and skyrmion-electron bound states in a néel antiferromagnet,” _Phys. Rev. B_ , vol. 107, p. 014406, Jan 2023. [Online]. Available: https://link.aps.org/doi/10.1103/PhysRevB.107.014406
  * (16) A. A. Belavin and A. M. Polyakov, “Metastable states of two-dimensional isotropic ferromagnets,” _Soviet Journal of Experimental and Theoretical Physics Letters_ , vol. 22, p. 245, Nov. 1975. [Online]. Available: http://www.jetpletters.ru/ps/1529/article_23383.shtml
  * (17) R. Rajaraman, _An Introduction to Solitons and Instantons in Quantum Field Theory_. North-Holland Publishing Company, 1982.
  * (18) A. N. Bogdanov and D. A. Yablonskiǐ, “Thermodynamically stable “vortices” in magnetically ordered crystals. the mixed state of magnets,” _Sov. Phys. JETP_ , vol. 68, p. 101, 01 1989. [Online]. Available: http://jetp.ras.ru/cgi-bin/e/index/e/68/1/p101?a=list
  * (19) N. Kulikov and V. V. Tugushev, “Reviews of topical problems: Spin-density waves and itinerant antiferromagnetism in metals,” _Physics-Uspekhi_ , vol. 27, pp. 954–976, 1984. [Online]. Available: https://ufn.ru/en/articles/1984/12/c/
  * (20) G. E. Volovik, “Linear momentum in ferromagnets,” _Journal of Physics C: Solid State Physics_ , vol. 20, no. 7, pp. L83–L87, mar 1987. [Online]. Available: https://doi.org/10.1088/0022-3719/20/7/003
  * (21) C. Kittel, _Quantum theory of solids_. New York: Wiley, 1963.
  * (22) L. Landau, E. Lifshitz, L. Pitaevskii, and V. Berestetskii, _Vol 4 - Quantum Electrodynamics_. Mir editions, 1982.
  * (23) L. H. Ryder, _Quantum Field Theory_ , 2nd ed. Cambridge University Press, Jun. 1996. [Online]. Available: http://www.worldcat.org/isbn/0521478146
  * (24) G. Tatara, H. Kohno, and J. Shibata, “Microscopic approach to current-driven domain wall dynamics,” _Physics Reports_ , vol. 468, no. 6, p. 213–301, Nov 2008. [Online]. Available: http://dx.doi.org/10.1016/j.physrep.2008.07.003
  * (25) G. Yin, Y. Liu, Y. Barlas, J. Zang, and R. K. Lake, “Topological spin hall effect resulting from magnetic skyrmions,” _Phys. Rev. B_ , vol. 92, p. 024411, Jul 2015. [Online]. Available: https://link.aps.org/doi/10.1103/PhysRevB.92.024411
  * (26) S. John and A. Golubentsev, “Spin flux and magnetic solitons in an interacting two-dimensional electron gas: Topology of two-valued wave functions,” _Phys. Rev. B_ , vol. 51, pp. 381–398, Jan 1995. [Online]. Available: https://link.aps.org/doi/10.1103/PhysRevB.51.381
  * (27) J. C. Tully, _Perspective on ”Zur Quantentheorie der Molekeln”_. Berlin, Heidelberg: Springer Berlin Heidelberg, 2001, pp. 173–176. [Online]. Available: https://doi.org/10.1007/978-3-662-10421-7_3
  * (28) L. D. Landau and L. M. Lifshitz, _Quantum Mechanics, Non-Relativistic Theory, Third Edition: Volume 3_ , 3rd ed. Butterworth-Heinemann, Jan. 1981. [Online]. Available: http://www.worldcat.org/isbn/0750635398
  * (29) S. Komineas, C. Melcher, and S. Venakides, “The profile of chiral skyrmions of small radius,” _Nonlinearity_ , vol. 33, no. 7, p. 3395–3408, May 2020\. [Online]. Available: http://dx.doi.org/10.1088/1361-6544/ab81eb
  * (30) V. P. Voronov, B. A. Ivanov, and A. M. Kosevich, “Two-dimensional dynamic topological solitons in ferromagnets,” _Soviet Journal of Experimental and Theoretical Physics_ , vol. 57, no. 6, p. 1303, Jun. 1983. [Online]. Available: http://jetp.ras.ru/cgi-bin/e/index/e/57/6/p1303?a=list
  * (31) S. Komineas, C. Melcher, and S. Venakides, “Chiral skyrmions of large radius,” _Physica D: Nonlinear Phenomena_ , vol. 418, p. 132842, 2021. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0167278920308435
  * (32) G. M. Luke, L. P. Le, B. J. Sternlieb, Y. J. Uemura, J. H. Brewer, R. Kadono, R. F. Kiefl, S. R. Kreitzman, T. M. Riseman, C. E. Stronach, M. R. Davis, S. Uchida, H. Takagi, Y. Tokura, Y. Hidaka, T. Murakami, J. Gopalakrishnan, A. W. Sleight, M. A. Subramanian, E. A. Early, J. T. Markert, M. B. Maple, and C. L. Seaman, “Magnetic order and electronic phase diagrams of electron-doped copper oxide materials,” _Phys. Rev. B_ , vol. 42, pp. 7981–7988, Nov 1990. [Online]. Available: https://link.aps.org/doi/10.1103/PhysRevB.42.7981
  * (33) P. K. Mang, S. Larochelle, A. Mehta, O. P. Vajk, A. S. Erickson, L. Lu, W. J. L. Buyers, A. F. Marshall, K. Prokes, and M. Greven, “Phase decomposition and chemical inhomogeneity in nd2-xcexcuo4±δ,” _Physical Review B_ , vol. 70, p. 094507, 2004. [Online]. Available: https://doi.org/10.1103/PhysRevB.70.094507
  * (34) E. M. Motoyama, G. Yu, I. M. Vishik, O. P. Vajk, P. K. Mang, and M. Greven, “Spin correlations in the electron-doped high-transition temperature superconductor nd2-xcexcuo4±δ,” _Nature_ , vol. 445, p. 186, 2007. [Online]. Available: https://doi.org/10.1038/nature05437
  * (35) H. Saadaoui, Z. Salman, H. Luetkens, T. Prokscha, A. Suter, W. A. MacFarlane, Y. Jiang, K. Jin, R. L. Greene, E. Morenzoni, and R. F. Kiefl, “The phase diagram of electron-doped la2-xcexcuo4-δ,” _Nat Commun_ , vol. 6, p. 6041, 2015. [Online]. Available: https://doi.org/10.1038/ncomms7041
  * (36) K. Yamada, K. Kurahashi, T. Uefuji, M. Fujita, S. Park, S.-H. Lee, and Y. Endoh, “Commensurate spin dynamics in the superconducting state of an electron-doped cuprate superconductor,” _Physical Review Letters_ , vol. 90, p. 137004, Apr 2003. [Online]. Available: https://doi.org/10.1103/PhysRevLett.90.137004
  * (37) H. J. Kang, P. Dai, H. A. Mook, D. N. Argyriou, V. Sikolenko, J. W. Lynn, Y. Kurita, S. Komiya, and Y. Ando, “Electronically competing phases and their magnetic field dependence in electron-doped nonsuperconducting and superconducting pr0.88lace0.12cuo4±δ,” _Physical Review B_ , vol. 71, p. 214512, 2005. [Online]. Available: https://doi.org/10.1103/PhysRevB.71.214512
  * (38) W. Yu, J. S. Higgins, P. Bach, and R. L. Greene, “Transport evidence of a magnetic quantum phase transition in electron-doped high-temperature superconductors,” _Physical Review B_ , vol. 76, p. 020503(R), 2007. [Online]. Available: https://doi.org/10.1103/PhysRevB.76.020503
  * (39) A. Dorantes, A. Alshemi, Z. Huang, A. Erb, T. Helm, and M. V. Kartsovnik, “Magnetotransport evidence for irreversible spin reorientation in the collinear antiferromagnetic state of underdoped nd2-xcexcuo4,” _Physical Review B_ , vol. 97, p. 054430, 2018. [Online]. Available: https://doi.org/10.1103/PhysRevB.97.054430
  * (40) R. Ramazashvili, P. D. Grigoriev, T. Helm, F. Kollmannsberger, M. Kunz, W. Biberacher, E. Kampert, H. Fujiwara, A. Erb, J. Wosnitza, R. Gross, and M. V. Kartsovnik, “Experimental evidence for zeeman spin–orbit coupling in layered antiferromagnetic conductors,” _npj Quantum Mater._ , vol. 6, Apr 2021. [Online]. Available: https://doi.org/10.1038/s41535-021-00309-6
  * (41) N. P. Armitage, P. Fournier, and R. L. Greene, “Progress and perspectives on electron-doped cuprates,” _Reviews of Modern Physics_ , vol. 82, p. 2421, 2010. [Online]. Available: https://doi.org/10.1103/RevModPhys.82.2421
  * (42) N. P. Armitage, F. Ronning, D. H. Lu, C. Kim, A. Damascelli, K. M. Shen, D. L. Feng, H. Eisaki, Z.-X. Shen, P. K. Mang, N. Kaneko, M. Greven, Y. Onose, Y. Taguchi, , and Y. Tokura, “Doping dependence of an $n$-type cuprate superconductor investigated by angle-resolved photoemission spectroscopy,” _Physical Review Letters_ , vol. 88, p. 257001, Feb 2002. [Online]. Available: https://doi.org/10.1103/PhysRevLett.88.257001
  * (43) J. He, C. R. Rotundu, M. S. Scheurer, Y. He, M. Hashimoto, K.-J. Xu, Y. Wang, E. W. Huang, T. Jia, S. Chen, B. Moritz, D. Lu, Y. S. Lee, T. P. Devereaux, and Z.-X. Shen, “Fermi surface reconstruction in electron-doped cuprates without antiferromagnetic long-range order,” _PNAS_ , vol. 116, p. 3449, Feb 2019. [Online]. Available: https://doi.org/10.1073/pnas.181612111
  * (44) P. De Gennes, _Superconductivity of Metals and Alloys_. Perseus, Westview Press, 1966.
  * (45) M. Tinkham, _Introduction to Superconductivity_ , 2nd ed. Dover Publications, Jun. 2004. [Online]. Available: http://www.worldcat.org/isbn/0486435032
  * (46) V. P. Kravchuk, O. Gomonay, D. D. Sheka, D. R. Rodrigues, K. Everschor-Sitte, J. Sinova, J. van den Brink, and Y. Gaididei, “Spin eigenexcitations of an antiferromagnetic skyrmion,” _Phys. Rev. B_ , vol. 99, p. 184429, May 2019\. [Online]. Available: https://link.aps.org/doi/10.1103/PhysRevB.99.184429
  * (47) S. M. Rezende, A. Azevedo, and R. L. Rodríguez-Suárez, “Introduction to antiferromagnetic magnons,” _Journal of Applied Physics_ , vol. 126, no. 15, p. 151101, 2019. [Online]. Available: https://doi.org/10.1063/1.5109132
  * (48) C. Kittel, “Theory of antiferromagnetic resonance,” _Phys. Rev._ , vol. 82, pp. 565–565, May 1951. [Online]. Available: https://link.aps.org/doi/10.1103/PhysRev.82.565
  * (49) R. M. White, _Quantum Theory of Magnetism_ , 3rd ed. Springer-Verlag, 2007.