# Charged magnons on the surface of a topological Insulator

I. Martinez-Berumen<EMAIL_ADDRESS>W. A. Coish
<EMAIL_ADDRESS>T. Pereg-Barnea<EMAIL_ADDRESS>Department
of Physics, McGill University, Montréal, Québec, Canada H3A 2T8

###### Abstract

We study a system of two-dimensional Dirac electrons (as is realized on the
surface of a 3D topological insulator) coupled to an array of localized spins.
The spins are coupled ferromagnetically to each other, forming an ordered
ground state with low-energy spin-wave excitations (magnons). The Dirac
electrons couple to the spins through a spin-dependent effective Zeeman field.
The out-of-plane effective Zeeman field therefore serves as a Dirac mass that
gaps the electronic spectrum. Once a spin is flipped, it creates a surrounding
domain in which the sign of the Dirac mass is opposite to that of the rest of
the sample. Therefore, an electronic bound state appears on the domain wall,
as predicted by Jackiw and Rebbi. However, in a quantum magnet, a localized
spin flip does not produce an eigenstate. Instead, the eigenstates correspond
to delocalized spin waves (magnons). As in the case of the single flipped
spin, the delocalized magnon also binds an in-gap electronic state. We name
this excitation a ‘Jackiw-Rebbi-Magnon’ (JRM) and study its signature in the
dynamic spin susceptibility. When the sample is tunnel-coupled to an
electronic reservoir, a magnon produced in a system without any electrons
hybridizes with a JRM (which binds a single electron), producing magnon-JRM
polaritons. For such a system, we identify a quantum phase transition when the
magnon-JRM polariton energy falls below that of the fully polarized
ferromagnetic ground state.

††preprint: APS/123-QED

##  I Introduction

The spectrum of a strong 3D topological insulator (3DTI) is gapped in the bulk
and exhibits gapless surface states described by a Dirac cone, protected by
time-reversal symmetry [1, 2, 3]. The addition of magnetic moments on the
surface breaks time-reversal symmetry; these magnetic moments produce a local
effective Zeeman field for the Dirac electrons rendering them massive [4, 5,
6]. Such moments can be realized as magnetic impurities/dopants in topological
insulators, producing magnetic topological insulators [7, 8, 9, 10].
Equivalently, these magnetic moments can be introduced as a separate system in
proximity to the sample surface [11, 12, 13]. The surface magnetic moments may
interact with each other directly or through the Ruderman–Kittel–Kasuya–Yosida
(RKKY) interaction mediated by the Dirac electrons on the surface. It has been
shown that the associated RKKY interaction is likely to be ferromagnetic [4,
14]. When the moments are ordered ferromagnetically, the electrons acquire a
finite Dirac mass everywhere on the surface [4, 5], resulting in a spectral
gap. However, when the ferromagnet exhibits disorder (e.g., domains or
skyrmions), the Dirac mass may change sign, resulting in the closure of the
spectral gap [15, 16, 4].

In the presence of magnetic domain walls where the Dirac mass changes sign,
in-gap chiral bound states known as Jackiw-Rebbi modes appear [17]. These
states are exponentially localized to the domain wall as no extended states
can exist at an energy within the surface gap [1, 4]. While localized to
domain walls, the Jackiw-Rebbi states may carry currents along the domain
walls [18]. Domain-wall chiral currents have been seen to contribute to
experimental measurements such as the hysteretic magnetoconductance in
magnetic topological insulator surfaces [18, 19, 20, 21, 22]. Theoretical
models have also been proposed to study the formation and dynamics of domain
walls in an applied magnetic field [16].

In this paper, we consider a 3D topological insulator having surface electrons
that are coupled to an array of localized spins. This coupling may arise,
e.g., from exchange or from a local Zeeman interaction. Previous studies of 3D
topological insulators coupled to insulating magnets have described the
magnetic system classically [16, 23, 24]. In contrast, in this work we
consider a quantum ferromagnet having low-energy spin-wave (magnon)
excitations. These excitations can be considered as delocalized domain walls
or as a superposition of domains in different configurations. The consequence,
as we show, is the formation of delocalized chiral bound states that are
attached to the magnon and we therefore name them ‘Jackiw-Rebbi Magnons’ (see
Fig. 1).

The rest of this paper is organized as follows: In Sec. II we introduce the
model for the surface states, the magnetic impurities, and the coupling
between them. In Sec. III we present the form of the elementary excitations,
the Jackiw-Rebbi Magnons (JRMs) and, as an example, we numerically solve a
lattice model having an appropriate long-wavelength limit. In Sec. IV we
propose and calculate an experimental observable, the dynamic transverse spin
susceptibility, that can be used to probe JRM excitations. Finally, in Sec. V
we summarize our results and discuss the limitations and scope of the model.

Figure 1: (a) Magnetic impurities in the ferromagnetic ground state in a
magnetic topological insulator do not allow for electronic bound states. (b) A
local mass change produced by an impurity spin flip allows for in-gap Jackiw-
Rebbi modes. (c) A spin wave produces a delocalized domain wall with attached
delocalized bound states.

##  II Model

We consider a model of ferromagnetically coupled magnetic impurities
interacting with itinerant electrons on the surface of a 3D topological
insulator. The system Hamiltonian is

$H_{S}=H_{m}+H_{em},$ (1)

where $H_{m}$ is the Hamiltonian for the magnetic system and where $H_{em}$
describes the electrons in contact with the magnetic impurities.

The magnetic system is described by the (ferromagnetic) Heisenberg Hamiltonian
(setting $\hbar=1$):

$H_{m}=-J\sum_{\braket{ij}}\mathbf{S}_{i}\cdot\mathbf{S}_{j}-b\sum_{i}S^{z}_{i},$
(2)

where $J>0$ is the exchange interaction between nearest neighbors,
$\mathbf{S}_{i}$ is the spin operator for the impurity at site $i$, and $b$ is
the effective Zeeman splitting. We assume that the magnetic impurities form a
Bravais lattice, not necessarily commensurate with the underlying crystalline
lattice. In this equation, and in the rest of this paper, we set $\hbar=1$,
but $\hbar$ will be restored in several formulas below for clarity.

In the absence of the electronic system, the Heisenberg Hamiltonian $H_{m}$
has a fully spin-polarized ferromagnetic ground state. If, without loss of
generality, we assume $b>0$, this ground state is

$\ket{\mathrm{FM}}\equiv\ket{S,S,\ldots,S},$ (3)

where $\ket{m_{1},m_{2},\ldots,m_{N}}$ is a simultaneous eigenstate of all
operators $\\{S^{z}_{i}\\}$ ($i=1,2,\ldots,N$) with eigenvalues $\\{m_{i}\\}$
($S$ is the total impurity spin such that $|m_{i}|\leq S$). The energy of the
ferromagnetic ground state is

$E_{\mathrm{FM}}=-\frac{J}{2}NzS^{2}-bNS,$ (4)

where $z$ is the coordination number. The low-energy excitations are the spin
waves (magnons) $\ket{\mathbf{q}}$:
$H_{m}\ket{\mathbf{q}}=E_{\mathbf{q}}\ket{\mathbf{q}}$, with

$\ket{\mathbf{q}}=\frac{1}{\sqrt{2S}}S_{\mathbf{q}}^{-}\ket{\mathrm{FM}},\quad
S^{\pm}_{\mathbf{q}}\equiv\frac{1}{\sqrt{N}}\sum_{j}e^{i\mathbf{q}\cdot\mathbf{R}_{j}}S^{\pm}_{j}.$
(5)

The energy of the spin wave with wavevector $\mathbf{q}$ is

$\displaystyle E_{\mathbf{q}}$ $\displaystyle=$ $\displaystyle
E_{\mathrm{FM}}+\hbar\omega_{\mathbf{q}},$ (6)
$\displaystyle\hbar\omega_{\mathbf{q}}$ $\displaystyle=$ $\displaystyle
b+2JS\sum_{\boldsymbol{\delta}}\sin^{2}\left(\frac{\mathbf{q}\cdot\boldsymbol{\delta}}{2}\right).$
(7)

Here, we have introduced $\\{\boldsymbol{\delta}\\}$, the set of vectors
locating nearest neighbors and we have assumed a Bravais lattice so that there
is a nearest neighbor at $-\boldsymbol{\delta}$ for every
$\boldsymbol{\delta}$.

The surface electrons are modeled by massive Dirac fermions, with a mass that
depends on the state of the impurity spins:

$H_{em}=\int
d^{2}r\psi^{\dagger}(\mathbf{r})\left[-iv_{\mathrm{F}}\boldsymbol{\sigma}\cdot\boldsymbol{\nabla}+m(\mathbf{r})\sigma_{z}\right]\psi(\mathbf{r}).$
(8)

Here,
$\psi(\mathbf{r})=A_{\mathrm{surf}}^{-1/2}\sum_{\mathbf{k}}e^{i\mathbf{k}\cdot\mathbf{r}}c_{\mathbf{k}}$,
where $A_{\mathrm{surf}}$ is the two-dimensional surface area,
$c_{\mathbf{k}}=(c_{\mathbf{k}\uparrow},\,c_{\mathbf{k}\downarrow})^{T}$, and
where $c_{\mathbf{k}s}^{\dagger}$ creates an electron with wavevector
$\mathbf{k}=(k_{x},k_{y})^{T}$ and spin $s$. The mass $m(\mathbf{r})$ is an
_operator_ , which depends on the impurity spins:

$m(\mathbf{r})=m_{0}+\sum_{i}\mathcal{J}(\mathbf{r}-\mathbf{R}_{i})S_{i}^{z},$
(9)

where we have introduced a local finite-ranged exchange coupling
$\mathcal{J}(\mathbf{r}-\mathbf{R}_{i})$ between each magnetic impurity $i$ at
position $\mathbf{R}_{i}$ and the electron spin density at position
$\mathbf{r}$. The constant background mass $m_{0}$ breaks time-reversal
symmetry; it may arise either from a global magnetic field or from some other
background magnetization.

The system Hamiltonian [Eq. (1)] includes terms that may arise from an out-of-
plane magnetic field ($\propto m_{0},\propto b$) through a Zeeman coupling to
the electrons, but we neglect the diamagnetic shifts and other orbital effects
that would arise from a magnetic vector potential. This is justified for small
impurity clusters of size $A_{\mathrm{surf}}$ and in a sufficiently weak
magnetic field $B$: $A_{\mathrm{surf}}B\ll\Phi_{0}$, where $\Phi_{0}=h/e$ is
the flux quantum. In addition, we neglect terms like
$\sim\mathcal{J}_{\parallel}S^{\nu}_{j}\psi^{\dagger}(\mathbf{r})\sigma_{\nu}\psi(\mathbf{r})$
[$\nu\in\left\\{x,y\right\\}$] that couple the in-plane components of the
impurity spins to the electron spins. These terms may be ignored in the
presence of a large exchange anisotropy,
$|\mathcal{J}(\mathbf{r})|\gg|\mathcal{J}_{\parallel}|$, or for a large
effective Zeeman term $|b|\gg\sqrt{2S}|\mathcal{J}_{\parallel}|$.

In the next section, we show that the elementary excitations of $H_{S}$ [Eq.
(1)] have the same form as the magnon states given by Eq. (5). The key
difference, as shown below, is that these new excitations come with an
electronic state localized to the vicinity of the lowered spin (see Fig. 1).
To create these excitations, starting from the ferromagnetic ground state, in
addition to a transverse magnetic field (which creates the spin excitation),
we couple our system to an electronic reservoir, which provides the electron
that occupies the localized electronic state. This is shown in Sec. IV (see
Fig. 3).

##  III Jackiw-Rebbi magnons

The total impurity spin in the $z$-direction,
$S_{\mathrm{tot}}^{z}=\sum_{i}S_{i}^{z}$, commutes with $H_{S}$. The
Hamiltonian $H_{S}$ for the combined electronic and magnetic system can
therefore be diagonalized independently in sectors where
$S^{z}_{\mathrm{tot}}$ is fixed. For the ferromagnetic ground state we have
$S^{z}_{\mathrm{tot}}=NS$, giving

$\braket{\mathrm{FM}}{H_{S}}{\mathrm{FM}}=E_{\mathrm{FM}}+H_{\mathrm{FM}},$
(10)

where the effective Hamiltonian
$H_{\mathrm{FM}}=\braket{\mathrm{FM}}{H_{em}}{\mathrm{FM}}$ is obtained from
Eq. (8) with the replacement
$m(\mathbf{r})\to\braket{\mathrm{FM}}{m(\mathbf{r})}{\mathrm{FM}}=m_{\mathrm{FM}}(\mathbf{r})$.
The mass profile $m_{\mathrm{FM}}(\mathbf{r})$ is

$m_{\mathrm{FM}}(\mathbf{r})=m_{0}+S\sum_{i}\mathcal{J}(\mathbf{r}-\mathbf{R}_{i}),$
(11)

and $E_{\mathrm{FM}}$ is the ferromagnetic ground-state energy, as given in
Eq. (4).

Provided that $\mathcal{J}(\mathbf{r})>0$ $\forall\mathbf{r}$ (and if
$m_{0}>0$), the mass $m_{\mathrm{FM}}(\mathbf{r})$ is finite and positive
everywhere, leading to a gapped single-particle energy spectrum for
$H_{\mathrm{FM}}$. The eigenstates of $H_{S}$, with $S_{\mathrm{tot}}^{z}=NS$,
are then given by a tensor product of electronic and magnetic states; the
electronic states are eigenstates of $H_{\mathrm{FM}}$ and the magnetic system
is described by the ferromagnetic ground state $\ket{\mathrm{FM}}$.

In the subspace of one spin lowered, $S^{z}_{\mathrm{tot}}=SN-1$, the magnetic
state can be constructed from the basis $\\{\ket{\mathbf{R}_{i}}\\}$,
$i=1,\ldots,N$, where

$\ket{\mathbf{R}_{i}}\equiv\frac{1}{\sqrt{2S}}S_{i}^{-}\ket{\mathrm{FM}}.$
(12)

Similar to $H_{\mathrm{FM}}$, we define $H_{i}$, the effective Hamiltonian for
the electrons when the magnetic system is frozen in the state
$\ket{\mathbf{R}_{i}}$:

$H_{i}=\braket{\mathbf{R}_{i}}{H_{em}}{\mathbf{R}_{i}},$ (13)

which is given by Eq. (8) with the replacement $m(\mathbf{r})\to
m_{i}(\mathbf{r})=\braket{\mathbf{R}_{i}}{m(\mathbf{r})}{\mathbf{R}_{i}}$:

$m_{i}(\mathbf{r})=m_{\mathrm{FM}}(\mathbf{r})-\mathcal{J}(\mathbf{r}-\mathbf{R}_{i}).$
(14)

This is the mass term produced by the impurities when the magnetic system is
in the state $\ket{\mathbf{R}_{i}}$ [see Fig. 2(a)].

###  III.1 Single-particle Jackiw-Rebbi magnons

The Hamiltonian $H_{i}$ has single-particle eigenstates $\ket{\alpha}_{i}$:

$H_{i}\ket{\alpha}_{i}=\varepsilon_{\alpha}\ket{\alpha}_{i},$ (15)

where the energies $\varepsilon_{\alpha}$ are independent of $i$ due to the
discrete translational invariance of the original Hamiltonian.

We now introduce the composite state made of the $\alpha^{\mathrm{th}}$
electronic eigenstate of $H_{i}$ and the spin state $\ket{\mathbf{R}_{i}}$,

$\ket{\alpha,\mathbf{R}_{i}}\equiv\ket{\alpha}_{i}\otimes\ket{\mathbf{R}_{i}}.$
(16)

The matrix elements of $H_{S}$ in this basis are:

$\begin{split}\braket{\alpha,\mathbf{R}_{i}}{H_{S}}{\beta,\mathbf{R}_{j}}=&\left(E_{\mathrm{FM}}+b+JSz+\varepsilon_{\alpha}\right)\delta_{ij}\delta_{\alpha\beta}\\\
&-\gamma_{\alpha\beta}JS\sum_{\boldsymbol{\delta}}\delta_{\mathbf{R}_{i},\mathbf{R}_{j}+\boldsymbol{\delta}}.\end{split}$
(17)

Above, we have introduced

$\gamma_{\alpha\beta}\equiv\prescript{}{i}{\braket{\alpha}{\beta}}_{j},\quad
i,j\in\braket{ij},$ (18)

the electronic overlap between the states $\ket{\alpha}_{i}$ and
$\ket{\beta}_{j}$. The notation $i,j\in\langle ij\rangle$ indicates that the
definition holds for all nearest neighbors $i,j$. For simplicity, we assume a
high-symmetry impurity lattice, so that all nearest-neighbor bonds are
equivalent.

Figure 2: (a) Mass profile, $m_{i}(\mathbf{r})$, given by Eq. (14) when
$\mathbf{R}_{i}=(35,35)$. Here, the exchange $\mathcal{J}(\mathbf{r})$ is
given by Eq. (31). The black line denotes the region where
$m_{i}(\mathbf{r})=0$, resulting in a sign change in the mass leading to a
domain-wall bound state. (b) Probability density (color scale) and local
velocity (red arrows) for the lowest positive-energy eigenstate,
$\ket{0}_{i}$, of $H^{\mathrm{latt}}_{i}$ [Eq. (30)]. (c) The spectrum of the
Hamiltonian $H_{i}^{\mathrm{latt}}$ is symmetric about $E=0$ as a consequence
of particle-hole symmetry. Isolated energies are found in the gap,
corresponding to bound states, localized to the domain wall. The values for
the parameters used are given in Sec. III.3. (d) Velocity profile for either
$\ket{\psi,\mathbf{q}}$ or $\ket{\Psi,\mathbf{q}}$, given by Eq. (34). The
black dots give the positions of magnetic impurities.

Jackiw-Rebbi bound states arise if the mass term $m_{i}(\mathbf{r})$ [Eq.
(14)] change sign in the vicinity of the flipped spin. In the rest of this
paper, we focus our attention on the lowest positive-energy single-particle
bound state. Although the analysis below also applies to any of the bound
states, we concentrate on the lowest positive-energy bound state. This choice
is motivated by the many-body case where we begin at half filling and the
lowest-energy unoccupied state is then populated. We label this bound state by
$\alpha=0$, such that $H_{i}\ket{0}_{i}=\varepsilon_{0}\ket{0}_{i}$ [see Fig.
2(c)]. The localized Jackiw-Rebbi state $\ket{0,\mathbf{R}_{i}}$ is coupled to
the states $\ket{\alpha,\mathbf{R}_{j}}$ with coupling strength
$JS\gamma_{0\alpha}$. If the separation in energy between this state and other
electronic states is large relative to the coupling, i.e.

$JS\gamma_{0\alpha}\ll|\varepsilon_{0}-\varepsilon_{\alpha}|,\quad(\alpha\neq
0),$ (19)

then we can apply degenerate perturbation theory and diagonalize the
Hamiltonian in the subspace spanned by the set of degenerate states
$\\{\ket{0,\mathbf{R}_{1}},\ldots,\ket{0,\mathbf{R}_{N}}\\}$. In this
subspace, the effective Hamiltonian is equivalent to a tight-binding model
with nearest-neighbor hopping proportional to the overlap
$\gamma\equiv\gamma_{00}$ [see Eqs. (17-18)]. The eigenstates of this model
are plane waves of the form

$\ket{\psi,\mathbf{q}}=\frac{1}{\sqrt{N}}\sum_{j}e^{i\mathbf{q}\cdot\mathbf{R}_{j}}\ket{0,\mathbf{R}_{j}},$
(20)

which we call Jackiw Rebbi Magnons (JRMs). These states differ from
conventional spin waves; in addition to carrying a single quantum of angular
momentum (spin), they carry a single charge associated with a bound electron.
Because these states are derived from a tight-binding model, they inherit a
simple tight-binding dispersion:

$\tilde{E}_{\mathbf{q}}=E_{\mathrm{FM}}+\hbar\tilde{\omega}_{\mathbf{q}},$
(21)

with

$\hbar\tilde{\omega}_{\mathbf{q}}=\varepsilon_{0}+b+JSz(1-\gamma)+2JS\gamma\sum_{\boldsymbol{\delta}}\sin^{2}\left(\frac{\mathbf{q}\cdot\boldsymbol{\delta}}{2}\right).$
(22)

###  III.2 Many-body Jackiw-Rebbi magnons

The analysis above is not restricted to single-particle states. For example,
the half-filled ($N_{e}$-particle) state resulting from filling the negative-
energy single-particle eigenstates of $H_{i}$ is

$\ket{\mathrm{hf},\mathbf{R}_{i}}\equiv\ket{\mathrm{hf}}_{i}\otimes\ket{\mathbf{R}_{i}};\quad\ket{\mathrm{hf}}_{i}=\prod_{\alpha\,(\varepsilon_{\alpha}<0)}c_{\alpha
i}^{\dagger}\ket{\mathrm{vac}}_{e},$ (23)

where $c_{\alpha i}^{\dagger}\ket{\mathrm{vac}}_{e}=\ket{\alpha}_{i}$ with
$\ket{\mathrm{vac}}_{e}$ the electronic vacuum. The minimum-energy
$N_{e}$-conserving excitation out of the half-filled state has a finite energy
cost $\Delta\varepsilon=2\varepsilon_{0}$ (the spectrum of $H_{i}$ is
symmetric around zero because it is particle-hole symmetric [25]). A general
$N_{e}$-particle eigenstate can be written as

$\begin{split}\ket{\boldsymbol{\alpha}_{N_{e}},\mathbf{R}_{j}}&\equiv\ket{\boldsymbol{\alpha}_{N_{e}}}_{j}\otimes\ket{\mathbf{R}_{j}},\\\
\ket{\boldsymbol{\alpha}_{N_{e}}}_{j}&=\left(c^{\dagger}_{\alpha_{1}j}\cdots
c^{\dagger}_{\alpha_{N_{e}}j}\right)\ket{\mathrm{vac}}_{e}.\end{split}$ (24)

If the minimum excitation energy, $2\varepsilon_{0}$, is large compared to the
matrix element that connects the states $\ket{\mathrm{hf},\mathbf{R}_{i}}$
with $\ket{\boldsymbol{\alpha}_{N_{e}},\mathbf{R}_{j}}$, that is

$JS\prescript{}{i}{\braket{\mathrm{hf}}{\boldsymbol{\alpha}_{N_{e}}}}_{j}\ll
2\varepsilon_{0},\quad i,j\in\braket{ij},$ (25)

then we can diagonalize the Hamiltonian in the degenerate subspace consisting
of
$\\{\ket{\mathrm{hf},\mathbf{R}_{1}},\ldots,\ket{\mathrm{hf},\mathbf{R}_{N}}\\}$
to obtain the states

$\ket{\mathrm{hf},\mathbf{q}}\equiv\frac{1}{\sqrt{N}}\sum_{j}e^{i\mathbf{q}\cdot\mathbf{R}_{j}}\ket{\mathrm{hf},\mathbf{R}_{j}}.$
(26)

Using the same argument, $H_{S}$ has approximate $(N_{e}+1)$-particle
eigenstates (many-body JRMs) of the form

$\ket{\Psi,\mathbf{q}}\equiv\frac{1}{\sqrt{N}}\sum_{j}e^{i\mathbf{q}\cdot\mathbf{R}_{j}}c_{0j}^{\dagger}\ket{\mathrm{hf},\mathbf{R}_{j}}.$
(27)

We use the term _many-body_ Jackiw-Rebbi magnons specifically for these
$(N_{e}+1)$-particle states which, in addition to carrying a quantum of
angular momentum (spin), also carry an additional electronic charge beyond
half filling. To differentiate this many-body state from the single-particle
state $\ket{\psi,\mathbf{q}}$ we use the capital letter $\Psi$. Note that the
hopping between $(N_{e}+1)$-particle states with a single spin-flip,
$c_{0j}^{\dagger}\ket{\mathrm{hf},\mathbf{R}_{j}}$, is now proportional to the
many-body overlap

$\gamma\to\gamma_{\Psi}=\braket{\mathbf{R}_{i},\mathrm{hf}}{c_{0i}c^{\dagger}_{0j}}{\mathrm{hf},\mathbf{R}_{j}},(i,j)\in\langle
ij\rangle,$ (28)

which may decrease rapidly as the number of particles increases.

###  III.3 Numerical Example

As a concrete example, we consider a square-lattice model that recovers
$H_{i}$ in the long-wavelength limit. This allows us to find lattice versions
of the eigenstates $\ket{0}_{i}$. The lattice model is obtained by first
rewriting the kinetic-energy term in Eq. (8) in $\mathbf{k}$-space:

$v_{\mathrm{F}}\int
d^{2}r\psi^{\dagger}(\mathbf{r})\left[-i\boldsymbol{\sigma}\cdot\boldsymbol{\nabla}\right]\psi(\mathbf{r})=v_{\mathrm{F}}\sum_{\mathbf{k}}c_{\mathbf{k}}^{\dagger}\mathbf{k}\cdot\boldsymbol{\sigma}c_{\mathbf{k}},$
(29)

then substituting $k_{\nu}\rightarrow a_{0}^{-1}\sin
k_{\nu}a_{0},\nu\in\left\\{x,y\right\\}$. To maintain a single Dirac cone (in
the limit $m_{0}\to 0,\,\mathcal{J}(\mathbf{r})\to 0$) at the $\Gamma$-point
and to avoid extra Dirac points at the Brillouin-zone corners, we furthermore
introduce an additional term $-2A\sigma_{z}(2-\cos k_{x}a_{0}-\cos
k_{y}a_{0})$. The equivalent term in a continuum model is
$-Aa_{0}^{2}\sigma_{z}k^{2}$ [26]. The resulting real-space lattice
Hamiltonian is therefore

$\begin{split}H_{i}^{\mathrm{latt}}=&\sum_{\mathbf{r}}c^{\dagger}_{\mathbf{r}}m_{i}(\mathbf{r})\sigma_{z}c_{\mathbf{r}}\\\
&+\frac{i}{2}\frac{\hbar
v_{\mathrm{F}}}{a_{0}}\sum_{\mathbf{r},\nu}\left(c^{\dagger}_{\mathbf{r}+a_{0}\mathbf{e}_{\nu}}\sigma_{\nu}c_{\mathbf{r}}-\mathrm{h.c.}\right)\\\
&+A\sum_{\mathbf{r},\nu}\left(c^{\dagger}_{\mathbf{r}+a_{0}\mathbf{e}_{\nu}}\sigma_{z}c_{\mathbf{r}}+\mathrm{h.c.}-2c_{\mathbf{r}}^{\dagger}\sigma_{z}c_{\mathbf{r}}\right),\end{split}$
(30)

where $\nu\in\left\\{x,y\right\\}$.

For the simulations (Fig. 2) we consider a $5\times 5$ impurity grid with
impurity spacing $a=14a_{0}$ on top of a crystalline lattice of size
$N_{x}=N_{y}=70$, where $N_{\nu}$ is the number of lattice sites in the
$\nu$-direction and consider periodic boundary conditions. Although here we
consider a square lattice of impurity spins for this particular numerical
example, note that the results obtained in Sec. III apply for any choice of
Bravais lattice. We use $A/(\hbar v_{\mathrm{F}}/a_{0})=1,m_{0}/(\hbar
v_{\mathrm{F}}/a_{0})=-1/2,S=1/2$. For
$\mathcal{J}(\mathbf{r}-\mathbf{R}_{i})$, we use the exponentially decaying
function

$\mathcal{J}(\mathbf{r}-\mathbf{R}_{i})=\mathcal{J}_{0}e^{-|\mathbf{r}-\mathbf{R}_{i}|/\lambda},$
(31)

with $\mathcal{J}_{0}/(\hbar v_{\mathrm{F}}/a_{0})=2,\lambda=a/2$, and set
$v_{\mathrm{F}}=a_{0}=1$. The parameters chosen here lead to a low-energy
bound state encircling the flipped impurity spin associated with a circulating
(chiral) current. The specific choice of parameters taken above is not
necessary to obtain this bound state. The only requirement is that there is a
gap closure in the vicinity of the flipped impurity spin. In the continuum
model, this requirement amounts to a change in sign of $m_{i}(\mathbf{r})$,
while in the lattice model this requires a change in sign of
$m_{i}(\mathbf{r})[m_{i}(\mathbf{r})-8A]$ in the vicinity of the lowered
impurity (see Eq. (3.91) in Ref. [26]). Even when the eigenstates vary on a
scale comparable to the lattice constant (far from the long-wavelength limit),
the qualitative physics of electronic states bound to magnetic domain walls
will be preserved.

Figure 2(a) shows the spatially dependent Dirac mass $m_{i}(\bm{r})$ [Eq.
(14)] resulting from a flipped impurity spin. Figure 2(b) shows the lowest-
energy Jackiw-Rebbi-like bound state in the vicinity of the flipped impurity
spin. In Fig 2(c) the electronic density of states of $H_{i}^{\mathrm{latt}}$
is shown, where bound states are clearly visible within the gap. The number
and energies of these in-gap bound states vary with the size of the negative-
mass domain [15]. The local velocity for $\ket{0}_{i}$ and
$\ket{\psi,\mathbf{q}}$ are shown in Figs. 2(b) and 2(d), respectively. The
local velocity is evaluated from the Heisenberg equation of motion,

$\mathbf{v}(t)\equiv\dot{\hat{\mathbf{r}}}(t)=\frac{i}{\hbar}[H_{i}^{\mathrm{latt}},\hat{\mathbf{r}}(t)],$
(32)

where
$\hat{\mathbf{r}}=\sum_{\mathbf{r}}\mathbf{r}c^{\dagger}_{\mathbf{r}}\sigma_{0}c_{\mathbf{r}}$
is the position operator. Therefore, the local velocity operator in the
$\nu$-direction is

$\begin{split}v_{\nu}(\mathbf{r})=\left(\frac{v_{\mathrm{F}}}{2}c^{\dagger}_{\mathbf{r}+\mathbf{e}_{\nu}}\sigma_{\nu}c_{\mathbf{r}}-iA\frac{a_{0}}{\hbar}c^{\dagger}_{\mathbf{r}+\mathbf{e}_{\nu}}\sigma_{z}c_{\mathbf{r}}+\mathrm{h.c.}\right).\end{split}$
(33)

The JRM [Eq. (20)] and the many-body JRM [Eq. (27)] have the same velocity
profile since

$\begin{split}\braket{\Psi,\mathbf{q}}{\mathbf{v}(\mathbf{r})}{\Psi,\mathbf{q}}&=\frac{1}{N}\sum_{i}\braket{\mathrm{hf},\mathbf{R}_{i}}{c_{0i}\mathbf{v}(\mathbf{r})c^{\dagger}_{0i}}{\mathrm{hf},\mathbf{R}_{i}}\\\
&=\frac{1}{N}\sum_{i}\prescript{}{i}{\braket{0}{\mathbf{v}(\mathbf{r})}{0}}_{i}.\end{split}$
(34)

The off-diagonal elements in the expression above,
$\braket{\mathrm{hf},\mathbf{R}_{i}}{\mathbf{v}(\mathbf{r})}{\mathrm{hf},\mathbf{R}_{j\neq
i}}$, vanish because $v_{\nu}(\mathbf{r})$ only acts on the electronic space.
The average velocity is then independent of the momentum $\mathbf{q}$ of the
Jackiw-Rebbi magnon. Moreover, the contribution from the negative-energy
states vanishes because the Hamiltonian $H_{i}$ has particle-hole symmetry and
the velocity operator is odd under the particle-hole transformation.

Including the orbital effect of a vector potential can break the particle-hole
symmetry and thus alter the calculation above. However, a weak out-of-plane
magnetic field will preserve the existence of in-gap domain-wall bound states
[15], leading to a similar qualitative picture to that presented here.

##  IV Spin Susceptibility

The dynamic transverse spin susceptibility can be used to probe the existence
of JRMs. In contrast with the case of a Heisenberg Hamiltonian [Eq. (2)], for
which a time-varying transverse magnetic field directly excites purely
magnetic excitations (magnons), the transverse magnetic field alone is not
sufficient to excite JRMs, which also carry a bound electronic charge. For
this reason, in addition to a time-varying transverse magnetic field, here we
also incorporate an electronic reservoir (see Fig. 3).

The Hamiltonian of our system (the topological-insulator surface electrons and
magnetic impurities) subjected to a transverse magnetic field and coupled to
an electronic reservoir is given by

$\displaystyle H(t)$ $\displaystyle=$ $\displaystyle H_{0}+V(t),$ (35)
$\displaystyle H_{0}$ $\displaystyle=$ $\displaystyle H_{S}+H_{R}+H_{T},$ (36)

where
$V(t)\propto\sum_{j}\mathbf{B}_{\perp}(\mathbf{R}_{j},t)\cdot\mathbf{S}_{j}$
contains the spatially and temporally varying transverse magnetic field
$\mathbf{B}_{\perp}=\left(B_{x},B_{y},0\right)^{T}$ (we neglect the electron-
field coupling as we consider the case in which the impurity-field coupling
dominates), $H_{S}$ is defined in Eq. (1), $H_{R}$ describes the electronic
reservoir, and $H_{T}$ describes tunneling to the reservoir:

$\displaystyle H_{R}$ $\displaystyle=$
$\displaystyle\sum_{\eta}\epsilon_{\eta}a_{\eta}^{\dagger}a_{\eta},$ (37)
$\displaystyle H_{T}$ $\displaystyle=$
$\displaystyle\sum_{\eta,j}\left(t_{\eta
j}c_{0j}^{\dagger}a_{\eta}+\mathrm{h.c.}\right)\ket{\mathbf{R}_{j}}\bra{\mathbf{R}_{j}}.$
(38)

The operator $a^{\dagger}_{\eta}$ creates an electron with energy
$\epsilon_{\eta}$ in the reservoir, where $\eta$ is a general quantum number.
We have included only those terms in $H_{T}$ that couple to the relevant mid-
gap bound states [Fig. 3(a)] under the assumption that all other levels are
far from resonance. In general, the tunnel coupling $t_{\eta j}$ associated
with $c_{0j}^{\dagger}$ will depend on the impurity $j$, but for simplicity
here we take $t_{\eta j}=t_{\eta}$ (independent of $j$), which can be realized
when the extended reservoir states couple equivalently to all impurity sites.
This choice guarantees that $H_{T}$ preserves the wavevector $\mathbf{q}$ in
converting a spin wave to a JRM:
${}_{e}\bra{\mathrm{vac}}\bra{\mathbf{q}}H_{T}\ket{\psi,\mathbf{q}^{\prime}}\propto\delta_{\mathbf{q},\mathbf{q}^{\prime}}$.

Figure 3: The electronic states available when the magnetic system is in the
ferromagnetic ground state $\ket{\mathrm{FM}}$ are not resonant with the
reservoir. The energy range of the reservoir is such that the spin-wave state,
with excitation energy $\omega_{\mathbf{q}}$ [Eq. (7)], may only be resonant
with the JRMs, with energy $\tilde{\omega}_{\mathbf{q}}-\epsilon_{\eta}$ [Eq
(22)] (here, $\epsilon_{\eta}$ is the energy of the reservoir electron that is
removed to populate the Jackiw-Rebbi state). The JRM excitations are only
accessible after a spin is flipped (magnon is created) through a transverse
magnetic field $\mathbf{B}_{\perp}(t)$. The tunneling Hamiltonian $H_{T}$
couples the electronic reservoir with the surface states of the 3D topological
insulator [see Eq. (38)].

The zero-temperature dynamic transverse susceptibility is

$\chi^{+-}(\mathbf{q},\omega)=\int dte^{i\omega t-0^{+}|t|}C(\mathbf{q},t).$
(39)

Here, the spin-spin correlation function is given by

$C(\mathbf{q},t)=-i\theta(t)\braket{\mathrm{gs}}{[\hat{S}_{\mathbf{q}}^{+}(t),S^{-}_{-\mathbf{q}}]}{\mathrm{gs}},$
(40)

where $\hat{S}^{+}_{\mathbf{q}}(t)=e^{iH_{0}t}S_{\mathbf{q}}^{+}e^{-iH_{0}t}$
is the time-dependent raising operator in the interaction picture and
$\ket{\mathrm{gs}}$ is the intially prepared ground state of $H_{0}$.

We now consider a special initial condition that allows for an exact solution.
In particular, we consider the fully spin-polarized ferromagnetic ground
state, a filled band of reservoir states, but no electrons occupying any of
the surface states of the topological insulator. Such an initial condition is
an eigenstate of the effective Hamiltonian $H_{0}$, but in practice it would
only be metastable, eventually collapsing due to non-resonant tunneling terms
that have been neglected in $H_{T}$. We discuss the more realistic case of a
filled surface-state valence band, below. Explicitly, the assumed ground state
is:

$\ket{\mathrm{gs}}=\ket{\mathrm{vac}}_{e}\otimes\ket{\mathrm{FM}}\otimes\ket{\mathrm{FB}}_{R},$
(41)

where $\ket{\mathrm{FB}}_{R}$ represents the filled reservoir band. For an
initially filled reservoir, the tunneling Hamiltonian only leads to
transitions within a restricted subspace of states, spanned by

$\displaystyle\ket{\psi(0)}$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{2S}}S_{-\mathbf{q}}^{-}\ket{\mathrm{gs}},$ (42)
$\displaystyle\ket{\eta}$ $\displaystyle=$
$\displaystyle\ket{\psi,\mathbf{-q}}\otimes a_{\eta}\ket{\mathrm{FB}}_{R}.$
(43)

Since $H_{T}$ preserves the wavevector $\mathbf{q}$, we have suppressed the
explicit $\mathbf{q}$-dependence in
$\left\\{\ket{\psi(0)},\,\ket{\eta}\right\\}$, above.

We can directly solve the time-dependent Schrödinger equation for this
problem,

$i\partial_{t}\ket{\psi(t)}=H_{0}\ket{\psi(t)},$ (44)

giving a solution (see Refs. 27, 28 for equivalent solutions):

$\ket{\psi{(t)}}=\psi_{0}(t)\ket{\psi(0)}+\sum_{\eta}\psi_{\eta}(t)\ket{\eta}.$
(45)

This solution then gives the correlation function $C(\mathbf{q},t)$ (and
hence, the susceptibility):

$C(\mathbf{q},t)=-i\theta(t)e^{iE_{\mathrm{gs}}t}2S\braket{\psi(0)}{\psi(t)}.$
(46)

Here, the ground state energy (defined by
$H_{0}\ket{\mathrm{gs}}=E_{\mathrm{gs}}\ket{\mathrm{gs}}$) is

$E_{\mathrm{gs}}\equiv E_{\mathrm{FM}}+\sum_{\eta}\epsilon_{\eta}.$ (47)

The exact susceptibility is then

$\chi^{+-}(\mathbf{q},\omega)=\frac{2S}{\omega-\omega_{\mathbf{q}}-\Sigma_{\mathbf{q}}(\omega)},$
(48)

where the self-energy is

$\Sigma_{\mathbf{q}}(\omega)=\sum_{\nu}\frac{|t_{\eta}^{2}|}{\omega-\left(\tilde{\omega}_{\mathbf{q}}-\epsilon_{\eta}\right)+i0^{+}}.$
(49)

In the equations above, $\omega_{\mathbf{q}}$ and
$\tilde{\omega}_{\mathbf{q}}-\epsilon_{\eta}$ correspond to the excitation
energy for a spin wave and a JRM with momentum $\mathbf{q}$, respectively [cf.
Eqs. (7) and (22)]. We additionally introduce phenomenological decay rates for
the JRMs ($\Gamma_{\mathrm{JRM}}$) and for the magnons ($\Gamma_{m}$) with the
replacements

$\displaystyle 0^{+}$ $\displaystyle\to$ $\displaystyle\Gamma_{\mathrm{JRM}},$
(50) $\displaystyle\Sigma_{\mathbf{q}}(\omega)$ $\displaystyle\to$
$\displaystyle\Sigma_{\mathbf{q}}(\omega)-i\Gamma_{m}.$ (51)

Evidence of JRMs can be found by analyzing peaks of the spectral function,

$A(\mathbf{q},\omega)=-2\,\mathrm{Im}\,\chi^{+-}(\mathbf{q},\omega).$ (52)

Due to the tunneling term $H_{T}$, neither the spin-wave $\ket{\psi(0)}\propto
S_{-\mathbf{q}}\ket{\mathrm{gs}}$, nor the JRM with a reservoir hole at energy
$\epsilon_{\eta}$, $\ket{\eta}$, is an eigenstate of $H_{0}$. As a result,
$A(\mathbf{q},\omega)$ shows multiple peaks associated with the excitation
energies for hybridized magnon-JRM polaritons (see the example illustrated in
Fig. 4). The characteristic double-peaked structure and avoided crossing is
distinct from the case of an isolated Heisenberg ferromagnet, having
susceptibility

$\chi_{0}^{+-}(\mathbf{q},\omega)=\frac{2S}{\omega-\omega_{\mathbf{q}}+i0^{+}}.$
(53)

Here, the corresponding spectral function
$A_{0}(\mathbf{q},\omega)=-2\,\mathrm{Im}\,\chi_{0}^{+-}(\mathbf{q},\omega)$
is given by a Dirac delta distribution centered at the spin-wave excitation
energy $\omega_{\mathbf{q}}$:
$A_{0}(\mathbf{q},\omega)\propto\delta(\omega-\omega_{\mathbf{q}})$.

Figure 4: Spectral function for a narrow-band reservoir ($W\to 0$) as a
function of the reservoir energy $\epsilon_{0}$. The two peaks correspond to
magnon-JRM polaritons. A quantum phase transition (QPT) occurs when the
excitation energy for the lowest-energy magnon-JRM polariton becomes negative.
In this plot we have chosen the parameters
$\Gamma=1,S=1/2,b=1/2,J=1/2,\varepsilon_{0}=5,\gamma=1/10,\Gamma_{m}=\Gamma_{\mathrm{JRM}}=1/10$,
where all energies are given in units of $\Gamma$.

As an illustrative example, we consider the reservoir band to form a quasi-
continuum of energies
$\epsilon_{\eta}\in\left[\epsilon_{0},\epsilon_{0}+W\right]$, as shown in Fig.
3(a). We additionally assume a constant density of states $\rho$ and constant
tunnel coupling $t_{\eta}=t$. In the narrow-band limit, $W\rightarrow 0$, the
self-energy is then

$\Sigma_{\mathbf{q}}(\mathbf{q},\omega)=\frac{\Gamma/\pi}{\omega-(\tilde{\omega}_{\mathbf{q}}-\epsilon_{0})+i0^{+}},$
(54)

where $\Gamma=\pi\rho|t|^{2}$. In this simple limit, the susceptibility has
only two poles corresponding to two magnon-JRM polaritons. The magnon-like
resonance with frequency $\sim\omega_{\mathbf{q}}$ couples to the JRM-like
resonance having frequency $\sim\tilde{\omega}_{\mathbf{q}}-\epsilon_{0}$,
with an avoided crossing of width $\sim\Gamma$ (Fig. 4).

A quantum phase transition (QPT) takes place if the lowest-energy magnon-JRM
polariton falls in energy below the ferromagnetic ground state. In the case
illustrated in Fig. 4, this corresponds to

$\omega_{\mathbf{q}}(\tilde{\omega}_{\mathbf{q}}-\epsilon_{0})<\left(\Gamma/\pi\right)^{2}.$
(55)

Above a critical value of the reservoir energy
$\epsilon_{0}>\epsilon_{\mathrm{crit}}=\tilde{\omega}_{\mathbf{q}}-\Gamma^{2}/(\omega_{\mathbf{q}}\pi^{2})$,
the fully spin-polarized ground state becomes unstable to the formation of
JRMs. The observation of such a transition due to a coupled reservoir could be
a robust indicator of the formation of JRMs. This transition can also be seen
in the spectral function (Fig. 4), where the lower magnon-JRM polariton peak
moves to negative frequency at large $\epsilon_{0}$.

In Eq. (41) we considered an initial ground state with empty bands for the
topological-insulator surface, but a filled band of reservoir states. Such an
initial condition is only metastable if the reservoir electrons have a finite
amplitude for tunneling to the surface. A more realistic (but computationally
more complex) situation is therefore realized when the topological insulator
surface is at half-filling, corresponding to a global chemical potential
positioned just above the top of the reservoir band
($\mu=\epsilon_{0}+W+0^{+}$ in Fig. 3). In this case, the ground state is
given by

$\ket{\mathrm{hf},\mathrm{FM}}\equiv\left(\prod_{v\,(\varepsilon_{v}<0)}c_{v}^{\dagger}\right)\ket{\mathrm{vac}}_{e}\otimes\ket{\mathrm{FM}}\otimes\ket{\mathrm{FB}}_{R}.$
(56)

Here, $c_{v}^{\dagger}$ creates an electron in a single-particle eigenstate of
$H_{\mathrm{FM}}$. To calculate the susceptibility, we compute the action of
$S_{-\mathbf{q}}^{-}$ on this state:

$\begin{split}S_{-\mathbf{q}}^{-}&\ket{\mathrm{hf},\mathrm{FM}}=\\\
&\sum_{j}e^{-i\mathbf{q}\cdot\mathbf{R}_{j}}\left(\prod_{v\,(\varepsilon_{v}<0)}c_{v}^{\dagger}\right)\ket{\mathrm{vac}}_{e}\otimes\ket{\mathbf{R}_{j}}\otimes\ket{\mathrm{FB}}_{R}.\end{split}$
(57)

Unlike the case of the initial state $\ket{\mathrm{gs}}$, here the action of
$S_{-\mathbf{q}}^{-}$ does not simply generate a magnon excitation independent
of the state of the surface electrons. Instead, to rewrite Eq. (57) in terms
of $H_{S}$ eigenstates, the creation operators $c_{v}^{\dagger}$ should be
expanded, for each $j$, in terms of electronic states that account for the
position of the spin flip:

$c_{v}^{\dagger}=\sum_{\alpha}U_{v,\alpha j}c_{\alpha j}^{\dagger},$ (58)

with some coefficients $U_{v,\alpha j}$.

Since the single-particle electronic orbitals must adapt to the magnetic
configuration, the character of the half-filled state will also change,
depending on the site of the spin flip. In addition to producing the half-
filled state for each $j$, the spin flip will generally scatter surface
electrons to many excited states, resulting in a complex dynamics, different
from the case described above for the initial state $\ket{\mathrm{gs}}$.
However, the results obtained for the initial state $\ket{\mathrm{gs}}$ could
still provide insight into the qualitative behavior for the half-filled ground
state considered above, provided there is a mechanism for rapid equilibration
of the surface electrons. To illustrate this point, consider the
$N_{e}$-electron state

$\ket{\mathrm{hf}}_{\mathrm{FM}}\equiv\left(\prod_{v(\varepsilon_{v}<0)}c_{v}^{\dagger}\right)\ket{\mathrm{vac}}_{e}.$
(59)

Rewriting this in the eigenbasis of $H_{j}$:

$\ket{\mathrm{hf}}_{\mathrm{FM}}=\prescript{}{j}{\braket{\mathrm{hf}}{\mathrm{hf}}}_{\mathrm{FM}}\ket{\mathrm{hf}}_{j}+\sum_{\boldsymbol{\alpha}}\prescript{}{j}{\braket{\boldsymbol{\alpha}_{N_{e}}}{\mathrm{hf}}}_{\mathrm{FM}}\ket{\boldsymbol{\alpha}_{N_{e}}}_{j}.$
(60)

The lowest-energy $N_{e}$-particle eigenstate of $H_{j}$ is
$\ket{\mathrm{hf}}_{j}$ [see Eq. (23)]. If the surface electrons interact with
a dissipative environment (e.g. phonons) that quickly relaxes the electronic
state $\ket{\mathrm{hf}}_{\mathrm{FM}}$ to $\ket{\mathrm{hf}}_{j}$ in Eq.
(57), the state $S_{-\mathbf{q}}^{-}\ket{\mathrm{hf},\mathrm{FM}}$ will relax
into the state $\ket{\mathrm{hf},-\mathbf{q}}\otimes\ket{\mathrm{FB}}_{R}$
[see Eq. (26)]. This eigenstate of $H_{S}+H_{R}$ couples only to the continuum
of many-body JRMs $\\{\ket{\Psi,-\mathbf{q}}\otimes
a_{\nu}\ket{\mathrm{FB}}_{R}\\}$, similar to the empty-band case considered
above. We can then obtain the susceptibility, which will have the same
qualitative behavior as that given in Eq. (48). There would, however, be
corrections to the short-time dynamics of the spin-spin correlation function
due to relaxation effects, modifying the spectral function at high frequency.
A precise microscopic accounting of this relaxation mechanism and the
associated effect on the spectral function is, however, beyond the scope of
this work.

##  V Conclusions

In this paper, we have introduced a model that can describe Dirac electrons on
the surface of a 3DTI coupled to magnetic impurities, which themselves
interact via a ferromagnetic Heisenberg exchange. In contrast to previous
work, where the magnetic system was treated classically [4, 5, 15, 16, 29, 6,
24, 23], here we have considered a quantum magnet, which together with the
Dirac electron give rise to new excitations - the Jackiw-Rebbi-Magnons (JRMs).
We have systematically derived the form of these excitations in two
situations: First, we have solved a single-particle problem, in which the
conduction and valence bands of the 3DTI surface are initially empty, and
second, we have solved a many-body problem, starting from a half-filled 3DTI
surface. While the low-energy theory is well captured by a continuum
description of the Dirac electrons, we have also constructed a lattice model
with the appropriate low energy limit. This lattice model allowed us to
numerically study the density of electronic states, the creation of chiral
domain-wall bound states, and the associated current and magnetic textures for
JRMs (shown in Fig. 2). Finally, we studied the dynamic transverse spin
susceptibility. In order to probe the JRM states, we coupled our system to an
electronic reservoir, allowing a charge to be transferred to the surface, thus
producing a JRM that can then contribute to the magnetic susceptibility. In
contrast with the conventional magnon response, which would be independent of
the reservoir energy, we find a response that disperses with the reservoir
energy due to the formation of magnon-JRM polaritons (Fig. 4). In fact, this
model shows a quantum phase transition at a critical value of the reservoir
energy, where the fully spin-polarized ground state will collapse in favor of
the production of magnon-JRM polaritons.

This paper is limited to the study of single JRMs that carry a single unit of
spin angular momentum and a single unit of charge. An interesting direction
for future work would be to consider an alternative class of JRMs involving
multiple spin flips, but still binding a single charge. Alternatively, it
would be interesting to consider multiply charged JRMs. In each case, new
interaction effects should be taken into account and these effects may lead to
a rich family of excitations. In addition to the conventional magnon-magnon
interactions present in a Heisenberg ferromagnet, JRMs produced from multiple
spin flips would be influenced by the molecular binding energy of a
delocalized bound electron, similar to the situation that may arise for
classical magnetic skyrmions on the surface of a 3DTI [29]. For multi-electron
JRMs (or for multiple JRMs, each with a single electronic charge), it will be
important to account for the Coulomb interaction, in addition to magnetic
interactions and the molecular binding energy. Furthermore, we only consider a
specific type of magnetic system: a Heisenberg ferromagnet in a Bravais
lattice with nearest-neighbor interactions. There is, however, a rich family
of magnetic systems where electronic Jackiw-Rebbi states could bind to
magnons, including ferromagnetic materials with next-nearest-neighbor
interactions, systems with multiple single-magnon bands coming from sublattice
degrees of freedom (on, e.g., a honeycomb lattice), and antiferromagnetic
materials. Extending the present analysis to this broader range of magnetic
systems would be another interesting direction for future investigation.

###### Acknowledgements.

IM-B acknowledges a fellowship from Consejo Nacional de Ciencia y Tecnología
(CONACYT) and Mitacs. TP-B and WAC acknowledge support from Natural Sciences
and Engineering Research Council of Canada (NSERC), Institut
Transdisciplinaire d’Information Quantique (INTRIQ), and Fonds de Recherche du
Québec–Nature et Technologies (FRQNT).

## References

  * Hasan and Kane [2010] M. Z. Hasan and C. L. Kane, Rev. Mod. Phys. 82, 3045 (2010).
  * Qi and Zhang [2011] X.-L. Qi and S.-C. Zhang, Rev. Mod. Phys. 83, 1057 (2011).
  * Fu  _et al._ [2007] L. Fu, C. L. Kane, and E. J. Mele, Physical review letters 98, 106803 (2007).
  * Liu  _et al._ [2009] Q. Liu, C.-X. Liu, C. Xu, X.-L. Qi, and S.-C. Zhang, Phys. Rev. Lett. 102, 156603 (2009).
  * Yu  _et al._ [2010] R. Yu, W. Zhang, H.-J. Zhang, S.-C. Zhang, X. Dai, and Z. Fang, Science 329, 61 (2010).
  * Tokura  _et al._ [2019] Y. Tokura, K. Yasuda, and A. Tsukazaki, Nat. Rev. Phys. 1, 126 (2019).
  * Lee  _et al._ [2015] I. Lee, C. K. Kim, J. Lee, S. J. Billinge, R. Zhong, J. A. Schneeloch, T. Liu, T. Valla, J. M. Tranquada, G. Gu, and J. C. Davis, Proc. Natl. Acad. Sci. USA 112, 1316 (2015).
  * Chen  _et al._ [2010] Y. L. Chen, J.-H. Chu, J. G. Analytis, Z. K. Liu, K. Igarashi, H.-H. Kuo, X. L. Qi, S. K. Mo, R. G. Moore, D. H. Lu, M. Hashimoto, T. Sasagawa, S. C. Zhang, I. R. Fisher, Z. Hussain, and Z. X. Shen, Science 329, 659 (2010).
  * Xu  _et al._ [2012] S. Y. Xu, M. Neupane, C. Liu, D. Zhang, A. Richardella, L. A. Wray, N. Alidoust, M. Leandersson, T. Balasubramanian, J. Sánchez-Barriga, O. Rader, G. Landolt, B. Slomski, J. H. Dil, J. Osterwalder, T. R. Chang, H. T. Jeng, H. Lin, A. Bansil, N. Samarth, and M. Z. Hasan, Nat. Phys. 8, 616 (2012).
  * Chang  _et al._ [2013] C.-Z. Chang, J. Zhang, X. Feng, J. Shen, Z. Zhang, M. Guo, K. Li, Y. Ou, P. Wei, L.-L. Wang, _et al._ , Science 340, 167 (2013).
  * Eremeev  _et al._ [2013] S. Eremeev, V. Men’Shov, V. Tugushev, P. M. Echenique, and E. V. Chulkov, Phys. Rev. B 88, 144430 (2013).
  * Vobornik  _et al._ [2011] I. Vobornik, U. Manju, J. Fujii, F. Borgatti, P. Torelli, D. Krizmancic, Y. S. Hor, R. J. Cava, and G. Panaccione, Nano Lett. 11, 4079 (2011).
  * Watanabe  _et al._ [2019] R. Watanabe, R. Yoshimi, M. Kawamura, M. Mogi, A. Tsukazaki, X. Yu, K. Nakajima, K. S. Takahashi, M. Kawasaki, and Y. Tokura, Appl. Phys. Lett. 115, 102403 (2019).
  * Abanin and Pesin [2011] D. Abanin and D. Pesin, Phys. Rev. Lett. 106, 136802 (2011).
  * Hurst  _et al._ [2015] H. M. Hurst, D. K. Efimkin, J. Zang, and V. Galitski, Phys. Rev. B 91, 060401 (2015).
  * Tiwari  _et al._ [2017] K. L. Tiwari, W. A. Coish, and T. Pereg-Barnea, Phys. Rev. B 96, 235120 (2017).
  * Jackiw and Rebbi [1976] R. Jackiw and C. Rebbi, Phys. Rev. D 13, 3398 (1976).
  * Yasuda  _et al._ [2017] K. Yasuda, M. Mogi, R. Yoshimi, A. Tsukazaki, K. Takahashi, M. Kawasaki, F. Kagawa, and Y. Tokura, Science 358, 1311 (2017).
  * Nakajima  _et al._ [2015] Y. Nakajima, P. Syers, X. Wang, R. Wang, and J. Paglione, Nat. Phys. 12, 213 (2015).
  * Checkelsky  _et al._ [2012] J. G. Checkelsky, J. Ye, Y. Onose, Y. Iwasa, and Y. Tokura, Nat. Phys. 8, 729 (2012).
  * Rosen  _et al._ [2017] I. T. Rosen, E. J. Fox, X. Kou, L. Pan, K. L. Wang, and D. Goldhaber-Gordon, npj Quantum Mater. 2, 69 (2017).
  * Wang  _et al._ [2016] W. Wang, C.-Z. Chang, J. S. Moodera, and W. Wu, npj Quantum Mater. 1, 16023 (2016).
  * Efimkin and Kargarian [2021] D. K. Efimkin and M. Kargarian, Phys. Rev. B 104, 075413 (2021).
  * Tserkovnyak and Loss [2012] Y. Tserkovnyak and D. Loss, Phys. Rev. Lett. 108, 187201 (2012).
  * Chiu  _et al._ [2016] C.-K. Chiu, J. C. Teo, A. P. Schnyder, and S. Ryu, Rev. Mod. Phys. 88, 035005 (2016).
  * Shen [2017] S.-Q. Shen,  _Topological insulators_ (Springer Singapore, 2017) 2nd edition.
  * Gurvitz [2017] S. Gurvitz, Front. Phys. 12, 1 (2017).
  * Gaveau and Schulman [1995] B. Gaveau and L. Schulman, J. Phys. A 28, 7359 (1995).
  * Tiwari  _et al._ [2019] K. L. Tiwari, J. Lavoie, T. Pereg-Barnea, and W. A. Coish, Phys. Rev. B 100, 125414 (2019).