# Defect bulk-boundary correspondence of topological skyrmion phases of matter

Shu-Wei Liu  Max Planck Institute for Chemical Physics of Solids, Nöthnitzer
Strasse 40, 01187 Dresden, Germany  Max Planck Institute for the Physics of
Complex Systems, Nöthnitzer Strasse 38, 01187 Dresden, Germany  Li-kun Shi
Max Planck Institute for the Physics of Complex Systems, Nöthnitzer Strasse
38, 01187 Dresden, Germany  Ashley M. Cook  Max Planck Institute for Chemical
Physics of Solids, Nöthnitzer Strasse 40, 01187 Dresden, Germany  Max Planck
Institute for the Physics of Complex Systems, Nöthnitzer Strasse 38, 01187
Dresden, Germany

###### Abstract

Unpaired Majorana zero-modes are central to topological quantum computation
schemes as building blocks of topological qubits, and are therefore under
intense experimental and theoretical investigation. Their generalizations to
parafermions and Fibonacci anyons are also of great interest, in particular
for universal quantum computation schemes. In this work, we find a different
generalization of Majorana zero-modes in effectively non-interacting systems,
which are zero-energy bound states that exhibit a cross structure—two
straight, perpendicular lines in the complex plane— composed of the complex
number entries of the zero-mode wavefunction on a lattice, rather than a
single straight line formed by complex number entries of the wavefunction on a
lattice as in the case of an unpaired Majorana zero-mode. These “cross” zero-
modes are realized for topological skyrmion phases under certain open boundary
conditions when their characteristic momentum-space spin textures trap
topological defects. They therefore serve as a second type of bulk-boundary
correspondence for the topological skyrmion phases. In the process of
characterizing this defect bulk-boundary correspondence, we develop recipes
for constructing physically-relevant model Hamiltonians for topological
skyrmion phases, efficient methods for computing the skyrmion number, and
introduce three-dimensional topological skyrmion phases into the literature.

## I Introduction

Discovery of the integer and fractional quantum Hall effects in the early
1980’s pushed condensed matter physics beyond the Landau paradigm Chiu  _et
al._ (2016) of symmetry-breaking phases of matter, towards additional
characterization of matter based on its topology. Topological phases of
matter, or those phases of matter with some observables linked to topological
invariants and therefore unaffected by sufficiently small perturbations, have
profound consequences. They realize quasiparticles useful for exceptionally
robust quantum computation schemes Lahtinen and Pachos (2017), and their
topology can manifest through bulk-boundary correspondence: a non-trivial
value for a topological invariant of a material bulk yields topologically-
protected, gapless boundary states Hasan and Kane (2010); Bernevig and Hughes
(2013); Qi and Zhang (2011) very promising for applications such as
spintronics He  _et al._ (2021) and topologically-protected quantum
computation Nayak  _et al._ (2008); Kitaev (2001).

Study of topological phases has greatly expanded since the discovery of the
quantum spin Hall insulator, a topological phase of matter realized without an
applied magnetic field or magnetic order, presenting the possibility of far
broader applications than topological phases studied earlier Bernevig  _et
al._ (2006). Efforts began in earnest to search for other previously-
unidentified topological phases of matter, leading to the three-dimensional
topological insulator protected by time-reversal symmetry Hasan and Moore
(2011), and later the ten-fold way classification scheme for topological
phases of matter Schnyder  _et al._ (2008a); Ryu  _et al._ (2010a). Much more
recently, efforts to expand the ten-fold way to include bosonic topological
phases Tiwari  _et al._ (2018), crystalline topological phases Fu (2011);
Hsieh  _et al._ (2012); Huang  _et al._ (2017); Song  _et al._ (2017);
Cornfeld and Carmeli (2021), the broader set of short-range and long-range
entangled topological phases Chen  _et al._ (2013a); Gong  _et al._ (2016),
topology of various quasiparticles Khanikaev  _et al._ (2012); Lan and Wen
(2014); Susstrunk2015; Vijay  _et al._ (2015); Dai  _et al._ (2020); Tokuda
_et al._ (2021) and driven Lindner  _et al._ (2011); Cayssol  _et al._ (2013);
Khemani  _et al._ (2016); Crowley  _et al._ (2019); Wintersperger  _et al._
(2020) and open Shen  _et al._ (2018); Gong  _et al._ (2018); Zheng and
Hofstetter (2018); Kawabata  _et al._ (2019) systems reveal the foundational
role of topology in condensed matter physics.

In all of these topological phases of matter, however, topology is associated
with the full set of degrees of freedom of a system. This reflects the
characterization of many topological phases based on projectors onto occupied
states Schnyder  _et al._ (2008a); Ryu  _et al._ (2010a). Recent work has
shown, however, that topological phases are also realized for subsets of the
degrees of freedom of the occupied states. Specifically, past work showed non-
trivial topology is also associated with just the spin degree of freedom (DOF)
of the ground-state in symmetry-protected topological phases even when the
ground-state is characterized by a larger set of degrees of freedom Cook
(2023), independent of the previously-known topology associated with the full
set of degrees of freedom of the ground-state.

These topological phases of the spin degree of freedom, known as topological
skyrmion phases of matter, serve as momentum-space duals of the widely-studied
skyrmion magnetic orders in real-space Rößler  _et al._ (2006); Mühlbauer  _et
al._ (2009); Romming  _et al._ (2013); Fert  _et al._ (2017). As these
skyrmions are realized in momentum-space and their symmetry-protection is
robust against crystal-field splitting, however, they are topological in a
stricter sense, exhibiting a bulk-boundary correspondence upon performing a
partial trace over all degrees of freedom of a system except spin. As these
topological phases are protected by a generalized particle-hole symmetry
present in centrosymmetric superconductors, such signatures may be important
for experiments and understanding of exotic mechanisms of superconductivity.
The consequences of this non-trivial topology are even more significant,
however, yielding novel topological phase transitions even in effectively non-
interacting systems while the symmetry protecting this spin topology is
preserved, as the topology is protected by a spin magnitude gap rather than an
energy gap. Thus, given that the topological skyrmion phases are
experimentally-relevant and generalize foundational concepts of topological
condensed matter physics, better understanding of these phases will yield
valuable insights.

In this work, we examine these topological skyrmion phases of matter,
revealing they exhibit a second kind of bulk-boundary correspondence due to
the momentum-space spin textures trapping defects, while also introducing the
three-dimensional topological skyrmion phases of matter into the literature.
These results agree with past work on topological classification of defects by
Teo and Kane Teo and Kane (2010), but there are surprising additional
consequences of topological defects in momentum-space. Past work shows
unpaired Majorana zero-modes may be realized for certain boundary conditions
when spin skyrmions are realized in momentum-space for the Hopf insulatorMoore
_et al._ (2008); Deng  _et al._ (2013); Liu  _et al._ (2017a); Yan  _et al._
(2017). We extend these results to the topological skyrmion phases of matter,
showing these momentum-space spin skyrmions realize more exotic topologically-
protected zero-modes. In identifying this second defect-based bulk-boundary
correspondence of the topological skyrmion phases of matter, we provide strong
evidence that the consequences of this physics, including the type-II
topological phase transition of the topological skyrmion phases, must be
explored in detail to better understand the foundations of topological
condensed matter physics.

The overview of the manuscript is as follows: in Sec. II, we review the
topological skyrmion phases and then expand on earlier results by presenting
methods for constructing toy models—including those corresponding to
Bogoliubov de Gennes Hamiltonians for superconductors—and highly-efficient
calculation of the topological invariant. We use these techniques to
characterize the models numerically, studying their phase diagrams and
momentum-space spin textures. In Sec. III, we then extend our discussion by
introducing three-dimensional topological skyrmion phases using the
construction methods of Sec. II. Using these methods, we construct toy models
for the three-dimensional skyrmion phase from the Hopf insulator Liu  _et al._
(2017a); Yan  _et al._ (2017). We then show these three-dimensional spin
skyrmions trap momentum-space defects, which realize remarkable
generalizations of unpaired Majorana zero-modes for certain boundary
conditions. We characterize these so-called cross zero-modes numerically and
analytically, showing how they relate to the unpaired Majorana zero-mode.
Finally, in Sec. IV, we provide some summarizing remarks and discuss possible
extensions for future research.

## II Chiral topological skyrmion phase in two dimensions

### II.1 Section overview

In this section, we first provide a brief overview of the topological skyrmion
phases and the standard method for computing the topological invariant
characterizing them. We then introduce the simplest, or “atomic”, symmetry-
allowed toy model Bloch Hamiltonians that can realize this non-trivial
topology. Based on these atomic toy models, we then present some general
design principles for constructing toy model Bloch Hamiltonians with four
bands that can realize the chiral phase. We then introduce a far more
efficient method for computing the topological invariant and discuss
situations where it is most useful. Finally, we compute phase diagrams to
characterize the topology of some canonical four-band toy models exhibiting
the chiral topological skyrmion phase.

### II.2 Review of the chiral topological skyrmion phase

The chiral topological skyrmion phase is a topological phase of matter
introduced in past work Cook (2023), which is characterized by formation of a
baby skyrmion in the ground state spin expectation value texture over the two-
dimensional Brillouin zone, even when systems have more degrees of freedom
than just spin. The phase is realized in Bloch Hamiltonians with an even
number of bands greater than two, and has previously been studied in a tight-
binding model relevant to superconducting Sr2RuO4. This non-trivial topology
results when the Bloch Hamiltonian is invariant under a generalized particle-
hole transformation corresponding to operator $\mathcal{C}^{\prime}$, which is
the product of the particle-hole operator $\mathcal{C}$ (with
$\mathcal{C}^{2}=-1$), and the spatial inversion operator, $\mathcal{I}$. This
symmetry ensures that the Hamiltonian serves as a mapping from the Brillouin
zone to a space acted on by the quotient $\mathrm{Sp}(2N)/\mathrm{U}(N)$,
where $\mathrm{Sp}(2N)$ is the compact symplectic group with $\mathrm{Sp}(2)$
($N=1$) isomorphic to the special unitary group $\mathrm{SU}(2)$. One physical
quantity acted upon by such a quotient is the spin angular momentum. As a
result, systems with $\mathcal{C}^{\prime}$ symmetry can realize topological
phases characterized by the skyrmion number $\mathcal{Q}$ computed as an
integral over the Brillouin zone of the normalized ground state spin
expectation value ${\boldsymbol{\hat{S}}}({\bf k})={\boldsymbol{S}}({\bf
k})/|{\boldsymbol{S}}({\bf k})|$, where ${\boldsymbol{S}}({\bf
k})=\sum_{n\in\text{filled}}\langle u_{n}({\bf k})|{\bf S}|u_{n}({\bf
k})\rangle$ and the spin operator $\bf S$ reads

$S_{x}=\tau_{3}\sigma_{1},\quad S_{y}=\tau_{0}\sigma_{2},\quad
S_{z}=\tau_{3}\sigma_{3}.$ (1)

where here $\\{\tau_{i}\\}$ ($\\{\sigma_{i}\\}$) ($i\in\\{0,1,2,3\\}$ is a set
of Pauli matrices acting in generalized particle-hole (spin-${1\over 2}$)
space. $\tau_{0}$ and $\sigma_{0}$ are $2\times 2$ identity matrices.

We may express the skyrmion number as

$\displaystyle\mathcal{Q}={1\over
4\pi}\int_{BZ}d^{2}k\langle\hat{\boldsymbol{S}}(\boldsymbol{k})\rangle\cdot\left(\partial_{k_{x}}\langle\hat{\boldsymbol{S}}(\boldsymbol{k})\rangle\times\partial_{k_{y}}\langle\hat{\boldsymbol{S}}(\boldsymbol{k})\rangle\right).$
(2)

The integrand here is one component of an object we call the skyrmion
curvature $\boldsymbol{\Omega}_{\mathcal{Q}}$, expressed as

$\displaystyle\boldsymbol{\Omega}_{\mathcal{Q}}=\\{\Omega_{\mathcal{Q},x},\Omega_{\mathcal{Q},y},\Omega_{\mathcal{Q},z}\\},$
(3)

where the $k$th component is defined as

$\displaystyle\Omega_{\mathcal{Q},k}=\epsilon_{ijk}\langle\hat{\boldsymbol{S}}(\boldsymbol{k})\rangle\cdot\left(\partial_{k_{i}}\langle\hat{\boldsymbol{S}}(\boldsymbol{k})\rangle\times\partial_{k_{j}}\langle\hat{\boldsymbol{S}}(\boldsymbol{k})\rangle\right).$
(4)

where $\epsilon_{ijk}$ is the Levi-Civita tensor and Einstein summation is
implicit.

### II.3 Methods for constructing toy model Hamiltonians realizing the 2D
chiral topological skyrmion phase

In order to construct toy models realizing the chiral topological skyrmion
phase, we require that they serve as mappings from the Brillouin zone to the
space of ground state spin expectation values corresponding to non-trivial
homotopy group $\pi_{2}\left(\mathrm{Sp}(2N)/\mathrm{U}(N)\right)=\mathbb{Z}$.
As discussed in earlier work, this requires the Hamiltonian be invariant under
a $\mathcal{C}^{\prime}$ transformation Liu  _et al._ (2017b). For a 2$N$-band
model, the $\mathcal{C}^{\prime}$ symmetry operates on fermion operators as
$\mathcal{C}^{\prime}f_{k}\mathcal{C}^{\prime-1}=Jf_{k}^{\dagger}$ where

$J=\begin{pmatrix}0&I_{N\times N}\\\ -I_{N\times N}&0\end{pmatrix}.$ (5)

The Hamiltonian matrix representation must therefore satisfy

$J^{-1}H(\textbf{k})J=-H^{*}(\textbf{k}).$ (6)

The form of the Hamiltonian for $N=2$ (corresponding to four bands) is
therefore restricted to include some or all of the ten checked terms in TABLE
1, expressed as Kronecker products of Pauli matrices
$\tau_{i}\sigma_{j}=\tau_{i}\otimes\sigma_{j}$:

$\begin{array}[]{|c|c|c|c|c|}\hline\cr&\sigma_{0}&\sigma_{1}&\sigma_{2}&\sigma_{3}\\\
\hline\cr\tau_{0}&\times&\times&\checkmark&\times\\\
\hline\cr\tau_{1}&\checkmark&\checkmark&\times&\checkmark\\\
\hline\cr\tau_{2}&\checkmark&\checkmark&\times&\checkmark\\\
\hline\cr\tau_{3}&\checkmark&\checkmark&\times&\checkmark\\\
\hline\cr\end{array}$ Table 1: Invariance of the Bloch Hamiltonian under the
generalized particle-hole transformation $\mathcal{C}^{\prime}$ constrains its
form. The 16 possible terms, each represented by a Kronecker product of Pauli
matrices $\tau_{i}\sigma_{j}=\tau_{i}\otimes\sigma_{j}$, are each represented
in a table, with the $\tau$ Pauli matrix of the term represented by a row in
the table and the $\sigma$ Pauli matrix of the term represented by a column.
Of the sixteen terms, ten are symmetry-allowed (corresponding to a
$\checkmark$ for the corresponding entry in the table) and six are symmetry-
forbidden (corresponding to $\times$ for the corresponding entry in the
table).

We further classify these ten terms based on how they affect the ground state
spin expectation value for half-filling. To simplify the problem, we first
consider each term individually, computing the spin expectation value for each
$\tau_{i}\sigma_{j}$ term, of the ten terms allowed by $\mathcal{C}^{\prime}$
symmetry, separately. We find that three of the ten terms yield non-zero
values for the $x$-, $y$\- or $z$-components of the spin when considered
alone: specifically, $\tau_{3}\sigma_{1}$ alone yields non-zero
$\langle\boldsymbol{S}_{x}\rangle$, $\tau_{0}\sigma_{2}$ alone yields non-zero
$\langle\boldsymbol{S}_{y}\rangle$, and $\tau_{3}\sigma_{3}$ alone yields non-
zero $\langle\boldsymbol{S}_{z}\rangle$. That is, there exists one term
$\tau_{\alpha}\sigma_{\beta}$ with four eigenvectors labeled by
$i\in\\{1,2,3,4\\}$ as $\\{\psi_{i,\alpha\beta}\\}$ satisfying the eigen
problem
$\tau_{\alpha}\sigma_{\beta}\psi_{i,\alpha\beta}=E_{i,\alpha\beta}\psi_{i,\alpha\beta}$
and eigenvalues $\\{E_{i,\alpha\beta}\\}$ satisfying $E_{1,\alpha\beta}\leq
E_{2,\alpha\beta}\leq E_{3,\alpha\beta}\leq E_{4,\alpha\beta}$, such that
$\langle\boldsymbol{S}_{\gamma}\rangle=\sum_{i\in\\{1,2\\}}\psi^{\dagger}_{i,\alpha\beta}\boldsymbol{S}_{\gamma}\psi^{\dagger}_{i,\alpha\beta}\neq
0$ for each $\gamma\in\\{x,y,z\\}$. As the chiral skyrmion phase requires
realization of a skyrmionic spin texture in the Brillouin zone and therefore
finite $x$-, $y$\- and $z$-components of the ground state spin texture
somewhere in the Brillouin zone, the simplest toy model for a non-trivial
chiral skyrmion phase must therefore include these three terms. This
furthermore suggests that the simplest way to construct various four-band
Bloch Hamiltonians restricted to these three terms is as

$\displaystyle\mathcal{H}(\boldsymbol{k})$
$\displaystyle=d_{1}(\boldsymbol{k})\tau_{3}\sigma_{1}+d_{2}(\boldsymbol{k})\tau_{0}\sigma_{2}+d_{3}(\boldsymbol{k})\tau_{3}\sigma_{3},$
(7)

where $\\{d_{i}(\boldsymbol{k})\\}$ are the momentum-dependent functions of a
two-band model for a Chern insulator expressed as
$\mathcal{H}=d_{1}(\boldsymbol{k})\sigma_{1}+d_{2}(\boldsymbol{k})\sigma_{2}+d_{3}(\boldsymbol{k})\sigma_{3}$.

The other seven symmetry-allowed terms may be classified depending on how they
alter the $x$-, $y$-, or $z$-components of the ground state spin texture.
Specifically, a spin component can be finite when one of these seven terms is
multiplied by the same constant as one of the three essential terms, or zero.
We can therefore construct TABLE 2, marking each of the additional seven terms
as yielding $F$ (finite) or $Z$ (zero) spin-component when it is multiplied by
the same prefactor as one of the three essential terms.

$\begin{array}[]{|c|c|c|c|c|c|c|c|}\hline\cr&\tau_{1}\sigma_{0}&\tau_{1}\sigma_{1}&\tau_{1}\sigma_{3}&\tau_{2}\sigma_{0}&\tau_{2}\sigma_{1}&\tau_{2}\sigma_{3}&\tau_{3}\sigma_{0}\\\
\hline\cr\tau_{3}\sigma_{1}&F&F&Z&F&F&Z&Z\\\
\hline\cr\tau_{0}\sigma_{2}&Z&F&F&Z&F&F&Z\\\
\hline\cr\tau_{3}\sigma_{3}&F&Z&F&F&Z&F&Z\\\ \hline\cr\end{array}$ Table 2: In
Eqn. (7), each component of the $d$-vector
$\boldsymbol{d}(\boldsymbol{k})=\langle
d_{1}(\boldsymbol{k}),d_{2}(\boldsymbol{k}),d_{3}(\boldsymbol{k})\rangle$
multiplies a single Kronecker product of Pauli matrices (
$\tau_{3}\sigma_{1}$, $\tau_{0}\sigma_{2}$, or $\tau_{3}\sigma_{3}$,
respectively). We may instead multiply each component of the $d$-vector by
multiple symmetry-allowed Kronecker product terms. At the simplest level, we
multiply instead by a sum of one of $\tau_{3}\sigma_{1}$,
$\tau_{0}\sigma_{2}$, or $\tau_{3}\sigma_{3}$ (which correspond to each row of
the table) and one of the remaining seven terms (which correspond to each
column of the table). When the row term and column term are each multiplied by
the same pre-factor, they yield either zero for one component of the ground
state spin vector, such that the magnitude of the ground state spin is zero
somewhere in the Brillouin zone (corresponding to $Z$ in the table), or finite
spin-component (corresponding to $F$ in the table). The latter scenario
ensures there is a choice of $d$-vector yielding finite ground-state spin
magnitude for each momentum $\boldsymbol{k}$ in the Brillouin zone, such that
the skyrmion number $\mathcal{Q}$ is well-defined.

This forms the basis for constructing toy model Hamiltonians realizing the 2D
chiral topological skyrmion phase, but does not exhaust all possibilities for
symmetry-allowed toy model Bloch Hamiltonians with four bands.

### II.4 Relation between the Chern number and the Skyrmion number in a
special model

Here we investigate the relation between the Chern number and the skyrmion
number in a special case.

The Hall conductivity in an ideal lattice evaluated directly from the Kubo
formula reads

$\displaystyle\sigma_{xy}(\omega)=\frac{e^{2}\hbar}{V}$
$\displaystyle\sum_{{\bf k},m\neq n}(f_{m{\bf k}}-f_{n{\bf k}})\times$
$\displaystyle\frac{{\rm Im}\,\big{[}\langle u_{m}({\bf k})|v_{x}|u_{n}({\bf
k})\rangle\langle u_{n}({\bf k})|v_{y}|u_{m}({\bf
k})\rangle\big{]}}{(\epsilon_{n{\bf k}}-\epsilon_{m{\bf k}})(\epsilon_{n{\bf
k}}-\epsilon_{m{\bf k}}-\hbar\omega-i\eta)},$ (8)

where $m,n$ are band indices, $\omega$ and $\eta$ are set to zero in the dc
clean limit, the velocity operator $v_{a}=\partial H({\bf k})/\hbar\partial
k_{a}$ $(a=x,y)$ in Bloch basis, and $V$ is the volume of the crystal in the
real space.

For a $4$-band model described by $H({\bf k})=\sum_{\alpha}d_{\alpha}({\bf
k})\Gamma^{\alpha}$, we can express the velocity operator explicitly as
$v_{a}=\sum_{\alpha}[\partial_{k_{a}}d_{\alpha}({\bf
k})]\Gamma^{\alpha}/\hbar$. The hermicity of $H({\bf k})$ mandates that
$\Gamma^{\alpha}\in\\{\tau_{i}\otimes\sigma_{j}~{}|~{}0\leq i,j\leq 3\\}$. And
it follows that the matrix element for the velocity operator reads

$\langle u_{n}({\bf k})|v_{a}|u_{m}({\bf
k})\rangle=\sum_{\alpha}[\partial_{a}d_{\alpha}({\bf
k})]\langle\Gamma^{\alpha}({\bf k})\rangle_{nm}/\hbar,\quad(a=x,y)$ (9)

where $\langle\Gamma^{\alpha}({\bf k})\rangle_{nm}\equiv\langle u_{n}({\bf
k})|\Gamma^{\alpha}|u_{m}({\bf k})\rangle$. Using the explicit form of
velocity matrix elements, the dc Hall conductivity reads

$\displaystyle\sigma_{xy}=\frac{e^{2}}{h}$
$\displaystyle\Big{[}\frac{1}{2\pi}\int_{\rm BZ}{\rm d}{\bf k}\sum_{m\neq
n}\sum_{\alpha\beta}(f_{m{\bf k}}-f_{n{\bf k}})\times$
$\displaystyle[\partial_{x}d_{\alpha}({\bf k})]\,[\partial_{y}d_{\beta}({\bf
k})]\frac{{\rm Im}\,\big{[}\langle\Gamma^{\alpha}({\bf
k})\rangle_{mn}\langle\Gamma^{\beta}({\bf
k})\rangle_{nm}\big{]}}{(\epsilon_{n{\bf k}}-\epsilon_{m{\bf k}})^{2}}\Big{]}$
$\displaystyle=\frac{e^{2}}{h}$ $\displaystyle{\cal C},$ (10)

with $\cal C$ the well-known Chern number. We note that the second equation is
valid only when the chemical potential is inside the gap.

Following Eqn. (2) as the expression for the skyrmion number and focusing on
the Hamiltonian Eqn. (7), we have the following properties:
(1)

$\displaystyle\langle u_{n}({\bf k})|{\bf S}|u_{n}({\bf k})\rangle=-{\bf
d}({\bf k})/d({\bf k}),\hskip 5.69054ptn\in\\{1,2\\}$ $\displaystyle\langle
u_{n}({\bf k})|{\bf S}|u_{n}({\bf k})\rangle=+{\bf d}({\bf k})/d({\bf
k}),\hskip 5.69054ptn\in\\{3,4\\}$ (11)

where ${\bf d}({\bf k})=\big{(}d_{1}({\bf k}),d_{2}({\bf k}),d_{3}({\bf
k})\big{)}$, $d({\bf k})=|{\bf d}({\bf k})|$, and $n=1,2,3,4$ refers to the
band index;
(2)

$\epsilon_{1{\bf k}}=\epsilon_{2{\bf k}}=-d({\bf k}),\quad\epsilon_{3{\bf
k}}=\epsilon_{4{\bf k}}=+d({\bf k});$ (12)

(3)

$\sum_{n=3,4}{\rm Im}\,\big{[}\langle S_{a}({\bf k})\rangle_{mn}\langle
S_{b}({\bf k})\rangle_{nm}\big{]}=-\epsilon_{abc}\frac{d_{c}({\bf k})}{d({\bf
k})},\quad(m=1,2)$ (13)

and

$\sum_{n=1,2}{\rm Im}\,\big{[}\langle S_{a}({\bf k})\rangle_{mn}\langle
S_{b}({\bf k})\rangle_{nm}\big{]}=+\epsilon_{abc}\frac{d_{c}({\bf k})}{d({\bf
k})}.\quad(m=3,4)$ (14)

With the above three properties (1), (2), and (3), at half filling (bands
$1,2$ are filled) we can rewrite the Chern number as:

$\displaystyle{\cal C}$ $\displaystyle=\frac{1}{2\pi}\int_{\rm BZ}{\rm d}{\bf
k}\,4\sum_{ab}[\partial_{x}d_{a}({\bf k})]\,[\partial_{y}d_{b}({\bf
k})]\frac{-\epsilon_{abc}d_{c}({\bf k})}{4\,d^{3}({\bf k})}$
$\displaystyle=-2\,\Big{[}\frac{1}{4\pi}\int_{\rm BZ}{\rm d}{\bf
k}\,{\boldsymbol{\hat{S}}}({\bf
k})\cdot\Big{(}\partial_{k_{x}}{\boldsymbol{\hat{S}}}({\bf
k})\times\partial_{k_{y}}{\boldsymbol{\hat{S}}}({\bf k})\Big{)}\Big{]}$
$\displaystyle=-2{\cal Q}.$ (15)

### II.5 Efficient method for computing two-dimensional skyrmion number using
skyrmion connection

While the skyrmion number may be computed using the expression given in Eq.
(2), we may also compute it using a far more efficient method based on
projectors, developed for computing the Chern number in past work Fukui  _et
al._ (2005); Chen  _et al._ (2013b). The idea of this earlier work was that
computing a Berry phase directly from summation of Berry curvature is more
costly than computing it as a sum of infinitesimal Berry phases, each computed
using a discretized form of the Berry connection. These past works found that
the latter approach converged to the correct, integer-quantized values for
Chern numbers using a far coarser mesh of $\boldsymbol{k}$-points in the
Brillouin zone than did the former. In analogy to this approach, we want to
compute the skyrmion number with a connection, rather than the skyrmion
curvature defined previously. We therefore first present an efficient method
for computing this connection starting from the ground-state spin expectation
value texture over the Brillouin zone.

To compute this skyrmion connection, we first construct an eigenstate of
effective two-band model from the numerically computed ground state spin
expectation value $\langle\hat{\boldsymbol{S}}(\boldsymbol{k})\rangle$ of our
four-band model. Writing
$\langle\hat{\boldsymbol{S}}(\boldsymbol{k})\rangle=\\{\langle\hat{\boldsymbol{S}}_{x}(\boldsymbol{k})\rangle,\langle\hat{\boldsymbol{S}}_{y}(\boldsymbol{k})\rangle,\langle\hat{\boldsymbol{S}}_{z}(\boldsymbol{k})\rangle\\}=\\{d_{S,x}(\boldsymbol{k}),d_{S,y}(\boldsymbol{k}),d_{S,z}(\boldsymbol{k})\\}$,
we may write the effective two-band model $\mathcal{H}_{\mathrm{eff}}$ as

$\displaystyle\mathcal{H}_{\mathrm{eff}}$
$\displaystyle=d_{S,x}(\boldsymbol{k})\sigma_{x}+d_{S,y}(\boldsymbol{k})\sigma_{y}+d_{S,z}(\boldsymbol{k})\sigma_{z}.$
(16)

Assuming half-filling of this effective two-band model, we may then write the
form of the occupied eigenstate $|\psi_{-}(\boldsymbol{k})\rangle$
analytically in terms of $\\{d_{S,i}\\}$ as

$\displaystyle|\psi_{-}(\boldsymbol{k})\rangle$
$\displaystyle=\begin{pmatrix}\sin\frac{\theta}{2}e^{-i\phi}\\\
-\cos\frac{\theta}{2}\end{pmatrix}$ (17)

where $\theta=\arccos\frac{d_{S,z}(\boldsymbol{k})}{d(\boldsymbol{k})}$ and
$\phi=\arctan\frac{d_{S,x}(\boldsymbol{k})}{d_{S,y}(\boldsymbol{k})}$.

$|\psi_{-}(\boldsymbol{k})\rangle$ may then be used to compute the skyrmion
number as the Chern number for the lower band $\mathcal{C}_{-}$ of
$\mathcal{H}_{\mathrm{eff}}$, given that the skyrmion number and Chern number
are equal in two-band models at half-filling.

We evaluate $\mathcal{C}_{-}$ using a discrete lattice corresponding to the
discrete values of the crystal momentum in a finite sample with periodic
boundary conditions. The Chern number of the lower band may then be evaluated
as

$\mathcal{C}_{-}={1\over
2\pi}\sum_{\boldsymbol{k}}\mathrm{Im}\ln\left(A^{n}_{\boldsymbol{k},\hat{x}}A^{n}_{\boldsymbol{k}+\hat{x},\hat{y}}A^{n}_{\boldsymbol{k}+\hat{x}+\hat{y}}A^{n}_{\boldsymbol{k}+\hat{y},-\hat{y}}\right),$
(18)

where
$A^{n}_{\boldsymbol{k},\mu}=\langle\psi_{-}(\boldsymbol{k})|\psi_{-}(\boldsymbol{k}+\mu\rangle$
is the Berry connection field of the two-band model and the $\\{\mu\\}$ are
the nearest-neighbor vectors of the square momentum-space lattice. In the
context of this effective two-band model, we call this the $i$th component of
the skyrmion connection $\boldsymbol{A}_{\mathcal{Q}}$, given that this
effective two-band model is constructed from the ground state spin texture of
a model for a chiral topological skyrmion insulator phase with more than two
bands.

This connection is more generally expressed, for the $-$ band, as

$\displaystyle\boldsymbol{A}_{\mathcal{Q}}=\\{A_{\mathcal{Q},x},A_{\mathcal{Q},y},A_{\mathcal{Q},z}\\},$
(19)

with the $i$th component expressed as

$\displaystyle A_{\mathcal{Q},i}=-i\langle\psi_{-}|{\partial\over\partial
k_{i}}|\psi_{-}\rangle.$ (20)

The corresponding skyrmion curvature, which is the curl of this skyrmion
connection, is the integrand of Eq. 3.

While this approach of relating the skyrmion number in a $2N$-band model to
the Chern number of a two-band model to make use of additional computation
methods is generally not needed at the four-band level for skyrmion textures
with topological charge sufficiently small in magnitude, it becomes invaluable
for higher-magnitude topological charges or skyrmion numbers in models with
more bands. As we show later, this effective connection may be used to compute
the three-dimensional skyrmion number as well.

### II.6 Characterization of canonical toy models

Here, we construct some example toy models for $N=2$ chiral topological
skyrmion phases in two dimensions. We first present some example Hamiltonians
to illustrate the decoupling of the total Chern number and skyrmion number at
the four-band level by computing phase diagrams. We then discuss specific
signatures of non-trivial skyrmion number at the four-band level beyond the
momentum-space skyrmionic spin texture itself: a contribution to the Hall
conductivity that may be expressed in terms of a quantized skyrmion number in
$\mathcal{C}^{\prime}$-invariant models, and finite spin current signatures
required for a system to possess non-trivial skyrmion number.

We first consider a model with quantized total Chern number and quantized
skyrmion number, corresponding to co-existing stable Chern insulator and
chiral topological skyrmion insulator phases, by restricting our Hamiltonian
to only the three special terms in TABLE 2. We then show that a model
constructed using three of the seven other terms in the table may realize a
Chern insulator phase characterized by a stable non-trivial total Chern number
at half-filling but be topologically unstable according to the skyrmion
number. After this, we consider toy models for superconductors written in
Bogoliubov de Gennes (BdG) formalism, discussing constraints on the
Hamiltonian required for it to describe topological skyrmion phases of matter
in this physical setting and consider a representative case with spin singlet
pairing and another case with spin-triplet pairing.

#### II.6.1 Toy model with coexisting Chern insulator and skyrmion insulator
phases

Figure 1: Phase diagrams for Hamiltonian Eq. (21) at half-filling as a
function of free parameters $\kappa$ and $\beta$ for the (a) total Chern
number $\mathcal{C}$, (b) skyrmion number $\mathcal{Q}$, (c) minimum direct
energy gap between the second-lowest and third-lowest bands, $d_{2-1}$, and
(d) minimum magnitude of the ground-state spin expectation value $|\langle
S\rangle|$.

As the three special terms mentioned above can each support a non-zero spin
component, we first consider the following Bloch Hamiltonian:

$\displaystyle\mathcal{H}(\boldsymbol{k})$
$\displaystyle=[1-\cos(k_{x}a)-\cos(k_{y}a)]\tau_{3}\sigma_{3}+\beta\sin(k_{x}a)\tau_{0}\sigma_{2}$
$\displaystyle+\kappa\sin(k_{y}a)\tau_{3}\sigma_{1}.$ (21)

Figure 2: Characterization of a representative point in the phase diagrams for
Bloch Hamiltonian Eq. (21) shown in Fig. 1 ($\beta=\kappa=1$) with non-trivial
skyrmion number $\mathcal{Q}$. Minimum direct gap between the bands second-
lowest and third-lowest in energy over the first Brillouin zone is shown in
(a). $x$, $y$, and $z$ components of the normalized ground-state spin
expectation value, $\langle S_{x}\rangle$, $\langle S_{y}\rangle$, and
$\langle S_{z}\rangle$, over the first Brillouin zone are shown in (b), (c),
and (d), respectively, and correspond to winding of the spin in agreement with
the non-trivial skyrmion number value shown in Fig. 1 (b).

To characterize the topology of this model, we compute the total Chern number
$\mathcal{C}$ and skyrmion number $\mathcal{Q}$ assuming half filling of the
bands, as well as the minimum direct gap between the bands second-lowest and
third-lowest in energy over the Brillouin zone, $d_{2-1}$, and the minimum
magnitude of the ground state spin expectation value over the Brillouin zone
$|\langle\boldsymbol{S}\rangle|_{min}$, to obtain the phase diagrams shown in
Fig. 1. The total Chern number $\mathcal{C}$, shown in Fig. 1(a), and the
skyrmion number $\mathcal{Q}$, shown in Fig. 1(b), are each integer-quantized.
We find a relationship between them of $\mathcal{C}=-2\mathcal{Q}$, in
agreement with Eqn. (15) for the atomic model. The minimum direct gap between
the second-lowest and third-lowest bands in energy, $d_{2-1}$, is shown in
Fig. 1(c), and indicates that all topological phase transitions are type-I:
they occur with a closing of a direct gap. Type-II topological phase
transitions are not expected at the four-band level, due to restrictions on
DOFs of the Bloch Hamiltonian other than the spin-$1/2$ DOF. We also find the
minimum magnitude of the ground state spin expectation value over the
Brillouin zone is finite everywhere in the phase diagram, as shown in Fig.
1(d), as required for well-defined and integer-quantized $\mathcal{Q}$.

We examine a point in the phase diagram corresponding to $\beta=\kappa=1$ in
greater detail, computing the minimum direct gap between the second-lowest
band in energy and the third-lowest band in energy, $d_{2-1}$, over the
Brillouin zone, as well as the momentum-space spin texture, as shown in Fig.
2. $d_{2-1}$ over the Brillouin zone is shown in Fig. 2 (a) and confirmed to
be finite everywhere the Brillouin zone. It also reflects the four-fold
rotational symmetry of the model. Each of the three components of the spin
expectation value are shown in Fig. 2(b), (c) and (d) respectively. They
collectively show the winding of the spin vector in momentum space
corresponding to a baby skyrmion.

#### II.6.2 Toy model with well-defined total Chern number and ill-defined
skyrmion number

To further demonstrate that the total Chern number and skyrmion number are
decoupled at the four-band level, we consider a model constructed from a two-
band Chern insulator Hamiltonian, which has zero magnitude for the ground
state spin expectation value everywhere in the Brillouin zone. We choose the
following Bloch Hamiltonian:

$\displaystyle\mathcal{H}(\boldsymbol{k})$
$\displaystyle=\left[1-\cos(k_{x}a)-\cos(k_{y}a)\right]\tau_{3}\sigma_{3}+\beta\sin(k_{x}a)\tau_{1}\sigma_{3}$
$\displaystyle+\kappa\sin(k_{y}a)\tau_{2}\sigma_{0},$ (22)

and characterize it at half-filling by computing phase diagrams for the total
Chern number $\mathcal{C}$, skyrmion number $\mathcal{Q}$, minimum direct gap
between the bands second and third-lowest in energy, $d_{2-1}$, and the
minimum ground state spin magnitude over the Brillouin zone, $|\langle
S\rangle|$ shown in Fig. 3.

Figure 3:  Phase diagrams for Hamiltonian Eq. (22) at half-filling as a
function of free parameters $\kappa$ and $\beta$ for the (a) total Chern
number $\mathcal{C}$, (b) skyrmion number $\mathcal{Q}$, (c) minimum direct
energy gap between the second-lowest and third-lowest bands, $d_{2-1}$, and
(d) minimum magnitude of the ground-state spin expectation value $|\langle
S\rangle|$.

The phase diagram for the total Chern number of this Hamiltonian (Eq. (22)) in
Fig. 3 (a) is identical to that of the total Chern number for the Hamiltonian
Eq. (21). However, the skyrmion number is uniformly zero as shown in Fig.
3(b), because the $x$ and $y$ spin components do not wind, as shown in Fig. 3
(b) and (c), respectively. Additionally, the skyrmion number is actually ill-
defined, as the minimum magnitude of the ground state spin expectation value
over the Brillouin zone is zero as shown in Fig. 3(d). Specifically, the
normalized $x$\- and $y$-components of the spin are ill-defined in the
Brillouin zone as the unnormalized components are zero, corresponding to noise
shown in Fig. 4(b) and (c). The $z$-component of the ground state spin
expectation value also goes to zero at certain $\boldsymbol{k}$ points as
shown in Fig. 4(d), resulting in zero magnitude for the spin expectation value
at these momenta.

Figure 4: Characterization of a representative point in the phase diagrams for
Bloch Hamiltonian Eq. (22) shown in Fig. 3 ($\beta=\kappa=1$) with ill-defined
skyrmion number $\mathcal{Q}$. Minimum direct gap between the bands second-
lowest and third-lowest in energy over the first Brillouin zone is shown in
(a). $x$, $y$, and $z$ components of the normalized ground-state spin
expectation value, $\langle S_{x}\rangle$, $\langle S_{y}\rangle$, and
$\langle S_{z}\rangle$, over the first Brillouin zone are shown in (b), (c),
and (d), respectively, and correspond to unstable skyrmion number shown in
Fig. 3 (b) due to the minimum ground state spin magnitude being zero as shown
in Fig. 3 (c).

#### II.6.3 Toy models for chiral skyrmion superconductor phases

As the generalized particle-hole operator $\mathcal{C}^{\prime}$ is the
product of the particle-hole operator $\mathcal{C}$ and the spatial inversion
operator $\mathcal{I}$, we naturally expect the skymion phases in
centrosymmetric superconductors. Indeed, this has already been confirmed for
tight-binding models of Sr2RuO4 Ueno  _et al._ (2013). To better understand
when topological skyrmion phases are possible in such physical settings, we
first consider a generic four-band model for a superconductor written in BdG
formalism. The Hamiltonian $\mathcal{H}$ takes the form

$\mathcal{H}=\sum_{\boldsymbol{k}}\Psi^{\dagger}_{\boldsymbol{k}}\mathcal{H}(\boldsymbol{k})\Psi_{\boldsymbol{k}}/2,$
(23)

where

$\mathcal{H}(\boldsymbol{k})=\begin{pmatrix}\mathcal{E}(\boldsymbol{k})&\Delta(\boldsymbol{k})\\\
\Delta^{\dagger}(\boldsymbol{k})&-\mathcal{E}^{T}(-\boldsymbol{k})\\\
\end{pmatrix}$ (24)

and
$\Psi_{\boldsymbol{k}}=\left(c_{\boldsymbol{k}s},c_{\boldsymbol{k}s^{\prime}},c^{\dagger}_{-\boldsymbol{k}s},c^{\dagger}_{-\boldsymbol{k}s^{\prime}}\right)^{t}$.
Here, $c_{\boldsymbol{k}s}$ is the annihilation operator of electrons with
momentum $\boldsymbol{k}=\left(k_{x},k_{y}\right)$ and spin angular momentum
DOF $s\in\\{\uparrow,\downarrow\\}$. $\mathcal{E}(\boldsymbol{k})$ is the
Hamiltonian of the normal state, and $\Delta(\boldsymbol{k})$ is the gap
function of the superconductor.

This model possesses $\mathcal{C}^{\prime}$ symmetry when
$-\epsilon^{T}(-\boldsymbol{k})=-\epsilon^{*}(\boldsymbol{k})$, which is
satisfied for $\epsilon(\boldsymbol{k})$ even in momentum $\boldsymbol{k}$.
Therefore, to construct some representative example chiral topological
skyrmion superconductor Hamiltonains, we take the normal state Hamiltonian to
describe a Chern insulator that is even in momentum. One suitable example
model is the following, from Sticlet  _et al._ Sticlet  _et al._ (2012):

$\epsilon(\boldsymbol{k})=\alpha\left[\cos(k_{x})\sigma_{x}+\cos(k_{y})\sigma_{y}\right]+\beta\cos(k_{x}+k_{y})\sigma_{z},$
(25)

which realizes topologically-distinct insulating phases characterized by Chern
numbers of $0$, $\pm 2$ for the lower band. We may then combine this with a
pairing matrix $\Delta(\boldsymbol{k})$ describing either spin singlet pairing
$\Delta_{s}(\boldsymbol{k})$ or spin-triplet pairing
$\Delta_{t}(\boldsymbol{k})$. We take

$\displaystyle\Delta_{s}(\boldsymbol{k})$
$\displaystyle=i\Delta_{0}\sigma_{y}$ (26)
$\displaystyle\Delta_{t}(\boldsymbol{k})$
$\displaystyle=i\Delta_{0}\left(\boldsymbol{d}(\boldsymbol{k})\cdot\boldsymbol{\sigma}\right)\sigma_{y},$
(27)

with the $\boldsymbol{d}$-vector of the spin-triplet pairing term here chosen
specifically to be
$\boldsymbol{d}(\boldsymbol{k})=\sin(k_{y})\hat{x}-\sin(k_{x})\hat{y}$ as one
example. This $\boldsymbol{d}$-vector choice has previously been proposed as
characterizing Sr2RuO4 in the high-field phase Ueno  _et al._ (2013).

Figure 5: Phase diagrams for BdG Hamiltonian with normal state Eq. (25)(with
$\beta=1$)

and spin-singlet pairing given by gap function Eq. (26) at half-filling as a
function of pairing strength $\Delta_{0}$ and hopping integral $\alpha$ for
the (a) total Chern number $\mathcal{C}$, (b) skyrmion number $\mathcal{Q}$,
(c) minimum direct superconducting gap $d_{2-1}$, and (d) minimum magnitude of
the ground-state spin expectation value $|\langle S\rangle|$.

We may now compute phase diagrams for these two models characterizing spin-
singlet and spin-triplet superconductivity, respectively. We first consider
the case of spin-singlet pairing. To characterize the topology of the model,
we compute the total Chern number $\mathcal{C}$ and skyrmion number
$\mathcal{Q}$ assuming two bands are filled, as well as the minimum direct gap
between the third-lowest and second-lowest bands, $d_{2-1}$, and the minimum
magnitude of the ground state spin expectation value over the Brillouin zone,
$|\langle\boldsymbol{S}\rangle|$. Each quantity is computed as a function of
parameter $\alpha$ controlling some of the kinetic terms and the pairing
strength $\Delta_{0}$. These results for the case of spin-singlet pairing are
shown in Fig. 5.

The total Chern number $\mathcal{C}$ phase diagram shown in Fig. 5 (a) may be
divided into three distinct parts: a topologically stable region with total
Chern number of $4$ shown in red, which has an hourglass shape, two triangular
regions where the Chern number is well-quantized at zero shown in green, and
regions where the total Chern number is highly unstable. In comparison, the
phase diagram for the skyrmion number in Fig. 5 (b) possesses a tan hourglass
region with $\mathcal{Q}=2$ where the total Chern number $\mathcal{C}=4$. This
is further evidence that the skyrmion number differs from the total Chern
number by a factor of $2$ at the four band level when they are each
topologically stable. There are then two triangular regions with the skyrmion
number very well-quantized at zero, and regions around the hourglass where the
skyrmion number appears unquantized but close to zero.

Figure 6: Phase diagrams for BdG Hamiltonian with normal state Eq. (25)(with
$\beta=1$)

and spin-triplet pairing given by gap function Eq. (27) at half-filling as a
function of pairing strength $\Delta_{0}$ and hopping integral $\alpha$ for
the (a) total Chern number $\mathcal{C}$, (b) skyrmion number $\mathcal{Q}$,
(c) minimum direct superconducting gap $d_{2-1}$, and (d) minimum magnitude of
the ground-state spin expectation value $|\langle S\rangle|$.

To better understand the model, we also examine the phase diagram for the
minimum direct gap $d_{2-1}$, shown in Fig. 5 (c). We see that it is finite
only for the hourglass region, explaining the instability observed in Fig. 5
(a) and (b). In addition, we compute the phase diagram for the minimum
magnitude of the ground state spin expectation value,
$|\langle\boldsymbol{S}\rangle|$. We see that it is finite in the hourglass
region as required for a topologically-stable value of the skyrmion number,
but zero in the region where the skyrmion number is not well-quantized. This
is as expected, since the skyrmion number is computed with the normalized
ground state spin expectation value, and is therefore unstable when the spin
is not normalizable somewhere in the Brillouin zone.

We now consider the case of spin-triplet pairing. By tuning the parameters
$\alpha$ and $\Delta_{0}$, we again compute the total Chern number
$\mathcal{C}$ and skyrmion number $\mathcal{Q}$ assuming half-filling of
bands, as well as the minimum direct gap over the Brillouin zone $d_{2-1}$ and
the minimum ground state spin expectation value over the Brillouin zone
$|\langle\boldsymbol{S}\rangle|$.

The total Chern number phase diagram, shown in Fig. 6 (a), is topologically-
stable almost everywhere in the phase diagram at values $\mathcal{C}=\pm 4$.
The skyrmion number phase diagram, shown in Fig. 6 (b), is similar in form,
yet the tan regions with $\mathcal{Q}=2$ appear far more stable than the blue
regions with $\mathcal{Q}\approx-2$. We see, again, that topological phase
transitions occur where the minimum direct gap over the Brillouin zone between
the two middle bands, $d_{2-1}$, goes to zero, as shown in Fig. 6 (c). We gain
insight into the difference between the $\mathcal{Q}=2$ regions and the
$\mathcal{Q}=-2$ regions by examining the phase diagram for the minimum ground
state spin expectation value over the Brillouin zone,
$|\langle\boldsymbol{S}\rangle|$: regions with $\mathcal{Q}=2$ have finite
spin everywhere in the BZ, while $|\langle\boldsymbol{S}\rangle|$ is zero in
magnitude somewhere the BZ for $\mathcal{Q}=-2$ regions. Thus these regions
are topologically unstable according to the skyrmion number while being
topologically stable according to the total Chern number, further indicating
that these two topological invariants are decoupled from one another at the
four band level. Remarkably, however, the ground state spin texture nearly
forms skyrmions even in these cases where the spin is somewhere zero in
magnitude in the BZ, similar to results for the mirror subsectors of the
previously-considered tight-binding model for Sr2RuO4 Ueno  _et al._ (2013).

Figure 7: Characterization of a representative point in the phase diagrams for
BdG Hamiltonian with normal state Hamiltonian Eq. (25) and spin-singlet gap
function Eq. (26) shown in Fig. 5 ($\Delta_{0}=1$ and $\alpha=2$). Minimum
direct superconducting gap $d_{2-1}$ over the first Brillouin zone is shown in
(a). $x$, $y$, and $z$ components of the normalized ground-state spin
expectation value, $\langle S_{x}\rangle$, $\langle S_{y}\rangle$, and
$\langle S_{z}\rangle$, over the first Brillouin zone are shown in (b), (c),
and (d), respectively.

Noticeably, only type-I topological phase transitions occur in each model,
that is, only topological phase transitions that occur with a closing of the
direct gap between the two middle bands. This results from the four-band
$\mathcal{C}^{\prime}$-invariant Hamiltonians lacking a suitable DOF, such as
orbital angular momentum, into which spin angular momentum can flow. Although
one could interpret the Bloch Hamiltonian as corresponding to a system with a
spin half DOF and a two-fold orbital DOF, the $\mathcal{C}^{\prime}$
invariance of the model restricts the orbital DOF such that only the type-I
topological phase transition occurs at the four-band level. In Hamiltonians
with six or more bands, however, topology of an orbital DOF is no longer
ensured by $\mathcal{C}^{\prime}$ symmetry as discussed in work introducing
the topological skyrmion phases of matter and observed for tight-binding
models describing superconducting Sr2RuO4 Ueno  _et al._ (2013), making the
type-II topological phase transition possible.

We fix a particular choice of parameters $\Delta_{0}=1,\alpha=2$ to study the
spin textures in the $\mathcal{Q}=-2$ skyrmion phase. As shown in Fig. 7 and
Fig. 8, the system is always gapped and each component of the spin expectation
is always well-defined and varies continuously in the BZ. In comparison with
the $\mathcal{Q}=-1$ case in Fig. 2, the $x$\- and $y$-component of the spin
expectation demonstrate very similar behaviour, while the $z$-component winds
once around the BZ for $\mathcal{Q}=-1$ and twice for $\mathcal{Q}=-2$,
demonstrating the nature of $\mathcal{Q}$ as a winding number.

Figure 8:  Characterization of a representative point in the phase diagrams
for BdG Hamiltonian with normal state Hamiltonian Eq. (25) and spin-singlet
gap function Eq. (27) shown in Fig. 6 ($\Delta_{0}=1$ and $\alpha=2$). Minimum
direct superconducting gap $d_{2-1}$ over the first Brillouin zone is shown in
(a). $x$, $y$, and $z$ components of the normalized ground-state spin
expectation value, $\langle S_{x}\rangle$, $\langle S_{y}\rangle$, and
$\langle S_{z}\rangle$, over the first Brillouin zone are shown in (b), (c),
and (d), respectively.

## III Three-dimensional chiral topological skyrmion phase

### III.1 Section overview

In this second section of the paper, we introduce the three-dimensional
counterpart of the two-dimensional chiral topological skyrmion phase, taking
advantage of lessons gained in constructing suitable models in two-dimensions,
which are discussed in Section II.

We first explain why such a three-dimensional counterpart phase is
anticipated, before explaining how toy models for the 3D phase are constructed
taking inspiration from the topological phase known as the Hopf insulator
Moore  _et al._ (2008). To do this, we explain the defining characteristic of
the 3D chiral skyrmion phase and how analogous physics defines the Hopf
insulator phase, before utilizing techniques from Section I to construct a toy
model Hamiltonian for the 3D chiral topological skyrmion phase from a two-band
Hamiltonian for the Hopf insulator and a two-band Hamiltonian for its
$\mathcal{C}^{\prime}$ partner.

We then introduce the topological invariant for characterization of the 3D
chiral skyrmion phase, and discuss its numerical computation. We then use this
invariant to characterize some canonical toy models for the 3D chiral skyrmion
phase. Finally, motivated by past work on realization of unpaired Majorana
zero-modes in two-dimensional systems with a defect that are characterized by
the Hopf insulator Hamiltonian, we show that the Hamiltonian for the 3D chiral
skyrmion phase with the same boundary conditions also realizes topologically-
protected zero-modes.

### III.2 Construction of models for the three-dimensional chiral skyrmion
phase

While a three-dimensional counterpart of the chiral topological skyrmion phase
is expected in Hamiltonians with $\mathcal{C}^{\prime}$ symmetry, as they
correspond to the non-trivial homotopy group
$\pi_{3}\left(\mathrm{Sp}(2N)/\mathrm{U}(N)\right)=\mathbb{Z}_{2}$, it is not
immediately clear how to construct a model realizing this topological phase.

We expect the three-dimensional phase to correspond to formation of a hedgehog
in the ground state spin expectation value over the three-dimensional
Brillouin zone, rather than a baby skyrmion as in the two-dimensional case.
Given that the atomic model for the two-dimensional chiral skyrmion phase is
Eq. 7, which may be constructed from a two-band Chern insulator and its
$\mathcal{C}^{\prime}$ partner, we can see that the atomic model for the
three-dimensional case would be constructed from some two-band Hamiltonian
that realizes a momentum-space hedgehog and its $\mathcal{C}^{\prime}$
partner. There is such a two-band model: Hamiltonians realizing the Hopf
insulator Moore  _et al._ (2008) phase of matter realize such hedgehogs in the
spin or pseudospin texture over the Brillouin zone when the topological
invariant characterizing the phase, the Hopf invariant, is non-trivial.

Figure 9: Phase diagrams as a function of $M$ and $\Delta_{0}$ for the (a)
minimum direct energy gap between the second- and third-lowest bands in
energy, $d_{2-1}$, (b) minimum spin magnitude $|\langle S\rangle|$ over the
Brillouin zone, and (c) skyrmion number $\mathcal{Q}$. (d) shows the hedgehog
shape formed by the ground-state spin expectation value texture over the
Brillouin zone for Hamiltonian Eq. (32) plus spin singlet pairing Eq. (26) at
half-filling in a plot of the ground state spin texture for two spherical
shells of momenta centred at the origin in the 3D Brillouin zone with
$|\boldsymbol{k}|=\pi$ (red arrows) and $|\boldsymbol{k}|=0.75\pi$ (blue
arrows)

.

The two-band Bloch Hamiltonian $h(\boldsymbol{k})$ for the Hopf insulator
takes the following form:

$h(\boldsymbol{k})=\sum_{i=1,2,3}d_{i}(\boldsymbol{k})\tau_{i},$ (28)

where

$d_{i}=z^{\dagger}\tau_{i}z,$ (29)

for which we take $z=(z_{1},z_{2})^{T}$ where

$\displaystyle z_{1}$ $\displaystyle=\sin k_{x}+i\sin k_{y},$ (30)
$\displaystyle z_{2}$ $\displaystyle=\sin k_{z}+i(\cos k_{x}+\cos k_{y}+\cos
k_{z}-M).$ (31)

The four-band atomic model for the three-dimensional topological skyrmion
phase of matter may then be written as

$\mathcal{H}(\boldsymbol{k})=\begin{pmatrix}h(\boldsymbol{k})&0\\\
0&-h^{*}(\boldsymbol{k})\end{pmatrix}.$ (32)

As in two-dimensions, we can add additional, symmetry-allowed terms to this
Bloch Hamiltonian, such as spin singlet pairing Eq. (26) or spin triplet
pairing Eq. (27). Phase diagrams are computed for these Hamiltonians with $M$
in Eq. (31) and $\Delta_{0}$ in Eq. (27) as varying parameters, and the
results are shown in Fig. 9 and 10, respectively. We show that the system
realizes a skyrmion phase for a range of parameters $\Delta_{0}$ and $M$ as
characterized by the skyrmion number.

### III.3 Topological invariant in three-dimensions

Here, we introduce the topological invariant of the three dimensional
topological skyrmion phases of matter. To provide context, we first review the
Hopf invariant. Although there are some similarities between the form of the
Hopf invariant and the topological invariant of the 3D topological skyrmion
phase, they are computed with different physical quantities and are therefore
in general distinct.

The Hopf invariant is the topological charge of a vector field defined over
the 3D Brillouin zone, where the vector is the projector onto occupied states.
Following the paper introducing the Hopf insulator protected by
$\mathcal{C}^{\prime}$ symmetry in systems with $2$ or more bands Liu  _et
al._ (2017b), we consider a Hamiltonian $\mathcal{H}(\boldsymbol{k})$ which
takes the following form when diagonalized,

$\mathcal{H}(\boldsymbol{k})=\mathrm{diag}\left(\lambda_{1}(\boldsymbol{k}),-\lambda_{1}(\boldsymbol{k}),...,\lambda_{N}(\boldsymbol{k}),-\lambda_{N}(\boldsymbol{k})\right).$
(33)

Continuously deforming all of the positive eigenvalues to $+1$, and all of the
negative eigenvalues to $-1$, the deformed Hamiltonian takes the form,

$\mathcal{H}(\boldsymbol{k})=U(\boldsymbol{k})I_{N,N}U^{\dagger}(\boldsymbol{k}),$
(34)

where $U(\boldsymbol{k})$ is the matrix of eigenvectors and
$I_{N,N}=I_{N}\otimes\sigma_{z}$, where $I_{N}$ is the $N\times N$ identity
matrix and $\sigma_{z}$ is the Pauli matrix $\sigma_{z}=\mathrm{diag}(1,-1)$.
We can then define the Hopf invariant in terms of the projectors onto occupied
states as in Liu  _et al._ Liu  _et al._ (2017b) as

$n_{\omega}={1\over 24\pi^{2}}\int d^{3}k\left(U^{-1}dU\right)^{3},$ (35)

and therefore denote it as a projector topological invariant. If the only
degree of freedom labeling the fermionic creation and annihilation operators
in a given $\boldsymbol{k}$-sector is spin, this projector also corresponds to
the ground-state spin expectation value. If the operators are labeled not just
by spin, but also by other degrees of freedom, such as an orbital degree of
freedom, or a particle-hole degree of freedom, the projector onto occupied
states becomes distinct from the ground-state spin expectation value vector.
In this case, one can still compute the topological charge for the projector
vector field, but also a second topological charge, for the spin expectation
value vector. In some cases, this spin topological charge quantizes and serves
as a second topological invariant distinct from the projector-based
topological invariant, which we define as a spin topological invariant. The 3D
topological skyrmion phase corresponds to the system realizing non-trivial,
quantized spin topological charge, when the spin expectation value vector is
distinct from the projector onto occupied states.

For these cases where the projector topological invariant and spin topological
invariant are computed with distinct physical quantities, we may then write
the topological invariant of the 3D chiral skyrmion phase as

$\displaystyle\mathcal{Q}_{3D}=-{1\over
8\pi}\int_{BZ}d^{3}k\boldsymbol{A}_{\mathcal{Q}}\cdot\boldsymbol{\Omega}_{\mathcal{Q}}$
(36)

where here $\boldsymbol{A}_{\mathcal{Q}}$ is the skyrmion connection defined
in Eq. 19, and $\boldsymbol{\Omega}_{\mathcal{Q}}$ is the skyrmion curvature
defined in Eq. 3.

In this work, we restrict ourselves to four band models constructed from a
two-band Hamiltonian and its $\mathcal{C}^{\prime}$ partner, which each
realize a topological phase characterized by the Hopf invariant. We may
therefore compute a Hopf invariant using projectors onto occupied states as
defined above, but also a separate 3D skyrmion number.

The 3D skyrmion number $\mathcal{Q}_{3D}$ may therefore be computed using Eq.
4 for the skyrmion curvature, and Eq. 20 for the skyrmion connection,
computing the derivatives numerically. Using this invariant
$\mathcal{Q}_{3D}$, we characterize the topology of the toy models presented
for the 3D chiral skyrmion phase.

A similar set of phase diagrams is computed for the Hamiltonian in Eq. (32)
with free parameter $M$ plus additional spin triplet pairing term Eq. (27)
with free parameter $\Delta_{0}$. The results are shown in Fig. 9 and 10. We
show that the system can be in the three-dimensional topological skyrmion
phase for a range of parameters $\Delta_{0}$ and $M$ as characterized by the
skyrmion number.

Figure 10: Phase diagrams as a function of $M$ and $\Delta_{0}$ for the (a)
minimum direct energy gap between the second- and third-lowest bands in
energy, $d_{2-1}$, (b) minimum spin magnitude $|\langle S\rangle|$ over the
Brillouin zone, and (c) skyrmion number $\mathcal{Q}$. (d) shows the hedgehog
shape formed by the ground-state spin expectation value texture over the
Brillouin zone for Hamiltonian Eq. (32) plus spin triplet pairing Eq. (27) at
half-filling in a plot of the ground state spin texture for two spherical
shells of momenta centred at the origin in the 3D Brillouin zone with
$|\boldsymbol{k}|=\pi$ (red arrows) and $|\boldsymbol{k}|=0.75\pi$ (blue
arrows)

.

### III.4 Bulk-boundary correspondence

We now explore the bulk-boundary correspondence associated with non-trivial
$\mathcal{Q}_{3D}$. Bulk-boundary correspondence is expected when a projector
topological invariant is non-trivial, and when symmetries protecting the
topological phase in the bulk are respected in the bulk and by open boundary
conditions Schnyder  _et al._ (2008a); Ryu  _et al._ (2010b). Here, we instead
investigate whether the bulk-boundary correspondence is realized for a
topological phase characterized by trivial projector topological invariant,
the Hopf invariant $n_{\omega}$, and non-trivial spin topological invariant
$\mathcal{Q}_{3D}$.

To provide context, we first review the bulk-boundary correspondence of the
Hopf insulator, as there are at least two different kinds of bulk-boundary
correspondence, which have been considered for this system and are relevant to
the present work. These two types of bulk-boundary correspondence are
distinguished by different open boundary conditions. The first kind of bulk-
boundary correspondence of the Hopf insulator is discussed in Moore  _et al._
Moore  _et al._ (2008) introducing the Hopf insulator phase. This corresponds
to open boundary conditions in the $\hat{z}$-direction, and periodic boundary
conditions in the $\hat{x}$ and $\hat{y}$ directions. Under these boundary
conditions, non-trivial Hopf invariant corresponds to gapless boundary modes
known as Fermi rings localized on the top and bottom surfaces.

Here, it is important to distinguish between the Hopf insulator phase realized
in a two-band system, and the Hopf insulator realized in systems with more
than two bands protected by generalized particle-hole symmetry
$\mathcal{C}^{\prime}$ Liu  _et al._ (2017a). The two-band Hopf insulator has
$\mathbb{Z}$ topological classification and Fermi ring surface states due to
bulk-boundary correspondence. The $\mathcal{C}^{\prime}$-protected Hopf
insulator with more than two bands instead has $\mathbb{Z}_{2}$ topological
classification. Furthermore, as the $\mathcal{C}^{\prime}$ operation involves
particle-hole conjugation in combination with spatial inversion, opening
boundary conditions in the $\hat{z}$-direction generically breaks
$\mathcal{C}^{\prime}$ near the surfaces. Gapless surface states are then not
topologically-protected. As our 3D chiral skyrmion phase is protected by
$\mathcal{C}^{\prime}$ symmetry, we do not expect to find this kind of bulk-
boundary correspondence.

A second kind of bulk-boundary correspondence has also been considered in the
context of the Hopf insulator Yan and Felser (2017). In this case, boundary
conditions are first opened in the $\hat{x}$\- and $\hat{y}$-directions, while
periodic boundary conditions are retained in the $\hat{z}$-direction. The
Hamiltonian is therefore defined in terms of two real-space coordinates $x$
and $y$, and one momentum component $k_{z}$. Then, one performs the
substitution $k_{z}\rightarrow\theta(x,y)$. This spatially-varying $\theta$ is
interpreted as an angle in the $x-y$ plane, which can characterize a zero-
dimensional defect.

We first review the results of Yan  _et al_ Yan  _et al._ (2017) on this
second kind of bulk-boundary correspondence, in particular discussing the
nature of the gapless boundary modes. The boundary conditions are shown in
Fig. 11 (a) for the Hopf insulator lattice model. A corresponding annulus
geometry for the low-energy theory is shown in Fig. 11 (b). This geometry
yields a bulk-boundary correspondence when the Hopf invariant is non-trivial,
and the Hopf insulator Bloch Hamiltonian realizes a hedgehog in its
(pseudo)spin texture centered at the origin of the three-dimensional Brillouin
zone. According to the topological classification of defects by Teo and Kane
Teo and Kane (2010), the defect can correspond to a non-trivial homotopy group
and realize a topologically non-trivial, odd-parity state. As the defect is a
zero-dimensional system within the Brillouin zone, it is, however, not
immediately clear how to realize unpaired Majorana zero-modes from this odd-
parity state. A key insight of Yan  _et al._ Yan  _et al._ (2017) was that re-
interpreting the Hopf insulator Bloch Hamiltonian as describing a two-
dimensional system with a defect means the zero-dimensional defect of the
three-dimensional Brillouin zone corresponds to an effectively one-dimensional
system along the line $\theta=0$ in the 2D system with a defect. Then,
similarly to one-dimensional wire proposals for realizing unpaired Majorana
zero-modes, an unpaired Majorana zero-mode is realized at each end of this
$\theta=0$ wire.

Figure 11: (a) The two-dimensional lattice with the defect located at $(0,0)$
and $\theta=\arctan(y/x)$. (b) Schematic illustration of the annulus geometry
considered for low-energy analytical theory of the defect bulk-boundary
correspondence. For the Hopf insulator, at $\theta=0$, one unpaired Majorana
zero-mode is localized on the inner boundary and one on the outer boundary of
the annulus, respectively, as shown in this schematic.

For such a finite-size, two-dimensional system described by the Hopf insulator
Hamiltonian with $k_{z}$ replaced by angle $\theta$ characterizing a defect at
the center of the system, we first numerically compute the spectrum and
eigenstates of the Hopf insulator Bloch Hamiltonian as is done in Yan  _et al_
Yan  _et al._ (2017), with a defect located at $(0,0)$ and
$\theta=\arctan(y/x)$. The results are shown in Fig. 12 for a finite-size
system with $40\times 40$ lattice sites. These results include the low-energy
spectrum, the particle and hole component probability densities of a zero-
energy eigenstate, and the probability densities for each unpaired Majorana
zero-mode, shown in Fig. 12 (a), (b), (c), and (d), respectively. As expected,
the unpaired Majorana zero-mode at the boundary occurs for $\theta=0$,
corresponding to the location of the zero-dimensional defect in the three-
dimensional Brillouin zone of the Hopf insulator.

We then additionally characterize the resultant unpaired Majorana zero-modes
in greater detail. Specifically, after taking linear combinations of the two
lowest energy eigenstates of the Hopf insulator Hamiltonian with these
boundary conditions to compute arrays with $40\times 40\times 2$ entries
representing the unpaired Majorana zero-mode state wavefunction, we then plot
each of these entries in the complex plane as a point, to better understand
later results for the 3D chiral skyrmion phase. These representations in the
complex plane of the unpaired Majorana bound states localized at the defect
and boundary are shown in Fig. 12 (e) and (f), respectively. The MZM at the
boundary appears as a straight line in the complex plane. This is expected for
an unpaired MZM, which has a purely real wavefunction up to some phase degree
of freedom. The structure of the MZM at the defect, in contrast, is polluted.
This complicated structure of the zero-mode state localized at the defect is
unexpected for an unpaired MZM: a key feature of unpaired MZMs is that they
should be straight lines i.e. they are purely real up to a phase degree of
freedomAlicea (2012); Read and Green (2000). The unpaired MZM localized at the
defect likely loses this feature because of coupling to higher energy states
at $\theta\neq 0$ due to the small radius of the defect.

Given that the 3D chiral skyrmion phase toy model presented here is
constructed from the Hopf insulator Hamiltonian, we anticipate a similar bulk-
boundary correspondence for the 3D chiral skyrmion phase, when describing a
two-dimensional system with open boundary conditions and a defect, as occurs
for the Hopf insulator Hamiltonian itself. We consider the following 3D chiral
skyrmion phase Hamiltonian,

$\displaystyle
H_{\textnormal{3D}}(\boldsymbol{k})=\mathcal{H}(\boldsymbol{k})+\Delta_{t}(\boldsymbol{k})$
(37)

where $\mathcal{H}(\boldsymbol{k})$ is the 4-band Hopf Hamiltonian in Eq. (32)
and $\Delta_{t}(\boldsymbol{k})$ is the spin triplet pairing matrix in Eq.
(27).

For this model, we then consider open boundary conditions in the $\hat{x}$\-
and $\hat{y}$-directions, and substitute an angle $\theta$ for momentum
component $k_{z}$. The results for 3D chiral skyrmion phase Hamiltonian Eq.
(37) under these conditions, when the 3D skyrmion number is $1$, are shown in
Fig. 13 with additional disorder-averaging to illustrate that the results are
independent of the specific form of $\Delta_{t}(\boldsymbol{k})$. The on-site
disorder term takes the form of
$\epsilon\sim\tau_{3}\sigma_{0}\cdot\epsilon_{0}\textnormal{U}(-1,1)$ where
$\tau_{3}\sigma_{0}$ is a $C^{\prime}$ conserving term, $\epsilon_{0}$ is the
magnitude of the disorder which is 5% of the difference between non-zero
energy levels (i.e. the lowest positive level and the highest negative level),
and $\textnormal{U}(-1,1)$ is a uniform distribution with values from -1 to 1.
The results shown are obtained from averaging 100 realizations.

While two zero-energy states are realized for the Hopf insulator under these
boundary conditions as shown in Fig. 12(a), four zero-energy states are
realized in the low-energy spectrum under the same conditions for the 3D
topological skyrmion phase, as shown in Fig. 13(a). The robustness of these
states against local disorder in the latter case has been numerically verified
by introducing random on-site perturbation at each lattice site.

Fig. 13(b) shows the probability density profile of the particle component
(denoted by $u$) and hole component (denoted by $v$) of a zero-energy
eigenstate realized by the 3D topological skyrmion phase Hamiltonian,
respectively. This can be compared directly with Fig. 12(b). Notably, we see
that, unlike in the case of the Hopf insulator, the zero-modes of the 3D
topological skyrmion phase outputted by numerical diagonalization are
particle-hole symmetric, but also localized. Fig. 13(c)(d) each show one
localized zero-mode for direct comparison with Fig. 12(c)(d). These two
localized zero-mode wavefunctions at the defect and boundary are then
visualized in the complex plane in Fig. 13(e) and (f), respectively, similarly
to Fig. 12(e) and (f). That is, each wavefunction is represented by an array
of $40\times 40\times 4$ complex numbers, and these complex numbers are
plotted in the complex plane as individual points. Similarly to the unpaired
MZM at the defect in the Hopf insulator, the zero-modes localized at the
defect in the 3D topological skyrmion phase have complex structure. The zero-
modes of the 3D topological skyrmion phase localized at the boundary, however,
reveal an interesting difference between the case of the Hopf insulator and
the 3D topological skyrmion phase. They appear as crosses, or two straight,
perpendicular lines, rather than as a single straight line as for a MZM, in
the complex plane. In agreement with the low-energy theory presented in the
following section and in the supplemental materials for the zero-modes on the
boundary, linear combinations of these two-cross states yield two unpaired
Majorana zero-modes (one purely real and one purely imaginary), each localized
at the boundary. Rather than combine to form a generic complex fermion bound
state, however, they are forbidden from hybridizing by a combination of
$C^{\prime}$ and $C$ symmetry to instead yield the cross structure.

Figure 12: For the Hopf Hamiltonian given by Eq. (28) with open boundary
conditions in the $x$ and $y$ directions and the $z$ component of momentum
$k_{z}$ replaced with $\theta$ characterizing a defect, (a) spectrum vs.
eigenvalue index $i$ in the vicinity of zero-energy, showing a two-fold
degeneracy at zero-energy. (b) is the probability density plot of particle
component $u$ with the hole component $v$ in the inset of one of the two zero-
energy states outputted by numerical diagonalization. (c) and (d) show the
strongly localized unpaired Majorana zero-modes computed as linear
combinations of the two zero-energy states yielded by numerical
diagonalization. (e) and (f) show each complex number in the eigenvector for
each of the unpaired Majorana bounds plotted in the complex plane. Figure 13:
For the disorder-averaged four-band Bloch Hamiltonian characterizing the 3D
topological skyrmion phase and given by Eq. (32), constructed from the Hopf
insulator Bloch Hamiltonian, its $\mathcal{C}^{\prime}$ partner, and spin
triplet pairing gap function given by Eq. (27) with open boundary conditions
in the $x$ and $y$ directions and the $z$ component of momentum $k_{z}$
replaced with $\theta$ characterizing a defect, (a) spectrum vs. eigenvalue
index $i$ in the vicinity of zero-energy, showing a four-fold degeneracy at
zero-energy. (b), (c) and (d) show the probability density for one zero-mode
localized at the defect and one localized at the boundary, respectively. (e)
and (f) show each complex number in the eigenvector plotted in the complex
plane for a zero-energy state localized at the defect and at the boundary,
respectively.

### III.5 Analytical characterization of topologically-protected zero-modes

To further characterize the cross-zero modes of the 3D topological skyrmion
phase, we derive a low-energy effective theory for the model Hamiltonian Eq.
(37) with open boundary conditions corresponding to an annulus geometry as in
Yan  _et al._ Yan  _et al._ (2017) and shown schematically in Fig. 11. Full
details are presented in the Supplementary Materials, and here we focus on the
resultant low-energy effective theory.

For a semi-infinite geometry, with the sample occupying the $x>0$ region the
resultant effective Hamiltonian written in the on-site basis of the full
Hamiltonian (generalized particle-hole DOF and spin DOF) is

$\displaystyle H_{\textnormal{eff}}$
$\displaystyle=\begin{pmatrix}0&\frac{3}{4}k_{y}+\frac{3}{4}i\lambda&-\frac{3}{8}\Delta_{0}&0\\\
\frac{3}{4}k_{y}-\frac{3}{4}i\lambda&0&0&\frac{3}{8}\Delta_{0}\\\
-\frac{3}{8}\Delta_{0}&0&0&-\frac{3}{4}k_{y}+\frac{3}{4}i\lambda\\\
0&\frac{3}{8}\Delta_{0}&-\frac{3}{4}k_{y}-\frac{3}{4}i\lambda&0\end{pmatrix}$
(38)

By taking $k_{y}\rightarrow-i\partial_{y}$ and $\lambda\rightarrow y/R$,
$H_{\textnormal{eff}}$ becomes

$\displaystyle H_{\textnormal{eff}}$
$\displaystyle=\frac{3}{4}\begin{pmatrix}0&-i\partial_{y}+i\frac{1}{R}y&-\frac{1}{2}\Delta_{0}&0\\\
-i\partial_{y}-i\frac{1}{R}y&0&0&\frac{1}{2}\Delta_{0}\\\
-\frac{1}{2}\Delta_{0}&0&0&i\partial_{y}+i\frac{1}{R}y\\\
0&\frac{1}{2}\Delta_{0}&i\partial_{y}-i\frac{1}{R}y&0\end{pmatrix}.$ (39)

Squaring both sides of the time-independent Schrödinger equation yields the
resultant equation

$\displaystyle\frac{9}{16}$
$\displaystyle\begin{pmatrix}\epsilon_{y,-}&0&0&-i\partial_{y}\Delta_{0}\\\
0&\epsilon_{y,+}&i\partial_{y}\Delta_{0}&0\\\
0&i\partial_{y}\Delta_{0}&\epsilon_{y,+}&0\\\
-i\partial_{y}\Delta_{0}&0&0&\epsilon_{y,-}\end{pmatrix}\psi=E^{2}\psi,$ (40)

where
$\epsilon_{y,\pm}=-\partial_{y}^{2}+\frac{1}{R^{2}}y^{2}\pm\frac{1}{R}+\frac{1}{4}\Delta_{0}^{2}$.

The eigenstates of this low-energy effective Hamiltonian may then be expressed
as $|\mu_{i}\rangle=\chi_{i}^{\top}|\psi_{i}\rangle$, where $\\{\chi_{i}\\}$
($i\in\\{1,2,3,4\\}$) are eigenvectors of $\tau_{2}\sigma_{2}$, and
$|\psi_{i}\rangle=|\psi_{i}(y)\rangle$ is a spatially-varying scalar function
satisfying the $i$th differential equation of the set of four differential
equations contained in the Hamiltonian. This yields $|\mu_{i}\rangle$ which
take the forms

$\displaystyle|\mu_{1}\rangle$
$\displaystyle=\frac{1}{\sqrt{2}}(1,0,0,1)^{\top}|\psi_{1}\rangle$
$\displaystyle|\mu_{2}\rangle$
$\displaystyle=\frac{1}{\sqrt{2}}(0,-1,1,0)^{\top}|\psi_{2}\rangle$ (41)
$\displaystyle|\mu_{3}\rangle$
$\displaystyle=\frac{1}{\sqrt{2}}(-1,0,0,1)^{\top}|\psi_{3}\rangle$
$\displaystyle|\mu_{4}\rangle$
$\displaystyle=\frac{1}{\sqrt{2}}(0,1,1,0)^{\top}|\psi_{4}\rangle$

These states serve as zero-energy eigenstates of the effective low-energy
Hamiltonian if the following differential equations are satisfied by the
$\\{\psi_{i}(y)\\}$ (omitting $\frac{1}{4}\Delta_{0}^{2}$ as higher order
terms):

$\displaystyle\left(-\partial_{y}^{2}+\frac{1}{R^{2}}y^{2}-\frac{1}{R}-i\partial_{y}\Delta_{0}\right)|\psi_{1}\rangle$
$\displaystyle=0$ (42)
$\displaystyle\left(-\partial_{y}^{2}+\frac{1}{R^{2}}y^{2}+\frac{1}{R}-i\partial_{y}\Delta_{0}\right)|\psi_{2}\rangle$
$\displaystyle=0$ (43)
$\displaystyle\left(-\partial_{y}^{2}+\frac{1}{R^{2}}y^{2}-\frac{1}{R}+i\partial_{y}\Delta_{0}\right)|\psi_{3}\rangle$
$\displaystyle=0$ (44)
$\displaystyle\left(-\partial_{y}^{2}+\frac{1}{R^{2}}y^{2}+\frac{1}{R}+i\partial_{y}\Delta_{0}\right)|\psi_{4}\rangle$
$\displaystyle=0.$ (45)

Note that by setting $\Delta_{0}=0$ we recover the differential equation in
the two-band Hopf insulator which only has a solution for $-\frac{1}{R}$, so
that we can take this solution, treat $-i\partial_{y}\Delta_{0}$ as a
perturbation term, apply first order perturbation theory (for degenerate
energy levels), and obtain the spatial distribution of the zero-energy modes,
$|\psi_{1}(y)\rangle$.

The zero-energy eigenstates at the inner boundary then take the form,

$\displaystyle|0\rangle_{R}$
$\displaystyle=|\psi_{{\color[rgb]{0,0,0}3}}(y)\rangle\begin{pmatrix}1\\\ 0\\\
0\\\ -1\\\ \end{pmatrix}$ (46) $\displaystyle|0\rangle_{I}$
$\displaystyle=|\psi_{{\color[rgb]{0,0,0}1}}(y)\rangle\begin{pmatrix}i\\\ 0\\\
0\\\ i\\\ \end{pmatrix}.$ (47)

The above analysis is carried out for $x>0$ i.e. the inner boundary. When
instead considering the other boundary corresponding to $x<0$, we reach the
same differential equations (S62) and (S64), this time with eigenstates
$|\mu_{2}\rangle$ and $|\mu_{4}\rangle$ in Eqn (S61). Hence, we similarly find
two additional zero-energy states of the form,

$\displaystyle|\bar{0}\rangle_{R}$
$\displaystyle=|\psi_{{\color[rgb]{0,0,0}1}}(y)\rangle\begin{pmatrix}1\\\ 0\\\
0\\\ -1\\\ \end{pmatrix}$ (48) $\displaystyle|\bar{0}\rangle_{I}$
$\displaystyle=|\psi_{{\color[rgb]{0,0,0}3}}(y)\rangle\begin{pmatrix}i\\\ 0\\\
0\\\ i\\\ \end{pmatrix},$ (49)

where $|\psi_{2}(y)\rangle$ is a corresponding spatial distribution for states
at the outer boundary.

These results indicate that $\mathcal{C}^{\prime}$ symmetry alone yields two
Majorana zero-modes at the defect and two at the boundary, but this is
insufficient to realize the cross zero-modes. $\mathcal{C}^{\prime}$ symmetry
in combination with $\mathcal{C}$ symmetry of the Hopf insulator Bloch
Hamiltonian is also required: the cross states are the most general complex
states permitted by $\mathcal{C}$ and $\mathcal{C}^{\prime}$ symmetry.
Specifically, certain linear combinations of these four cross states satisfy
constraints imposed by $\mathcal{C}$, that the four zero-modes may be
considered as two pairs of MZMs, with one MZM in each pair localized at the
defect and the other localized at the boundary, but may also instead be
considered as two other pairings of MZMs as required by $\mathcal{C}^{\prime}$
symmetry, with one pair localized at the boundary and one pair localized at
the defect. This difference stems from the relationship between $\mathcal{C}$
and $\mathcal{C}^{\prime}$, that $\mathcal{C}^{\prime}$ acts like
$\mathcal{C}$ upon exchange of position and momentum variables. This point
was, to our knowledge, first observed in work by Liu  _et al._ Liu  _et al._
(2017b): under $\mathcal{C}$ transformation, the position variables are
invariant, but all of the momenta variables pick up a sign. However, under the
$\mathcal{C}^{\prime}$ transformation, the momenta do not change, but all of
the position variables pick up a sign. Therefore, for all intents and
purposes, we can replace $\mathcal{C}^{\prime}$ with $\mathcal{C}$ as long as
we reverse the roles of the position and momenta variables.

We may also understand the non-trivial topological state from the perspective
of the homotopy group
$\pi_{3}\left(\mathrm{Sp}(2N)/\mathrm{U}(N)\right)=\mathbb{Z}_{2}$ being
relevant both to mappings from the Brillouin zone to the space of flat-band
Hamiltonians, but also to mappings from the Brillouin zone to the space of
ground state spin expectation value yielding the topological momentum-space,
ground-state spin textures discussed here. This could correspond to a
$\mathbb{Z}_{2}\times\mathbb{Z}_{2}$ topological classification for the system
overall. As a result, even if the topological invariant associated with
mappings from the Brillouin zone to the space of flat-band Hamiltonians is
even and trivial, the system may still harbor topologically-protected zero-
modes due to the co-existing topological skyrmion phase.

This Hamiltonian realizing the cross zero-modes may therefore be interpreted
as a generalized BdG Hamiltonian for a centrosymmetric superconductor, where
$\mathcal{C}$ and $\mathcal{C}^{\prime}$ do not transform the normal state
Bloch Hamiltonian in the same way because it is not even in momentum. For such
cases, the Hilbert space is enlarged with both a particle-hole DOF, and an
additional generalized particle-hole DOF for $\mathcal{C}^{\prime}$. We can
therefore also interpret the full Hamiltonian as a spinless superconductor,
written in terms of the Bloch Hamiltonian as
$\mathcal{H}=\sum_{\boldsymbol{k}}\Psi^{\dagger}_{\boldsymbol{k}}\mathcal{H}(\boldsymbol{k})\Psi_{\boldsymbol{k}}/2$,
where
$\Psi_{\boldsymbol{k}}=\left(c_{\boldsymbol{k}},c^{\dagger}_{-\boldsymbol{k}},c^{\dagger}_{\boldsymbol{k}},c_{-\boldsymbol{k}}\right)^{t}$.
Here, $c_{\boldsymbol{k}}$ is the annihilation operator of electrons with
momentum $\boldsymbol{k}=\left(k_{x},k_{y},k_{z}\right)$.

It is interesting that the cross zero-modes are strongly localized only at one
boundary or the other: this is due to only two linearly-independent cross
states being possible at a given point in real-space. The cross-states
therefore must additionally be orthogonal as a result of their spatial
distribution. Therefore two cross states are strongly localized at the inner
boundary, and the remaining two are strongly localized at the outer boundary
to form a four-fold degenerate ground state manifold. Importantly, this
indicates the cross zero-modes should be more robust against finite-size
effects than unpaired Majorana zero-modes. This is supported by numerics,
which show well-defined cross zero-modes at the outer boundary in smaller
systems (see 13 (f)).

In summary, in the topologically non-trivial regime of the Hamiltonian
presented here realizing the 3D topological skyrmion phase, corresponding to a
skyrmion forming in the ground-state spin expectation value texture at half-
filling, there are two cross states localized on the outer boundary, and two
other cross states localized at the defect. These cross states are a
manifestation of a bulk-boundary correspondence of topological skyrmion phases
due to topology of defects trapped by the momentum-space spin textures.

## IV Discussion & Conclusion

In this work, we present toy models for recently-identified topological phases
of matter in two dimensions, characterized by the topological charge of a
momentum-space spin texture, a spin topological invariant rather than the
topological charge of a momentum-space texture of the projector onto occupied
states, a projector topological invariant, as done in the ten-fold way
classification scheme Schnyder  _et al._ (2008b). We realize these topological
phases in systems which can therefore be more fully characterized in terms of
both a projector topological invariant and a spin topological invariant.
Extending our results to introduce three-dimensional topological skyrmion
phases corresponding to trivial projector topological invariant and non-
trivial spin topological invariant, we furthermore find evidence of a bulk-
boundary correspondence of topological skyrmion phases. As bulk-boundary
correspondence is generally expected only for non-trivial projector
topological invariant, these results serve as an important signature and
consequence of topological skyrmion phases.

The bulk-boundary correspondence considered here, a defect bulk-boundary
correspondence, is realized by opening boundary conditions of the 3D system in
two directions $x$ and $y$, and replacing the momentum component in a third
$z$ direction by an angle $\theta(x,y)$, which can characterize a zero-
dimensional defect in the $xy$-plane. For these open boundary conditions, we
characterize the system’s topology in terms of a bulk projector invariant, the
Hopf invariant, and a bulk spin invariant, the 3D skyrmion number. Under these
open boundary conditions, the topological skyrmion phases corresponding to
trivial projector invariant and non-trivial spin invariant possess
topologically-protected “cross” zero-modes, which generalize unpaired Majorana
zero-modes: the zero-mode wavefunctions in lattice-based systems consist of a
set of complex numbers that form a cross-like structure when each of these
complex numbers is plotted in the complex plane, rather than forming a
straight line corresponding to a purely real state up to some phase degree of
freedom as in the case of an unpaired Majorana zero-mode.

In order to demonstrate this defect bulk-boundary correspondence of the
topological skyrmion phases, we first present a method for constructing Bloch
Hamiltonian models for topological skyrmion phases with four bands and
characterize myriad examples. We furthermore show BdG Hamiltonians for
centrosymmetric superconductors naturally possess the $\mathcal{C}^{\prime}$
symmetry in addition to particle-hole symmetry $\mathcal{C}$ when the normal
state Hamiltonian is even in momentum. More generally, the topological
skyrmion phases are described by generalized BdG Hamiltonians corresponding to
artificially enlarging the Hilbert space to possess a generalized particle-
hole DOF in addition to a particle-hole DOF.

We introduce an efficient method for computing skyrmion numbers based on past
work by Fukui  _et al._ Fukui  _et al._ (2005). At the four-band level, we
find the skyrmion number computed using the ground state spin expectation
value at half-filling is no longer equal to the total Chern number. Instead,
we find that the total Chern number $\mathcal{C}$ and the skyrmion number
$\mathcal{Q}$ are related by $\mathcal{C}=-2\mathcal{Q}$ analytically for the
atomic model—and numerically for all other models considered—of the chiral
topological skyrmion phase characterized by Bloch Hamiltonians with four
bands, when each of $\mathcal{C}$ and $\mathcal{Q}$ are topologically stable.
It is also possible for $\mathcal{C}$ to be topologically stable and
$\mathcal{Q}$ to be topologically unstable, due to the magnitude of the ground
state spin being zero somewhere in the Brillouin zone. We do not find any
instances where the reverse case occurs, in which $\mathcal{Q}$ is
topologically stable and $\mathcal{C}$ is topologically unstable, and will
explore this further in future work.

We do not explore bulk-boundary correspondence of two-dimensional topological
skyrmion phases in this work. However, since completion of this work, bulk-
boundary correspondence for the 2D topological skyrmion phases has also been
further characterized in terms of the observable-enriched partial trace and
observable-enriched entanglement spectrum introduced by Winter  _et al._ Cook
(2023); Winter  _et al._ . That is, bulk-boundary correspondence associated
with non-trivial skyrmion number $\mathcal{Q}$, when boundary conditions are
opened in one direction, corresponds to $\mathcal{Q}$ chiral modes in the
observable-enriched entanglement spectrum bulk gap.

We use knowledge from constructing various toy models for the chiral
topological skyrmion phase in two dimensions to introduce the three-
dimensional topological skyrmion phase, constructing it from the Hopf
insulator Hamiltonian. We characterize the topology of these models using both
a projector topological invariant and a spin topological invariant. We then
show that the three-dimensional spin skyrmion in momentum space traps a
defect, which can yield a bulk-boundary correspondence for the three-
dimensional skyrmion phase: taking a momentum component to be an angle
$\theta$ characterizing a real-space defect in a two-dimensional system with
open boundary conditions in the other two directions, we find two zero-modes
are localized at the defect and two are localized at the boundary when the
system is in the three-dimensional topological skyrmion phase. These zero-
modes are realized when the projector topological invariant, the Hopf
invariant, is trivial and the spin topological invariant, the 3D skyrmion
number, is non-trivial. Therefore, while bulk-boundary correspondence has
previously been associated strictly with non-trivial projector topological
invariants, we show bulk-boundary correspondence is also associated with non-
trival spin topological invariant. Such a zero-mode exhibits a cross-like
structure when each of the complex numbers characterizing the zero-mode
wavefunction is plotted in the complex plane. Such a cross zero-mode is
protected by a combination of particle-hole and generalized particle-hole
symmetries and is robust against disorder. They generalize unpaired Majorana
zero-modes, and may be thought of as an unpaired Majorana zero-mode and its
$\mathcal{C}^{\prime}$ partner combined in a single, localized state, yet
isolated from one another by a $\pi/2$ phase difference in the complex plane.

In conclusion, these results indicate that the topological skyrmion phases
exhibit a bulk-boundary correspondence by trapping defects that can yield
topologically-protected, gapless boundary states. This serves as confirmation
of the topological skyrmion phases and their implications, including their
type-II topological phase transition contradicting the flat-band limit
assumption. As the flat-band limit assumption is used in the construction of
the ten-fold way topological classification scheme of Schnyder et al. Schnyder
_et al._ (2008b) and Ryu et al. Ryu  _et al._ (2010b), the type-II topological
phase transition of the topological skyrmion phases Cook (2023) indicates
topological skyrmion phases are a new topological class outside the
classification of Schnyder et al. Schnyder  _et al._ (2008b) and Ryu et al.
Ryu  _et al._ (2010b), rather than additional topological characterization of
the current topological classification of free fermions. In addition, these
gapless boundary states due to defect bulk-boundary correspondence of the
topological skyrmion phases may interact with the boundary states of other
topological phases. Such interplay between myriad topological phases will be
explored in future work, especially given the potential for the topological
skyrmion phase and the Chern insulator phase to co-exist. Finally, the
properties and consequences of the cross zero-modes of the three-dimensional
topological skyrmion phase will also be explored further in future work, given
the significance of unpaired Majorana zero-modes for topological quantum
computation and the fundamental generalization of unpaired Majorana zero-modes
presented by the cross zero-modes discovered here.

Acknowledgements \- The authors thank B. Braunecker and B. Doucot for helpful
discussions and comments on the manuscript.

## References

  * Chiu  _et al._ (2016) C.-K. Chiu, J. C. Y. Teo, A. P. Schnyder, and S. Ryu, Rev. Mod. Phys. 88, 035005 (2016).
  * Lahtinen and Pachos (2017) V. T. Lahtinen and J. K. Pachos, SciPost Phys. 3, 021 (2017).
  * Hasan and Kane (2010) M. Z. Hasan and C. L. Kane, Rev. Mod. Phys. 82, 3045 (2010).
  * Bernevig and Hughes (2013) B. A. Bernevig and T. L. Hughes,  _Topological Insulators and Topological Superconductors_ (Princeton University Press, 2013).
  * Qi and Zhang (2011) X.-L. Qi and S.-C. Zhang, Rev. Mod. Phys. 83, 1057 (2011).
  * He  _et al._ (2021) Q. L. He, T. L. Huges, N. P. Armitage, Y. Tokura, and K. L. Wang, Nature Materials (2021).
  * Nayak  _et al._ (2008) C. Nayak, S. H. Simon, A. Stern, M. Freedman, and S. Das Sarma, Rev. Mod. Phys. 80, 1083 (2008).
  * Kitaev (2001) A. Y. Kitaev, Physics-Uspekhi 44, 131 (2001).
  * Bernevig  _et al._ (2006) B. A. Bernevig, T. L. Hughes, and S.-C. Zhang, Science 314, 1757 (2006), http://science.sciencemag.org/content/314/5806/1757.full.pdf .
  * Hasan and Moore (2011) M. Z. Hasan and J. E. Moore, Ann. Rev. Condensed Matter Physics 2, 55 (2011).
  * Schnyder  _et al._ (2008a) A. P. Schnyder, S. Ryu, A. Furusaki, and A. W. W. Ludwig, Phys. Rev. B 78, 195125 (2008a).
  * Ryu  _et al._ (2010a) S. Ryu, A. P. Schnyder, A. Furusaki, and A. W. W. Ludwig, New Journal of Physics 12, 065010 (2010a).
  * Tiwari  _et al._ (2018) A. Tiwari, X. Chen, K. Shiozaki, and S. Ryu, Phys. Rev. B , 245133 (2018).
  * Fu (2011) L. Fu, Phys. Rev. Lett. 106, 106802 (2011).
  * Hsieh  _et al._ (2012) T. H. Hsieh, H. Lin, J. Liu, W. Duan, A. Bansil, and L. Fu, Nature Communications 3, 982 EP (2012).
  * Huang  _et al._ (2017) S.-J. Huang, H. Song, Y.-P. Huang, and M. Hermele, Phys. Rev. B 96, 205106 (2017).
  * Song  _et al._ (2017) H. Song, S.-J. Huang, L. Fu, and M. Hermele, Phys. Rev. X 7, 011020 (2017).
  * Cornfeld and Carmeli (2021) E. Cornfeld and S. Carmeli, Phys. Rev. Research 3, 013052 (2021).
  * Chen  _et al._ (2013a) X. Chen, Z.-C. Gu, Z.-X. Liu, and X.-G. Wen, Phys. Rev. B 87, 155114 (2013a).
  * Gong  _et al._ (2016) Z.-X. Gong, M. F. Maghrebi, A. Hu, M. L. Wall, M. Foss-Feig, and A. V. Gorshkov, Phys. Rev. B 93, 041102 (2016).
  * Khanikaev  _et al._ (2012) A. B. Khanikaev, S. Hossein Mousavi, W.-K. Tse, M. Kargarian, A. H. MacDonald, and G. Shvets, Nature Materials 12, 233 EP (2012).
  * Lan and Wen (2014) T. Lan and X.-G. Wen, Phys. Rev. B 90, 115119 (2014).
  * Vijay  _et al._ (2015) S. Vijay, J. Haah, and L. Fu, Phys. Rev. B 92, 235136 (2015).
  * Dai  _et al._ (2020) Y. Dai, Z. Zhou, A. Ghosh, A. Kubo, C.-B. Huang, and H. Petek, Nature 588, 616 (2020).
  * Tokuda  _et al._ (2021) S. Tokuda, S. Souma, K. Segawa, T. Takshashi, Y. Ando, T. Nakanishi, and T. Sato, Commun Phys 4, 170 (2021).
  * Lindner  _et al._ (2011) N. H. Lindner, G. Refael, and V. Galitski, Nature Physics 7, 490 (2011).
  * Cayssol  _et al._ (2013) J. Cayssol, B. Dóra, F. Simon, and R. Moessner, Physica Status Solidi (RRL) – Rapid Research Letters 7, 101 (2013), https://onlinelibrary.wiley.com/doi/pdf/10.1002/pssr.201206451 .
  * Khemani  _et al._ (2016) V. Khemani, A. Lazarides, R. Moessner, and S. L. Sondhi, Phys. Rev. Lett. 116, 250401 (2016).
  * Crowley  _et al._ (2019) P. J. D. Crowley, I. Martin, and A. Chandran, Phys. Rev. B 99, 064306 (2019).
  * Wintersperger  _et al._ (2020) K. Wintersperger, C. Braun, F. N. Ünal, A. Eckardt, M. Di Liberto, N. Goldman, I. Bloch, and M. Aidelsburger, Nature Physics 16, 1058–1063 (2020).
  * Shen  _et al._ (2018) H. Shen, B. Zhen, and L. Fu, Phys. Rev. Lett. 120, 146402 (2018).
  * Gong  _et al._ (2018) Z. Gong, Y. Ashida, K. Kawabata, K. Takasan, S. Higashikawa, and M. Ueda, Phys. Rev. X 8, 031079 (2018).
  * Zheng and Hofstetter (2018) J.-H. Zheng and W. Hofstetter, Phys. Rev. B 97, 195434 (2018).
  * Kawabata  _et al._ (2019) K. Kawabata, S. Higashikawa, Z. Gong, Y. Ashida, and M. Ueda, Nature Communications 10, 297 (2019).
  * Cook (2023) A. M. Cook, Journal of Physics: Condensed Matter 35, 184001 (2023).
  * Rößler  _et al._ (2006) U. K. Rößler, A. N. Bogdanov, and C. Pfleiderer, Nature 442, 797 (2006).
  * Mühlbauer  _et al._ (2009) S. Mühlbauer, B. Binz, F. Jonietz, C. Pfleiderer, A. Rosch, A. Neubauer, R. Georgii, and P. Böni, Science 323, 915 (2009).
  * Romming  _et al._ (2013) N. Romming, C. Hanneken, M. Menzel, J. E. Bickel, B. Wolter, K. von Bergmann, A. Kubetzka, and R. Wiesendanger, Science 341, 636 (2013).
  * Fert  _et al._ (2017) A. Fert, N. Reyren, and V. Cros, Nature Reviews Materials 2, 17031 (2017).
  * Teo and Kane (2010) J. C. Y. Teo and C. L. Kane, Phys. Rev. B 82, 115120 (2010).
  * Moore  _et al._ (2008) J. E. Moore, Y. Ran, and X.-G. Wen, Phys. Rev. Lett. 101, 186805 (2008).
  * Deng  _et al._ (2013) D.-L. Deng, S.-T. Wang, C. Shen, and L.-M. Duan, Phys. Rev. B 88, 201105 (2013).
  * Liu  _et al._ (2017a) C. Liu, F. Vafa, and C. Xu, Phys. Rev. B 95, 161116 (2017a).
  * Yan  _et al._ (2017) Z. Yan, R. Bi, and Z. Wang, Phys. Rev. Lett. 118, 147003 (2017).
  * Liu  _et al._ (2017b) C. Liu, F. Vafa, and C. Xu, Phys. Rev. B 95, 161116 (2017b).
  * Fukui  _et al._ (2005) T. Fukui, Y. Hatsugai, and H. Suzuki, Journal of the Physical Society of Japan 74, 1674 (2005), https://doi.org/10.1143/JPSJ.74.1674 .
  * Chen  _et al._ (2013b) Y. Chen, D. L. Bergman, and A. A. Burkov, Phys. Rev. B 88, 125110 (2013b).
  * Ueno  _et al._ (2013) Y. Ueno, A. Yamakage, Y. Tanaka, and M. Sato, Phys. Rev. Lett. 111, 087002 (2013).
  * Sticlet  _et al._ (2012) D. Sticlet, F. Piéchon, J.-N. Fuchs, P. Kalugin, and P. Simon, Phys. Rev. B 85, 165456 (2012).
  * Ryu  _et al._ (2010b) S. Ryu, A. P. Schnyder, A. Furusaki, and A. W. W. Ludwig, New Journal of Physics 12, 065010 (2010b).
  * Yan and Felser (2017) B. Yan and C. Felser, Annual Review of Condensed Matter Physics 8, 337 (2017), https://doi.org/10.1146/annurev-conmatphys-031016-025458 .
  * Alicea (2012) J. Alicea, Reports on Progress in Physics 75, 076501 (2012).
  * Read and Green (2000) N. Read and D. Green, Phys. Rev. B 61, 10267 (2000).
  * Schnyder  _et al._ (2008b) A. P. Schnyder, S. Ryu, A. Furusaki, and A. W. W. Ludwig, Phys. Rev. B 78, 195125 (2008b).
  * (55) J. H. Winter, B. Braunecker, and A. M. Cook, in progress. .

Supplemental material for “Defect bulk-boundary correspondence of topological
skyrmion phases of matter”
Shu-Wei Liu,1 Li-kun Shi,1 and Ashley M. Cook1,2,∗

1Max Planck Institute for Chemical Physics of Solids, Nöthnitzer Strasse 40,
01187 Dresden, Germany
2Max Planck Institute for the Physics of Complex Systems, Nöthnitzer Strasse
38, 01187 Dresden, Germany
∗Electronic address<EMAIL_ADDRESS>
(Dated: )

## S1 Low-energy, analytical theory of the 3D topological skyrmion phase

Here we derive a low-energy theory for the four-band, three-dimensional
skyrmion phase with Hamiltonian defined in Eq. (34), for a semi-infinite
geometry with the sample occupying the $x>0$ region. Translational symmetry in
the $\hat{x}$-direction is therefore broken while momentum component $k_{y}$
remains as a good quantum number. Thus, a real-space coordinate is used in the
$\hat{x}$-direction to label lattice sites. The wave functions and energy
eigenvalues can be obtained by solving the following equation

$\begin{pmatrix}H_{0}&H_{1}&H_{2}&0&0&0&\cdots\\\
H_{1}^{\dagger}&H_{0}&H_{1}&H_{2}&0&0&\cdots\\\
H_{2}^{\dagger}&H_{1}^{\dagger}&H_{0}&H_{1}&H_{2}&0&\cdots\\\
0&H_{2}^{\dagger}&H_{1}^{\dagger}&H_{0}&H_{1}&H_{2}&\cdots\\\
\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\ddots\\\
\end{pmatrix}\begin{pmatrix}\psi_{1}\\\ \psi_{2}\\\ \psi_{3}\\\ \psi_{4}\\\
\vdots\end{pmatrix}=E\begin{pmatrix}\psi_{1}\\\ \psi_{2}\\\ \psi_{3}\\\
\psi_{4}\\\ \vdots\end{pmatrix}$ (S1)

where $\psi_{j}$ is the $j^{th}$ entry of an eigenstate for this equation
corresponding to the $j$th lattice site in the $\hat{x}$-direction. Each
$\psi_{j}$ has four components corresponding to the generalized particle-hole
and spin degrees of freedom at each site in real-space, and the non-zero terms
of the Hamiltonian matrix representation on the left-hand side of Eq. S1 take
the following forms: the term on the diagonal is

$H_{0}(k_{y},\lambda,\Delta_{0})=\begin{pmatrix}h_{0}&0\\\ 0&-h_{0}^{*}\\\
\end{pmatrix}+\Delta_{s},$ (S2)

and terms off the diagonal are

$H_{1}(k_{y},\lambda,\Delta_{0})=\begin{pmatrix}h_{1}&0\\\ 0&-h_{1}^{*}\\\
\end{pmatrix}+\begin{pmatrix}0&0&-\frac{1}{2}\Delta_{0}&0\\\
0&0&0&-\frac{1}{2}\Delta_{0}\\\ \frac{1}{2}\Delta_{0}&0&0&0\\\
0&\frac{1}{2}\Delta_{0}&0&0\end{pmatrix},$ (S3)

and

$\displaystyle H_{2}$
$\displaystyle=\begin{pmatrix}-\frac{1}{2}&-\frac{1}{2}&0&0\\\
\frac{1}{2}&\frac{1}{2}&0&0\\\ 0&0&\frac{1}{2}&\frac{1}{2}\\\
0&0&-\frac{1}{2}&-\frac{1}{2}\end{pmatrix},$ (S4)

where $\lambda$ is a parameter substituted for the $\hat{z}$-component of the
momentum, $k_{z}$, characterizing a defect. Here, $\lambda=n\theta(i,j)$ with
$\theta$ being the relative angle between the defect and site$(i,j)$, $n$ is
an integer taken to be $1$ in this study and $\Delta_{0}$ is the strength of
the spin-triplet pairing.

Here, $h_{0}$ and $h_{1}$ are momentum- and $\lambda$-dependent terms taking
the forms

$\displaystyle h_{0}(k_{y},\lambda)$
$\displaystyle=\begin{pmatrix}A(k_{y},\lambda)&B(k_{y},\lambda)\\\
B(k_{y},\lambda)&-A(k_{y},\lambda)\end{pmatrix}$ (S5) $\displaystyle
h_{1}(k_{y},\lambda)$
$\displaystyle=\begin{pmatrix}\frac{3}{2}-\cos{k_{y}}-\cos{\lambda}&-i\sin{\lambda}+\sin{k_{y}}-(\cos{k_{y}}+\cos{\lambda}-\frac{3}{2})\\\
-i\sin{\lambda}+\sin{k_{y}}+(\cos{k_{y}}+\cos{\lambda}-\frac{3}{2})&-\frac{3}{2}+\cos{k_{y}}+\cos{\lambda}\end{pmatrix}$

where here

$\displaystyle A(k_{y},\lambda)$
$\displaystyle=-\frac{13}{4}-\cos{2k_{y}}+3(\cos{k_{y}}+\cos{\lambda})-2\cos{k_{y}}\cos{\lambda}$
(S6) $\displaystyle B(k_{y},\lambda)$
$\displaystyle=2i\sin{k_{y}}\sin{\lambda}+(\sin{2k_{y}}+2\sin{k_{y}}\cos{\lambda}-3\sin{k_{y}}).$

$\Delta_{s}$ is also momentum-dependent, taking the form

$\displaystyle\Delta_{s}(k_{y},\Delta_{0})$
$\displaystyle=\begin{pmatrix}0&0&-\Delta_{0}\sin{k_{y}}&0\\\
0&0&0&\Delta_{0}\sin{k_{y}}\\\ -\Delta_{0}\sin{k_{y}}&0&0&0\\\
0&\Delta_{0}\sin{k_{y}}&0&0\end{pmatrix}.$ (S7)

To find zero-energy solutions in analogy to Yan  _et al._ Yan  _et al._
(2017), we first set $k_{y}=0$, $\lambda=0$ and $\Delta_{0}=0$. The non-zero
entries of the Hamiltonian in Eq. S1 then take the following forms:

$\displaystyle H_{0}(0,0,0)$
$\displaystyle=\begin{pmatrix}-\frac{1}{4}&0&0&0\\\ 0&\frac{1}{4}&0&0\\\
0&0&\frac{1}{4}&0\\\ 0&0&0&-\frac{1}{4}\end{pmatrix}$ (S8) $\displaystyle
H_{1}(0,0,0)$ $\displaystyle=\begin{pmatrix}-\frac{1}{2}&-\frac{1}{2}&0&0\\\
\frac{1}{2}&\frac{1}{2}&0&0\\\ 0&0&\frac{1}{2}&\frac{1}{2}\\\
0&0&-\frac{1}{2}&-\frac{1}{2}\end{pmatrix}=H_{2}$

We then consider eigenstates of

$\tau_{0}\otimes\sigma_{1}=\begin{pmatrix}0&1&0&0\\\ 1&0&0&0\\\ 0&0&0&1\\\
0&0&1&0\end{pmatrix},$ (S9)

which are

$\displaystyle|\nu_{1}\rangle$
$\displaystyle=\frac{1}{\sqrt{2}}(0,0,-1,1)^{\top}$ (S10)
$\displaystyle|\nu_{2}\rangle$
$\displaystyle=\frac{1}{\sqrt{2}}(-1,1,0,0)^{\top}$
$\displaystyle|\nu_{3}\rangle$
$\displaystyle=\frac{1}{\sqrt{2}}(0,0,1,1)^{\top}$
$\displaystyle|\nu_{4}\rangle$
$\displaystyle=\frac{1}{\sqrt{2}}(1,1,0,0)^{\top}$

where $|\nu_{1}\rangle$ and $|\nu_{2}\rangle$ have eigenvalue
$\lambda_{-}=-1$, and $|\nu_{3}\rangle$ and $|\nu_{4}\rangle$ have eigenvalue
$\lambda_{+}=+1$. Operating the matrices $H_{i}$’s on the $|\nu_{i}\rangle$’s,
they produce:

$\displaystyle H_{0}\begin{pmatrix}|\nu_{1}\rangle\\\ |\nu_{2}\rangle\\\
|\nu_{3}\rangle\\\ |\nu_{4}\rangle\\\ \end{pmatrix}$
$\displaystyle=\frac{1}{4}\begin{pmatrix}-|\nu_{3}\rangle\\\
|\nu_{4}\rangle\\\ -|\nu_{1}\rangle\\\ |\nu_{2}\rangle\\\ \end{pmatrix}$ (S11)
$\displaystyle H_{1}\begin{pmatrix}|\nu_{1}\rangle\\\ |\nu_{2}\rangle\\\
|\nu_{3}\rangle\\\ |\nu_{4}\rangle\\\ \end{pmatrix}$
$\displaystyle=\begin{pmatrix}0\\\ 0\\\ -|\nu_{1}\rangle\\\ |\nu_{2}\rangle\\\
\end{pmatrix}$ $\displaystyle H_{1}^{\dagger}\begin{pmatrix}|\nu_{1}\rangle\\\
|\nu_{2}\rangle\\\ |\nu_{3}\rangle\\\ |\nu_{4}\rangle\\\ \end{pmatrix}$
$\displaystyle=\begin{pmatrix}-|\nu_{3}\rangle\\\ |\nu_{4}\rangle\\\ 0\\\ 0\\\
\end{pmatrix}$

The zero-energy wavefunctions take the form

$|\Psi_{i}\rangle=\sum_{j}a_{i,j}|j\rangle\otimes|\nu_{i}\rangle$ (S12)

where $a_{i,j}$ is a function of $j$ that normalizes $|\Psi_{i}\rangle$ and
$|j\rangle$ indicates localization on the $j^{th}$ site:

$|j\rangle=(0,0,\cdots,0,1,0,\cdots)^{\top}$ (S13)

i.e. having $1$ in the $j^{th}$ element and $0$ everywhere else. We can
explore different possibilities of $|\Psi_{i}\rangle$ by substituting this
expression in Eq.(S1) and consider the $j^{th}$ element of the equation. In
the case of $i=1$:

$H_{2}^{\dagger}a_{1,j-2}|\nu_{1}\rangle+H_{1}^{\dagger}a_{1,j-1}|\nu_{1}\rangle+H_{0}a_{1,j}|\nu_{1}\rangle+H_{1}a_{1,j+1}|\nu_{1}\rangle+H_{2}a_{1,j+2}|\nu_{1}\rangle=E=0$
(S14)

which simplifies to

$(-a_{1,j-2}-a_{1,j-1}-\frac{1}{4}a_{1,j}+0+0)|\nu_{3}\rangle=0$ (S15)

and therefore obtaining the following recurrence relation in $a_{1,j}$:

$a_{1,j}=-4(a_{1,j-2}+a_{1,j-1})$ (S16)

which solves to give the general formula:

$a_{1,j}=(-2)^{j}C_{1}+(-2)^{j}jC_{2}$ (S17)

but there is no normalizable solution for $C_{1}$ or $C_{2}$.
In the case of $i=2$:

$H_{2}^{\dagger}a_{2,j-2}|\nu_{2}\rangle+H_{1}^{\dagger}a_{2,j-1}|\nu_{2}\rangle+H_{0}a_{2,j}|\nu_{2}\rangle+H_{1}a_{2,j+1}|\nu_{2}\rangle+H_{2}a_{2,j+2}|\nu_{2}\rangle=E=0$
(S18)

which simplifies to

$(a_{2,j-2}+a_{2,j-1}+\frac{1}{4}a_{1,j}+0+0)|\nu_{4}\rangle=0$ (S19)

and therefore obtaining the following recurrence relation in $a_{2,j}$:

$a_{2,j}=-4(a_{2,j-2}+a_{2,j-1}).$ (S20)

which is identical to the $i=1$ case and is unnormalizable.

In the case of $i=3$:

$H_{2}^{\dagger}a_{3,j-2}|\nu_{3}\rangle+H_{1}^{\dagger}a_{3,j-1}|\nu_{3}\rangle+H_{0}a_{3,j}|\nu_{3}\rangle+H_{1}a_{3,j+1}|\nu_{3}\rangle+H_{2}a_{3,j+2}|\nu_{3}\rangle=E=0$
(S21)

which simplifies to

$(0+0-\frac{1}{4}a_{3,j}-a_{2,j+1}-a_{2,j+2})|\nu_{1}\rangle=0$ (S22)

and therefore obtaining the following recurrence relation in $a_{3,j}$:

$a_{3,j}=-4(a_{3,j+1}+a_{3,j+2}).$ (S23)

which solves to give the general formula:

$a_{3,j}=(-\frac{1}{2})^{j}C_{1}+(-\frac{1}{2})^{j}jC_{2}$ (S24)

which now has normalizable solution for $C_{1}$ and $C_{2}$. Note that each of
the two terms in $a_{3,j}$ can form an individual solution, so we denote the
two choices of $a_{3,j}$ as $\alpha_{j}$ and $\beta_{j}$, and thus obtain the
following expression:

$\displaystyle\alpha_{j}$ $\displaystyle=\sqrt{3}(-\frac{1}{2})^{j}$ (S25)
$\displaystyle\beta_{j}$
$\displaystyle=-\sqrt{\frac{27}{20}}(-\frac{1}{2})^{j-1}(j-1)$

and in summary:

$\displaystyle|\Psi_{\alpha}\rangle$
$\displaystyle=\sum_{j}\alpha_{j}|j\rangle\otimes|\nu_{3}\rangle$ (S26)
$\displaystyle|\Psi_{\beta}\rangle$
$\displaystyle=\sum_{j}\beta_{j}|j\rangle\otimes|\nu_{3}\rangle$

such that
$\langle\Psi_{\alpha}|\Psi_{\alpha}\rangle=\langle\Psi_{\beta}|\Psi_{\beta}\rangle=1.$

In the case of $i=4$:

$H_{2}^{\dagger}a_{4,j-2}|\nu_{4}\rangle+H_{1}^{\dagger}a_{4,j-1}|\nu_{4}\rangle+H_{0}a_{4,j}|\nu_{2}\rangle+H_{1}a_{4,j+1}|\nu_{4}\rangle+H_{2}a_{4,j+2}|\nu_{4}\rangle=E=0$
(S27)

which simplifies to

$(0+0+\frac{1}{4}a_{4,j}+a_{2,j+1}+a_{2,j+2})|\nu_{2}\rangle=0$ (S28)

and therefore obtaining the following recurrence relation in $a_{4,j}$:

$a_{4,j}=-4(a_{4,j+1}+a_{4,j+2}).$ (S29)

which is an identical form to the $i=3$ case so the solutions are

$\displaystyle\gamma_{j}$ $\displaystyle=\sqrt{3}(-\frac{1}{2})^{j}$ (S30)
$\displaystyle\delta_{j}$
$\displaystyle=-\sqrt{\frac{27}{20}}(-\frac{1}{2})^{j-1}(j-1)$

and

$\displaystyle|\Psi_{\gamma}\rangle$
$\displaystyle=\sum_{j}\gamma_{j}|j\rangle\otimes|\nu_{4}\rangle$ (S31)
$\displaystyle|\Psi_{\delta}\rangle$
$\displaystyle=\sum_{j}\delta_{j}|j\rangle\otimes|\nu_{4}\rangle$

such that
$\langle\Psi_{\gamma}|\Psi_{\gamma}\rangle=\langle\Psi_{\delta}|\Psi_{\delta}\rangle=1.$
The next step in normalization is to notice $|\Psi_{\alpha}\rangle$ and
$|\Psi_{\beta}\rangle$ are both constructed from $|\nu_{3}\rangle$ so
$\langle\Psi_{\alpha}|\Psi_{\beta}\rangle$ is not automatically 0. We apply
the Gram-Schmidt orthogonalization as the following:

$\displaystyle|\Psi_{\alpha}^{o}\rangle$ $\displaystyle=|\Psi_{\alpha}\rangle$
(S32) $\displaystyle|\Psi_{\beta}^{o}\rangle$
$\displaystyle=N(|\Psi_{\beta}\rangle-|\Psi_{\alpha}\rangle\langle\Psi_{\alpha}|\Psi_{\beta}\rangle)$

so that $\langle\Psi_{\alpha}^{o}|\Psi_{\beta}^{o}\rangle$ is conveniently 0
and the normalization constant $N$ can be worked out from first determining

$\displaystyle\langle\Psi_{\alpha}|\Psi_{\beta}\rangle$
$\displaystyle=\sum_{j}\alpha^{*}_{j}\beta_{j}$ (S33)
$\displaystyle=\sum_{j}-\sqrt{\frac{81}{20}}(-\frac{1}{2})^{2j-1}(j-1)$
$\displaystyle=-\sqrt{\frac{81}{20}}(-\frac{2}{9})$
$\displaystyle=\frac{1}{\sqrt{5}}$

and then by requiring

$\displaystyle\langle\Psi_{\beta}^{o}|\Psi_{\beta}^{o}\rangle$
$\displaystyle=N^{2}(\langle\Psi_{\beta}|-\frac{1}{\sqrt{5}}\langle\Psi_{\alpha}|)(|\Psi_{\beta}\rangle-\frac{1}{\sqrt{5}}|\Psi_{\alpha}\rangle)=1$
(S34) $\displaystyle N^{2}(1-\frac{1}{5}-\frac{1}{5}+\frac{1}{5})$
$\displaystyle=1$ $\displaystyle N$ $\displaystyle=\frac{\sqrt{5}}{2}$

and therefore

$\displaystyle|\Psi_{\alpha}^{o}\rangle$ $\displaystyle=|\Psi_{\alpha}\rangle$
(S35) $\displaystyle|\Psi_{\beta}^{o}\rangle$
$\displaystyle=\frac{\sqrt{5}}{2}|\Psi_{\beta}\rangle-\frac{1}{2}|\Psi_{\alpha}\rangle$

and similarly

$\displaystyle|\Psi_{\gamma}^{o}\rangle$ $\displaystyle=|\Psi_{\gamma}\rangle$
(S36) $\displaystyle|\Psi_{\delta}^{o}\rangle$
$\displaystyle=\frac{\sqrt{5}}{2}|\Psi_{\delta}\rangle-\frac{1}{2}|\Psi_{\gamma}\rangle.$

Next, we consider the Hamiltonian in the neighbourhood of
$k_{y}=0,\lambda=0,\Delta_{0}=0$ by expanding $H$ to first order of these
three variables i.e. $H_{i}(k_{y},\lambda,\Delta_{0})=H_{i}(0,0,0)+\Delta
H_{i}(k_{y},\lambda,\Delta_{0})$ where

$\Delta H=\begin{pmatrix}\Delta H_{0}&\Delta H_{1}&0&0&0&\cdots\\\ \Delta
H_{1}^{\dagger}&\Delta H_{0}&\Delta H_{1}&0&0&\cdots\\\ 0&\Delta
H_{1}^{\dagger}&\Delta H_{0}&\Delta H_{1}&0&\cdots\\\
\vdots&\vdots&\vdots&\vdots&\vdots&\ddots\\\ \end{pmatrix}$ (S37)

with

$\displaystyle\Delta H_{0}$ $\displaystyle=\begin{pmatrix}0&k_{y}&0&0\\\
k_{y}&0&0&0\\\ 0&0&0&-k_{y}\\\ 0&0&-k_{y}&0\end{pmatrix}$ (S38)
$\displaystyle\Delta H_{1}$
$\displaystyle=\begin{pmatrix}0&-i\lambda+k_{y}&-\frac{1}{2}\Delta_{0}&0\\\
-i\lambda+k_{y}&0&0&-\frac{1}{2}\Delta_{0}\\\
\frac{1}{2}\Delta_{0}&0&0&-i\lambda-k_{y}\\\
0&\frac{1}{2}\Delta_{0}&-i\lambda-k_{y}&0\end{pmatrix},$ (S39)

and the low-energy effective Hamiltonian is

$\displaystyle H_{\textnormal{eff}}$
$\displaystyle=\begin{pmatrix}\langle\Psi_{\alpha}^{o}|\Delta
H|\Psi_{\alpha}^{o}\rangle&\langle\Psi_{\alpha}^{o}|\Delta
H|\Psi_{\beta}^{o}\rangle&\langle\Psi_{\alpha}^{o}|\Delta
H|\Psi_{\gamma}^{o}\rangle&\langle\Psi_{\alpha}^{o}|\Delta
H|\Psi_{\delta}^{o}\rangle\\\ \langle\Psi_{\beta}^{o}|\Delta
H|\Psi_{\alpha}^{o}\rangle&\langle\Psi_{\beta}^{o}|\Delta
H|\Psi_{\beta}^{o}\rangle&\langle\Psi_{\beta}^{o}|\Delta
H|\Psi_{\gamma}^{o}\rangle&\langle\Psi_{\beta}^{o}|\Delta
H|\Psi_{\delta}^{o}\rangle\\\ \langle\Psi_{\gamma}^{o}|\Delta
H|\Psi_{\alpha}^{o}\rangle&\langle\Psi_{\gamma}^{o}|\Delta
H|\Psi_{\beta}^{o}\rangle&\langle\Psi_{\gamma}^{o}|\Delta
H|\Psi_{\gamma}^{o}\rangle&\langle\Psi_{\gamma}^{o}|\Delta
H|\Psi_{\delta}^{o}\rangle\\\ \langle\Psi_{\delta}^{o}|\Delta
H|\Psi_{\alpha}^{o}\rangle&\langle\Psi_{\delta}^{o}|\Delta
H|\Psi_{\beta}^{o}\rangle&\langle\Psi_{\delta}^{o}|\Delta
H|\Psi_{\gamma}^{o}\rangle&\langle\Psi_{\delta}^{o}|\Delta
H|\Psi_{\delta}^{o}\rangle\end{pmatrix}.$ (S40)

The overlap integral in the eigenbasis of $|\nu_{i}\rangle$’s is:

$\displaystyle\langle\nu_{3}|\Delta H_{1}^{\dagger}|\nu_{3}\rangle$
$\displaystyle=i\lambda-k_{y}$ $\displaystyle\langle\nu_{3}|\Delta
H_{1}^{\dagger}|\nu_{4}\rangle$ $\displaystyle=-\frac{1}{2}\Delta_{0}$
$\displaystyle\langle\nu_{4}|\Delta H_{1}^{\dagger}|\nu_{4}\rangle$
$\displaystyle=i\lambda+k_{y}$ (S41) $\displaystyle\langle\nu_{3}|\Delta
H_{0}|\nu_{3}\rangle$ $\displaystyle=-k_{y}$
$\displaystyle\langle\nu_{3}|\Delta H_{0}|\nu_{4}\rangle$ $\displaystyle=0$
$\displaystyle\langle\nu_{4}|\Delta H_{0}|\nu_{4}\rangle$
$\displaystyle=k_{y}$ (S42) $\displaystyle\langle\nu_{3}|\Delta
H_{1}|\nu_{3}\rangle$ $\displaystyle=-i\lambda-k_{y}$
$\displaystyle\langle\nu_{3}|\Delta H_{1}|\nu_{4}\rangle$
$\displaystyle=\frac{1}{2}\Delta_{0}$ $\displaystyle\langle\nu_{4}|\Delta
H_{1}|\nu_{4}\rangle$ $\displaystyle=-i\lambda+k_{y}.$ (S43)

And it is useful to note the following sums (and that $\gamma_{j}=\alpha_{j}$
and $\delta_{j}=\beta_{j}$) :

$\displaystyle\sum_{j}\alpha_{j}\alpha_{j}=\sum_{j}\beta_{j}\beta_{j}=1$
$\displaystyle\sum_{j}\alpha_{j}\beta_{j}=\frac{1}{\sqrt{5}}$ (S44)
$\displaystyle\sum_{j}\alpha_{j}\alpha_{j+1}=-\frac{1}{2}$
$\displaystyle\sum_{j}\beta_{j}\beta_{j+1}=-\frac{4}{5}$ (S45)
$\displaystyle\sum_{j}\alpha_{j}\beta_{j+1}=-\frac{2}{\sqrt{5}}$
$\displaystyle\sum_{j}\beta_{j}\alpha_{j+1}=-\frac{1}{2\sqrt{5}}.$ (S46)

The matrix elements are thus calculated as:

$\displaystyle\langle\Psi_{\alpha}^{o}|\Delta H|\Psi_{\alpha}^{o}\rangle$
$\displaystyle=\sum_{j}\alpha_{j}\langle\nu_{3}|(\alpha_{j-1}\Delta
H_{1}^{\dagger}+\alpha_{j}\Delta H_{0}+\alpha_{j+1}\Delta
H_{1})|\nu_{3}\rangle$ $\displaystyle=\sum_{j}\alpha_{j}\alpha_{j+1}(i\lambda-
k_{y})+\alpha_{j}\alpha_{j}(-k_{y})+\alpha_{j}\alpha_{j+1}(-i\lambda-k_{y})$
$\displaystyle=(-k_{y})\sum_{j}\alpha_{j}\alpha_{j}+2\alpha_{j}\alpha_{j+1}$
$\displaystyle=(-k_{y})(1-2\times\frac{1}{2})$ $\displaystyle=0$ (S47)
$\displaystyle\langle\Psi_{\alpha}^{o}|\Delta H|\Psi_{\beta}^{o}\rangle$
$\displaystyle=\langle\Psi_{\alpha}|\Delta
H(\frac{\sqrt{5}}{2}|\Psi_{\beta}\rangle-\frac{1}{2}|\Psi_{\alpha}\rangle)$
$\displaystyle=\frac{\sqrt{5}}{2}\sum_{j}\alpha_{j}\langle\nu_{3}|(\beta_{j-1}\Delta
H_{1}^{\dagger}+\beta_{j}\Delta H_{0}+\beta_{j+1}\Delta
H_{1})|\nu_{3}\rangle-\frac{1}{2}\times 0$
$\displaystyle=\frac{\sqrt{5}}{2}\sum_{j}\alpha_{j+1}\beta_{j}(i\lambda-
k_{y})+\alpha_{j}\beta_{j}(-k_{y})+\alpha_{j}\beta_{j+1}(-i\lambda-k_{y})$
$\displaystyle=\frac{\sqrt{5}}{2}\left[(-\frac{1}{2\sqrt{5}})(i\lambda-
k_{y})+\frac{1}{\sqrt{5}}(-k_{y})+(-\frac{2}{\sqrt{5}})(-i\lambda-
k_{y})\right]$
$\displaystyle=-\frac{1}{4}i\lambda+\frac{1}{4}k_{y}-\frac{1}{2}k_{y}+i\lambda+k_{y}$
$\displaystyle=\frac{3}{4}k_{y}+\frac{3}{4}i\lambda$ (S48)

and

$\displaystyle\langle\Psi_{\alpha}^{o}|\Delta H|\Psi_{\gamma}^{o}\rangle$
$\displaystyle=\sum_{j}\alpha_{j}\langle\nu_{3}|(\alpha_{j-1}\Delta
H_{1}^{\dagger}+\alpha_{j}\Delta H_{0}+\alpha_{j+1}\Delta
H_{1})|\nu_{4}\rangle$
$\displaystyle=\sum_{j}\alpha_{j+1}\alpha_{j}(-\frac{1}{2}\Delta_{0})+\alpha_{j}\alpha_{j}(0)+\alpha_{j}\alpha_{j+1}(\frac{1}{2}\Delta_{0})$
$\displaystyle=0$ (S49) $\displaystyle\langle\Psi_{\alpha}^{o}|\Delta
H|\Psi_{\delta}^{o}\rangle$ $\displaystyle=\langle\Psi_{\alpha}|\Delta
H(\frac{\sqrt{5}}{2}|\Psi_{\delta}\rangle-\frac{1}{2}|\Psi_{\gamma}\rangle)$
$\displaystyle=\frac{\sqrt{5}}{2}\sum_{j}\alpha_{j}\langle\nu_{3}|(\beta_{j-1}\Delta
H_{1}^{\dagger}+\beta_{j}\Delta H_{0}+\beta_{j+1}\Delta
H_{1})|\nu_{4}\rangle-\frac{1}{2}\times 0$
$\displaystyle=\frac{\sqrt{5}}{2}\sum_{j}\alpha_{j+1}\beta_{j}(-\frac{1}{2}\Delta_{0})+\alpha_{j}\beta_{j}(0)+\alpha_{j}\beta_{j+1}(\frac{1}{2}\Delta_{0})$
$\displaystyle=\frac{\sqrt{5}}{2}\left[(-\frac{1}{2\sqrt{5}})(-\frac{1}{2}\Delta_{0})+(-\frac{2}{\sqrt{5}})(\frac{1}{2}\Delta_{0})\right]$
$\displaystyle=-\frac{3}{8}\Delta_{0}$ (S50)

Moving on to the next line in $H_{\textnormal{eff}}$:

$\displaystyle\langle\Psi_{\beta}^{o}|\Delta H|\Psi_{\beta}^{o}\rangle$
$\displaystyle=(\frac{\sqrt{5}}{2}\langle\Psi_{\beta}|-\frac{1}{2}\langle\Psi_{\alpha}|)\Delta
H(\frac{\sqrt{5}}{2}|\Psi_{\beta}\rangle-\frac{1}{2}|\Psi_{\alpha}\rangle)$
$\displaystyle=\frac{5}{4}\langle\Psi_{\beta}|\Delta
H|\Psi_{\beta}\rangle-\frac{\sqrt{5}}{4}\langle\Psi_{\alpha}|\Delta
H|\Psi_{\beta}\rangle-\frac{\sqrt{5}}{4}\langle\Psi_{\beta}|\Delta
H|\Psi_{\alpha}\rangle+\frac{1}{4}\langle\Psi_{\alpha}|\Delta
H|\Psi_{\alpha}\rangle$
$\displaystyle=\frac{5}{4}\sum_{j}\beta_{j}\langle\nu_{3}|(\beta_{j-1}\Delta
H_{1}^{\dagger}+\beta_{j}\Delta H_{0}+\beta_{j+1}\Delta
H_{1})|\nu_{3}\rangle-\frac{1}{2}(\frac{3}{4}k_{y}+\frac{3}{4}i\lambda)-\frac{1}{2}(\frac{3}{4}k_{y}-\frac{3}{4}i\lambda)+0$
$\displaystyle=-\frac{3}{4}k_{y}+\frac{5}{4}\sum_{j}\beta_{j}\beta_{j+1}(i\lambda-
k_{y})+\beta_{j}\beta_{j}(-k_{y})+\beta_{j}\beta_{j+1}(-i\lambda-k_{y})$
$\displaystyle=-\frac{3}{4}k_{y}+\frac{5}{4}\left[(-\frac{4}{5})(i\lambda-
k_{y})+(-k_{y})+(-\frac{4}{5})(-i\lambda-k_{y})\right]$ $\displaystyle=0$
(S51) $\displaystyle\langle\Psi_{\beta}^{o}|\Delta H|\Psi_{\gamma}^{o}\rangle$
$\displaystyle=(\frac{\sqrt{5}}{2}\langle\Psi_{\beta}|-\frac{1}{2}\langle\Psi_{\alpha}|)\Delta
H|\Psi_{\gamma}\rangle$
$\displaystyle=\frac{\sqrt{5}}{2}\langle\Psi_{\beta}|\Delta
H|\Psi_{\gamma}\rangle-\frac{1}{2}\langle\Psi_{\alpha}|\Delta
H|\Psi_{\gamma}\rangle$
$\displaystyle=\frac{\sqrt{5}}{2}\sum_{j}\beta_{j}\langle\nu_{3}|(\alpha_{j-1}\Delta
H_{1}^{\dagger}+\alpha_{j}\Delta H_{0}+\alpha_{j+1}\Delta
H_{1})|\nu_{4}\rangle-0$
$\displaystyle=\frac{\sqrt{5}}{2}\sum_{j}\beta_{j+1}\alpha_{j}(-\frac{1}{2}\Delta_{0})+\beta_{j}\alpha_{j}(0)+\beta_{j}\alpha_{j+1}(\frac{1}{2}\Delta_{0})$
$\displaystyle=\frac{\sqrt{5}}{2}\left[(-\frac{2}{\sqrt{5}})(-\frac{1}{2}\Delta_{0})+(-\frac{1}{2\sqrt{5}})(\frac{1}{2}\Delta_{0})\right]$
$\displaystyle=\frac{3}{8}\Delta_{0}$ (S52)
$\displaystyle\langle\Psi_{\beta}^{o}|\Delta H|\Psi_{\delta}^{o}\rangle$
$\displaystyle=(\frac{\sqrt{5}}{2}\langle\Psi_{\beta}|-\frac{1}{2}\langle\Psi_{\alpha}|)\Delta
H(\frac{\sqrt{5}}{2}|\Psi_{\delta}\rangle-\frac{1}{2}|\Psi_{\gamma}\rangle)$
$\displaystyle=\frac{5}{4}\langle\Psi_{\beta}|\Delta
H|\Psi_{\delta}\rangle-\frac{\sqrt{5}}{4}\langle\Psi_{\beta}|\Delta
H|\Psi_{\gamma}\rangle-\frac{\sqrt{5}}{4}\langle\Psi_{\alpha}|\Delta
H|\Psi_{\delta}\rangle+\frac{1}{4}\langle\Psi_{\alpha}|\Delta
H|\Psi_{\gamma}\rangle$
$\displaystyle=\frac{5}{4}\sum_{j}\beta_{j}\langle\nu_{3}|(\beta_{j-1}\Delta
H_{1}^{\dagger}+\beta_{j}\Delta H_{0}+\beta_{j+1}\Delta
H_{1})|\nu_{4}\rangle-\frac{1}{2}(\frac{3}{8}\Delta_{0})-\frac{1}{2}(-\frac{3}{8}\Delta_{0})+0$
$\displaystyle=\frac{5}{4}\sum_{j}\beta_{j+1}\beta_{j}(-\frac{1}{2}\Delta_{0})+\beta_{j}\beta_{j}(0)+\beta_{j}\beta_{j+1}(\frac{1}{2}\Delta_{0})$
$\displaystyle=0$ (S53)

Moving on to the next line of $H_{\textnormal{eff}}$:

$\displaystyle\langle\Psi_{\gamma}^{o}|\Delta H|\Psi_{\gamma}^{o}\rangle$
$\displaystyle=\sum_{j}\alpha_{j}\langle\nu_{4}|(\alpha_{j-1}\Delta
H_{1}^{\dagger}+\alpha_{j}\Delta H_{0}+\alpha_{j+1}\Delta
H_{1})|\nu_{4}\rangle$
$\displaystyle=\sum_{j}\alpha_{j+1}\alpha_{j}(i\lambda+k_{y})+\alpha_{j}\alpha_{j}(k_{y})+\alpha_{j}\alpha_{j+1}(-i\lambda+k_{y})$
$\displaystyle=k_{y}\sum_{j}\alpha_{j}\alpha_{j}+2\alpha_{j}\alpha_{j+1}$
$\displaystyle=k_{y}(1-2\times\frac{1}{2})$ $\displaystyle=0$ (S54)
$\displaystyle\langle\Psi_{\gamma}^{o}|\Delta H|\Psi_{\delta}^{o}\rangle$
$\displaystyle=\langle\Psi_{\gamma}|\Delta
H(\frac{\sqrt{5}}{2}|\Psi_{\delta}\rangle-\frac{1}{2}|\Psi_{\gamma}\rangle)$
$\displaystyle=\frac{\sqrt{5}}{2}\sum_{j}\alpha_{j}\langle\nu_{4}|(\beta_{j-1}\Delta
H_{1}^{\dagger}+\beta_{j}\Delta H_{0}+\beta_{j+1}\Delta
H_{1})|\nu_{4}\rangle-\frac{1}{2}\times 0$
$\displaystyle=\frac{\sqrt{5}}{2}\sum_{j}\alpha_{j+1}\beta_{j}(i\lambda+k_{y})+\alpha_{j}\beta_{j}(k_{y})+\alpha_{j}\beta_{j+1}(-i\lambda+k_{y})$
$\displaystyle=\frac{\sqrt{5}}{2}\left[(-\frac{1}{2\sqrt{5}})(i\lambda+k_{y})+\frac{1}{\sqrt{5}}(k_{y})+(-\frac{2}{\sqrt{5}})(-i\lambda+k_{y})\right]$
$\displaystyle=-\frac{1}{4}i\lambda-\frac{1}{4}k_{y}+\frac{1}{2}k_{y}+i\lambda-
k_{y}$ $\displaystyle=-\frac{3}{4}k_{y}+\frac{3}{4}i\lambda$ (S55)
$\displaystyle\langle\Psi_{\delta}^{o}|\Delta H|\Psi_{\delta}^{o}\rangle$
$\displaystyle=(\frac{\sqrt{5}}{2}\langle\Psi_{\delta}|-\frac{1}{2}\langle\Psi_{\gamma}|)\Delta
H(\frac{\sqrt{5}}{2}|\Psi_{\delta}\rangle-\frac{1}{2}|\Psi_{\gamma}\rangle)$
$\displaystyle=\frac{5}{4}\langle\Psi_{\delta}|\Delta
H|\Psi_{\delta}\rangle-\frac{\sqrt{5}}{4}\langle\Psi_{\delta}|\Delta
H|\Psi_{\gamma}\rangle-\frac{\sqrt{5}}{4}\langle\Psi_{\gamma}|\Delta
H|\Psi_{\delta}\rangle+\frac{1}{4}\langle\Psi_{\gamma}|\Delta
H|\Psi_{\gamma}\rangle$
$\displaystyle=\frac{5}{4}\sum_{j}\beta_{j}\langle\nu_{4}|(\beta_{j-1}\Delta
H_{1}^{\dagger}+\beta_{j}\Delta H_{0}+\beta_{j+1}\Delta
H_{1})|\nu_{4}\rangle-\frac{1}{2}(-\frac{3}{4}k_{y}-\frac{3}{4}i\lambda)-\frac{1}{2}(-\frac{3}{4}k_{y}+\frac{3}{4}i\lambda)+0$
$\displaystyle=\frac{3}{4}k_{y}+\frac{5}{4}\sum_{j}\beta_{j}\beta_{j+1}(i\lambda+k_{y})+\beta_{j}\beta_{j}(k_{y})+\beta_{j}\beta_{j+1}(-i\lambda+k_{y})$
$\displaystyle=\frac{3}{4}k_{y}+\frac{5}{4}\left[(-\frac{4}{5})(i\lambda+k_{y})+(k_{y})+(-\frac{4}{5})(-i\lambda+k_{y})\right]$
$\displaystyle=0$ (S56)

Collecting all terms together, the effective Hamiltonian is

$\displaystyle H_{\textnormal{eff}}$
$\displaystyle=\begin{pmatrix}0&\frac{3}{4}k_{y}+\frac{3}{4}i\lambda&-\frac{3}{8}\Delta_{0}&0\\\
\frac{3}{4}k_{y}-\frac{3}{4}i\lambda&0&0&\frac{3}{8}\Delta_{0}\\\
-\frac{3}{8}\Delta_{0}&0&0&-\frac{3}{4}k_{y}+\frac{3}{4}i\lambda\\\
0&\frac{3}{8}\Delta_{0}&-\frac{3}{4}k_{y}-\frac{3}{4}i\lambda&0\end{pmatrix}$
(S57)

By taking
$H_{\textnormal{eff}}(k_{y}\rightarrow-i\partial_{y},\lambda\rightarrow y/R)$,
$H_{\textnormal{eff}}$ becomes

$\displaystyle H_{\textnormal{eff}}$
$\displaystyle=\frac{3}{4}\begin{pmatrix}0&-i\partial_{y}+i\frac{1}{R}y&-\frac{1}{2}\Delta_{0}&0\\\
-i\partial_{y}-i\frac{1}{R}y&0&0&\frac{1}{2}\Delta_{0}\\\
-\frac{1}{2}\Delta_{0}&0&0&i\partial_{y}+i\frac{1}{R}y\\\
0&\frac{1}{2}\Delta_{0}&i\partial_{y}-i\frac{1}{R}y&0\end{pmatrix}$ (S58)

and squaring both sides of the time-independent Schrödinger equation yields

$\displaystyle a_{11}$
$\displaystyle=(0,-i\partial_{y}+i\frac{1}{R}y,-\frac{1}{2}\Delta_{0},0)\cdot(0,-i\partial_{y}-i\frac{1}{R}y,-\frac{1}{2}\Delta_{0},0)$
(S59)
$\displaystyle=(-\partial_{y}+\frac{1}{R}y)(\partial_{y}+\frac{1}{R}y)+\frac{1}{4}\Delta_{0}^{2}$
$\displaystyle=-\partial_{y}^{2}-\frac{1}{R}\left(y\partial_{y}-\partial_{y}y\right)+\frac{1}{R^{2}}y^{2}+\frac{1}{4}\Delta_{0}^{2}$
$\displaystyle=-\partial_{y}^{2}-\frac{1}{R}+\frac{1}{R^{2}}y^{2}+\frac{1}{4}\Delta_{0}^{2}$
$\displaystyle a_{14}$
$\displaystyle=(0,-i\partial_{y}+i\frac{1}{R}y,-\frac{1}{2}\Delta_{0},0)\cdot(0,\frac{1}{2}\Delta_{0},i\partial_{y}+i\frac{1}{R}y,0)$
$\displaystyle=\frac{1}{2}\Delta_{0}(-i\partial_{y}+i\frac{1}{R}y-i\partial_{y}-i\frac{1}{R}y)$
$\displaystyle=-i\Delta_{0}\partial_{y}$

And the remaining elements are calculated in a similar way. The resultant
equation is

$\displaystyle\frac{9}{16}$
$\displaystyle\begin{pmatrix}\epsilon_{y,-}&0&0&-i\partial_{y}\Delta_{0}\\\
0&\epsilon_{y,+}&i\partial_{y}\Delta_{0}&0\\\
0&i\partial_{y}\Delta_{0}&\epsilon_{y,+}&0\\\
-i\partial_{y}\Delta_{0}&0&0&\epsilon_{y,-}\end{pmatrix}\psi=E^{2}\psi,$ (S60)

where
$\epsilon_{y,\pm}=-\partial_{y}^{2}+\frac{1}{R^{2}}y^{2}\pm\frac{1}{R}+\frac{1}{4}\Delta_{0}^{2}$.

The eigenstates of this low-energy effective Hamiltonian may then be expressed
as $|\mu_{i}\rangle=\chi_{i}^{\top}|\psi_{i}\rangle$, where $\chi_{i}$
satisfies the constraints of the Hamiltonian matrix representation, and
$|\psi_{i}\rangle=|\psi_{i}(y)\rangle$ is a spatially-varying scalar function
satisfying each of the differential equations contained in the Hamiltonian. We
first identify $\\{\chi_{i}\\}$ ($i\in\\{1,2,3,4\\}$) as eigenvectors of
$\tau_{2}\sigma_{2}$, yielding $|\mu_{i}\rangle$ which take the forms

$\displaystyle|\mu_{1}\rangle$
$\displaystyle=\frac{1}{\sqrt{2}}(1,0,0,1)^{\top}|\psi_{1}\rangle$
$\displaystyle|\mu_{2}\rangle$
$\displaystyle=\frac{1}{\sqrt{2}}(0,-1,1,0)^{\top}|\psi_{2}\rangle$ (S61)
$\displaystyle|\mu_{3}\rangle$
$\displaystyle=\frac{1}{\sqrt{2}}(-1,0,0,1)^{\top}|\psi_{3}\rangle$
$\displaystyle|\mu_{4}\rangle$
$\displaystyle=\frac{1}{\sqrt{2}}(0,1,1,0)^{\top}|\psi_{4}\rangle$

These states serve as zero-energy eigenstates of the effective low-energy
Hamiltonian if the following differential equations are satisfied by the
$\\{\psi_{i}(y)\\}$ (omitting $\frac{1}{4}\Delta_{0}^{2}$ as higher order
terms):

$\displaystyle\left(-\partial_{y}^{2}+\frac{1}{R^{2}}y^{2}-\frac{1}{R}-i\partial_{y}\Delta_{0}\right)|\psi_{1}\rangle$
$\displaystyle=0$ (S62)
$\displaystyle\left(-\partial_{y}^{2}+\frac{1}{R^{2}}y^{2}+\frac{1}{R}-i\partial_{y}\Delta_{0}\right)|\psi_{2}\rangle$
$\displaystyle=0$ (S63)
$\displaystyle\left(-\partial_{y}^{2}+\frac{1}{R^{2}}y^{2}-\frac{1}{R}+i\partial_{y}\Delta_{0}\right)|\psi_{3}\rangle$
$\displaystyle=0$ (S64)
$\displaystyle\left(-\partial_{y}^{2}+\frac{1}{R^{2}}y^{2}+\frac{1}{R}+i\partial_{y}\Delta_{0}\right)|\psi_{4}\rangle$
$\displaystyle=0.$ (S65)

Note that by setting $\Delta_{0}=0$ we recover the differential equation in
the two-band Hopf insulator which only has solution for $-\frac{1}{R}$, so
that we can take this solution and treat $\pm i\partial_{y}\Delta_{0}$ as a
perturbation term so that the resultant wavefunction is still strongly
localized at $y=0$, as shown in Fig. S14(a) and (b). On the other hand, for
the $+\frac{1}{R}$ part which has no normalizable solution in the two-band
Hopf insulator case, there are still no normalizable solutions with the
addition of the $\pm i\partial_{y}\Delta_{0}$ term. The solution instead grows
exponentially with increasing $y$ as shown in Fig. S14(c) and (d).

Figure S14: Numerical solutions for Eq. (S62)-Eq.(S65) with $R=50$ and
$\Delta_{0}=0.1$. In (a), the real and imaginary components of normalized
$\psi_{1}$ and $\frac{\partial\psi_{1}}{\partial y}$ are plotted. Similarly,
(b), (c) and (d) show the real and imaginary components $\psi_{3}$,
$\psi_{2}$, $\psi_{4}$, and their correponding partial derivatives,
respectively. (a) and (b) show that Eq. (S62) and Eq. (S64), which are the
$-\frac{1}{R}$ equations with $\pm i\partial_{y}\Delta_{0}$ term, have
localized solutions at $y=0$, and (c)(d), which are the $+\frac{1}{R}$
equations with $\pm i\partial_{y}\Delta_{0}$ term, show that Eq. (S63) and Eq.
(S65) have no localized solutions as the solution grows exponentially with
increasing $y$.

Therefore, there are zero-energy eigenstates at the inner boundary of the
form,

$\displaystyle|0\rangle_{R}$
$\displaystyle=|\psi_{3}(y)\rangle\begin{pmatrix}1\\\ 0\\\ 0\\\ -1\\\
\end{pmatrix}$ (S66) $\displaystyle|0\rangle_{I}$
$\displaystyle=|\psi_{1}(y)\rangle\begin{pmatrix}i\\\ 0\\\ 0\\\ i\\\
\end{pmatrix}.$ (S67)

Note that the differential equations (S62)-(S65) are obtained from
$|\mu_{1,2,3,4}\rangle$, and would be identical if $i|\mu_{1,2,3,4}\rangle$
were used instead. Therefore, by including a factor of $i$ in
$|\mu_{1}\rangle$, we obtain $|0\rangle_{R(I)}$ as the purely real
(imaginary), $\mathcal{C}^{\prime}$-invariant Majorana zero-modes. More
generally, from numerical diagonalization, we expect two states
$|0\rangle_{\pm}=|0\rangle_{R}\pm|0\rangle_{I}$. Correspondingly,
$|0\rangle_{R(I)}=|0\rangle_{+}+(-)|0\rangle_{-}$.

The above analysis is carried out for $x>0$ i.e. the inner boundary. As for
$x<0$ (the outer boundary), the essential difference is that we swap $H_{1}$
with $H_{1}^{\dagger}$ and $H_{2}$ with $H_{2}^{\dagger}$ in Eqn (S1):

$\begin{pmatrix}H_{0}&H_{1}^{\dagger}&H_{2}^{\dagger}&0&0&0&\cdots\\\
H_{1}&H_{0}&H_{1}^{\dagger}&H_{2}^{\dagger}&0&0&\cdots\\\
H_{2}&H_{1}&H_{0}&H_{1}^{\dagger}&H_{2}^{\dagger}&0&\cdots\\\
0&H_{2}&H_{1}&H_{0}&H_{1}^{\dagger}&H_{2}^{\dagger}&\cdots\\\
\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\ddots\\\
\end{pmatrix}\begin{pmatrix}\psi_{-1}\\\ \psi_{-2}\\\ \psi_{-3}\\\
\psi_{-4}\\\ \vdots\end{pmatrix}=E\begin{pmatrix}\psi_{-1}\\\ \psi_{-2}\\\
\psi_{-3}\\\ \psi_{-4}\\\ \vdots\end{pmatrix}$ (S68)

and repeat the same calculation. This time we get

$\displaystyle\frac{9}{16}$
$\displaystyle\begin{pmatrix}\epsilon_{y,-}&0&0&i\Delta_{0}\partial_{y}\\\
0&\epsilon_{y,+}&-i\Delta_{0}\partial_{y}&0\\\
0&-i\Delta_{0}\partial_{y}&\epsilon_{y,+}&0\\\
i\Delta_{0}\partial_{y}&0&0&\epsilon_{y,-}\end{pmatrix}\psi=E^{2}\psi,$ (S69)

where $\epsilon_{y,\pm}=-\partial_{y}^{2}+\frac{1}{R^{2}}y^{2}\pm\frac{1}{R}$.

Subsequently we will reach the same differential equations (S62) and (S64),
again with eigenstates $|\mu_{1}\rangle$ and $|\mu_{3}\rangle$ in Eqn (S61).
Hence, we similarly find two additional zero-energy states of the form,

$\displaystyle|\bar{0}\rangle_{R}$
$\displaystyle=|\psi_{1}(y)\rangle\begin{pmatrix}1\\\ 0\\\ 0\\\ -1\\\
\end{pmatrix}$ (S70) $\displaystyle|\bar{0}\rangle_{I}$
$\displaystyle=|\psi_{3}(y)\rangle\begin{pmatrix}i\\\ 0\\\ 0\\\ i\\\
\end{pmatrix}.$ (S71)

Thus, if the system has $\mathcal{C}^{\prime}$ symmetry alone, there will be
no localized Majorana zero modes which are equal-weight superposition of
particle and hole components.

In the case of the skyrmion phase in the main text, $|0\rangle_{R(I)}$ is
further constrained by $|\bar{0}\rangle_{I(R)}$, as the linear combinations
$|0\rangle_{R(I)}\pm|\bar{0}\rangle_{I(R)}$ must also respect this particle-
hole symmetry. These constraints are severe enough to enforce a cross-like
structure on each of the four localized eigenstates of the Hamiltonian at
approximately zero energy, such that each consists of two Majorana zero-modes
which are now isolated from one another by a $\pi/2$ phase in the complex
plane, rather than being isolated from one another by separation over real-
space.