# Effective models for dense vortex lattices in the Kitaev honeycomb model

David J. Alspaugh  Sorbonne Université, CNRS, Laboratoire de Physique
Théorique de la Matière Condensée, LPTMC, F-75005 Paris, France  Jean-Noël
Fuchs  Sorbonne Université, CNRS, Laboratoire de Physique Théorique de la
Matière Condensée, LPTMC, F-75005 Paris, France  Anna Ritz-Zwilling  Sorbonne
Université, CNRS, Laboratoire de Physique Théorique de la Matière Condensée,
LPTMC, F-75005 Paris, France  Julien Vidal  Sorbonne Université, CNRS,
Laboratoire de Physique Théorique de la Matière Condensée, LPTMC, F-75005
Paris, France

###### Abstract

We introduce low-energy effective models for dense configurations of vortices
in the Kitaev honeycomb model. Specifically, we consider configurations of
vortices in which vortex-free plaquettes form triangular lattices against a
vortex-full background. Depending on the vortex density, these “dual”
configurations belong to either one of two families classified by translation
and inversion symmetry. As a function of a time-reversal symmetry breaking
term, one family exhibits gapped phases with even Chern numbers separated by
extended gapless phases, while the other exhibits gapped phases with even or
odd Chern numbers, separated by critical points. We construct an effective
model for each family, determine the parameters of these models by fitting the
integrated density of states, and reproduce energy spectra and Chern numbers
of the Kitaev honeycomb model. We also derive phase diagrams and determine
these models’ validity.

##  I Introduction

The characterization of free fermion phases via topological invariants
represents a major achievement in the understanding of quantum matter [1, 2,
3, 4]. Given the spatial dimension, these invariants may be determined through
an analysis of global symmetries such as time-reversal (TRS), particle-hole,
and chiral symmetries, in a classification scheme known as the tenfold way. In
2006, Kitaev introduced an exactly solvable spin-1/2 model on a honeycomb
lattice that can be mapped to a free Majorana-fermion problem in a static
$\mathbb{Z}_{2}$ gauge field [5]. The gauge field results in the possibility
of static $\pi$ fluxes, called vortices, in each plaquette and the problem
separates into independent vortex sectors. In the isotropic limit and in the
vortex-free sector where the ground state lies (see below for details), the
spectrum is gapless, but acquires a gap and a nonzero Chern number $\nu$ when
TRS is broken [5]. To break TRS while preserving the integrability of the
model, Kitaev introduced a three-spin term that results in a next-nearest
neighbor hopping of the Majorana fermions. The corresponding model, closely
related to Haldane’s model [6], has $\nu=\pm 1$ and has been the subject of
many studies (see, e.g., Refs. [7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17]).

From the tenfold way, this model is a class-D topological superconductor in
two dimensions, exhibiting an integer-valued Chern number $\nu\in\mathbb{Z}$
when the energy spectrum is gapped. Notably, Kitaev found that the bulk
excitations of the model behave either as Abelian (even $\nu$) or non-Abelian
(odd $\nu$) anyons, whose properties depend solely on $\nu\mod 16$ [5, 18].
These sixteen topological phases are suggested to occur in the fractional
quantum Hall effect at filling factors with even denominator (such as the
famous $5/2$) [19], and recently proposed honeycomb materials such as
$\alpha$-Li2IrO3 and $\alpha$-RuCl3 have been suggested to exhibit Kitaev-like
quantum spin liquid phases [20].

Figure 1: Examples of vortex arrangements with triangular symmetry: (a)
direct, (b) dual even, and (c) dual odd. The vortex density is $\rho=1/n$ for
the direct configuration and $\rho=(n-1)/n$ for the dual configurations. Grey
hexagons indicate plaquettes that host vortices ($w_{p}=-1$), while white
hexagons indicate their absence ($w_{p}=+1$). Vectors
($\bm{A}_{1}$,$\bm{A}_{2}$) define the geometric unit cell (GUC), whereas
($\bm{a}_{1}$,$\bm{a}_{2}$) define the Hamiltonian unit cell (HUC) compatible
with the periodicity of the $u_{jk}$’s.

In the Kitaev model, the Chern number $\nu$ corresponds to the number of
chiral Majorana edge modes propagating along the boundaries of the system and
is associated with quantized thermal transport [5]. A simple argument due to
Volovik shows that if $\nu$ is odd, an isolated vortex traps an unpaired
Majorana zero mode (MZM) [21]. In contrast, if $\nu$ is even, it traps a
complex fermionic mode at finite energy, similar to the lowest Caroli-de
Gennes-Matricon state in a $\nu=0$ superconductor [22]. This finite-energy
mode can also be described as a pair of coupled MZMs.

Previous studies have analyzed the topological phases arising from triangular
vortex configurations introduced on the honeycomb lattice [12, 13, 14, 15].
These configurations are labeled by a vortex density $\rho=1/n$ with $n$ being
an integer. Because of geometrical constraints, not every integer $n$ is
accessible in a triangular lattice. For instance, $\rho=1/3$ is allowed [see
Fig. 1(a)], but $\rho=1/5$ cannot be realized [8, 15]. In the dilute limit,
i.e., when $n\gg 1$, the spectra of these triangular vortex configurations
only exhibit gapped phases with even Chern numbers, separated by gapless
critical points as the three-spin coupling is varied. These numerical results
are explained by an effective model that describes the hopping of Majorana
fermions on a triangular lattice, whose sites correspond to the locations of
the isolated vortices [12, 13]. Indeed, in the vortex-free sector, the Chern
number of the Kitaev honeycomb model is $\nu=\pm 1$, which implies that an
isolated vortex traps an unpaired MZM [5]. In the following, we will refer to
these triangular vortex configurations as “direct.”

The “duals” to these triangular vortex configurations correspond to dense
vortex configurations, with density $\rho=(n-1)/n$, where the plaquettes that
do not contain vortices form a triangular pattern against a background of
plaquettes that all contain vortices. These dense vortex configurations have
recently been studied numerically in Refs. [14, 15], and can be divided into
two families according to their behavior with respect to translation and
inversion symmetry. We shall refer to these two families as “dual even” (if
$n$ is even such as $\rho=3/4,11/12,15/16,...$) and “dual odd” (if $n$ is odd
such as $\rho=2/3,8/9,11/12,...$) in the following [see Fig. 1(b) and (c)] [8,
15]. The dual even family exhibits gapped phases with even Chern numbers,
separated by extended gapless phases (spin metals with a Majorana-Fermi
surface [23, 24]), whereas the dual odd family has gapped phases with even or
odd Chern numbers, separated by gapless points. These dual configurations are
not describable by the effective model introduced in Refs. [12, 13].

In this article, we construct two different effective models describing the
dual even and dual odd families. These models describe the hopping of pairs of
coupled Majorana fermions on a triangular lattice, where the lattice sites
correspond to the absences of vortices within a vortex sea. These effective
models allow us to understand the variety of Chern numbers and the nature of
gapless phases observed. The article is organized as follows: in Sec. II we
briefly introduce the Kitaev honeycomb model in the Majorana fermion language
and discuss the behaviors of the vortex-free and vortex-full sectors. In Sec.
III we identify three distinct symmetry families for direct and dual vortex
configurations based on translation and inversion symmetries. In Sec. IV we
show that when the low-energy bands of the Kitaev honeycomb model are well-
separated from the high-energy bands, they may be accurately described by an
effective tight-binding model for all three of these families. We first review
the construction of the low-energy model discussed in Refs. [12, 13] for the
direct case, and then introduce new effective models for the dual even and
dual odd cases. To demonstrate the accuracy of these models, we directly
compare their integrated density of states (IDOS) with that of the original
Kitaev honeycomb model. Finally, in Sec. V, we discuss our results, conclude
and give perspectives. Appendices provide details on specific points such as
the validity of the effective models (App. A) and some possible dispersion
relations of vortex bands (App. B).

##  II Kitaev honeycomb model

Figure 2: Standard gauge for the Hamiltonian defined in Eq. (1). Here a black arrow ($J$ hopping) pointing from site $k$ to site $l$ indicates that $u_{lk}=+1$, and consequently $u_{kl}=-1$. The two triangular sublattices of the honeycomb lattice are denoted by black and white colored sites. The dashed pink arrows ($\kappa$ hopping) represent the sign of the product $u_{jl}u_{lk}$ whenever the triplet $(k,l,j)$ is oriented clockwise: Whenever $(k,l,j)$ is a clockwise triplet of sites, and $u_{jl}u_{lk}=-1$, then the dashed pink arrow points from site $j$ to site $k$. |  $\kappa=0$ |  $\kappa\ll J$ |  $\kappa\gg J$ |  $J=0$
---|---|---|---|---
Vortex-free ($\rho=0$) Two bands |  $\Delta=0$ |  $\Delta=6\sqrt{3}\kappa$ $\nu=+1$ |  $\Delta=2J$ $\nu=+1$ |  $\Delta=0$
Vortex-full ($\rho=1$) Four bands |  $\Delta=0$ |  $\Delta\simeq 2\sqrt{3}\kappa$ $\nu=+2$ |  $\Delta\simeq 2\sqrt{3}\kappa$ $\nu=-2$ |  $\Delta=2\sqrt{3}\kappa$ $\nu=-2$
Table 1: Band gap $\Delta$ and Chern number $\nu$ (when $\Delta>0$) of the
vortex-free and vortex-full sectors (backgrounds) of Eq. (1) in the $\kappa\ll
J$ and $\kappa\gg J$ limits.

We begin by recalling some essential facts about Kitaev’s honeycomb model. The
original spin-1/2 model can be reformulated in terms of an effective quadratic
fermionic Hamiltonian given by [5]

$H=\dfrac{i}{4}\sum_{j,k}A_{jk}c_{j}c_{k}.$ (1)

Here $c_{i}$ are Majorana operators defined on the vertices of the honeycomb
lattice such that $c_{i}^{\dagger}=c_{i}$ and
$\\{c_{i},c_{j}\\}=2\delta_{ij}$. The summation runs over all sites of the
lattice, and $A$ is a real skew-symmetric matrix obeying the following rules:

$A_{jk}=\begin{cases}2Ju_{jk},&\text{if $j,k$ are nearest neighbors}\\\
2\kappa u_{jl}u_{lk},&\text{if $(k,l,j)$ is oriented clockwise}\\\
0,&\text{otherwise}\end{cases}.$ (2)

Here the $u_{jk}$ parameters are $\mathbb{Z}_{2}$ gauge variables defined such
that $u_{jk}=-u_{kj}=\pm 1$ for each link of the lattice. We assume here that
the coupling $J$ is the same in all three directions, corresponding to the
isotropic limit of the original spin model. This Hamiltonian describes
noninteracting Majorana fermions coupled to a static $\mathbb{Z}_{2}$ gauge
field. The term $J$ describes the nearest-neighbor hopping, whereas $\kappa$
denotes the strength of the three-spin term that breaks TRS, corresponding to
next-nearest neighbor hopping (see Fig. 2). In the following, without loss of
generality, we will assume that $J$ and $\kappa$ are nonnegative.

For each plaquette $p$ of the honeycomb lattice, we define the
$\mathbb{Z}_{2}$ variable $w_{p}=\prod_{(j,k)\in p}u_{jk}={\rm
e}^{i\phi_{p}}$, where $\phi_{p}$ is the flux threading the plaquette. Here,
for each plaquette, the pair of sites $(j,k)$ is taken such that $k$ always
belongs to one sublattice of the honeycomb lattice, and $j$ belongs to the
other. A value of $w_{p}=-1$ ($+1$) indicates a flux $\phi_{p}=\pi$ ($0$) in
the plaquette $p$, meaning the presence (absence) of a vortex. Two sets of
$\mathbb{Z}_{2}$ link variables $u_{jk}$ defined on the lattice are said to be
equivalent if they correspond to the same configuration of plaquette
variables, or “vortex configuration.” In analyzing Eq. (1), one can then
restrict to a particular vortex sector labelled by the pattern of $w_{p}$’s.
In Fig. 1 we show several examples of vortex configurations, where the
$w_{p}=-1$ plaquettes are shaded grey.

Direct vortex configurations [see Fig. 1 (a)] are defined by two vectors
$\bm{A}_{1}$ and $\bm{A}_{2}$ that indicate the presence of a vortex at
location $\bm{r}=m\,\bm{A}_{1}+n\,\bm{A}_{2}$ for $(m,n)\in\mathbb{Z}^{2}$.
For dual configurations [see Fig. 1 (b) and (c)], we define the vectors
$\bm{A}_{1}$ and $\bm{A}_{2}$ to indicate the absence of a vortex at position
$\bm{r}$. In the following we shall refer to the parallelogram built from the
vectors $\bm{A}_{1}$ and $\bm{A}_{2}$ as the geometric unit cell (GUC). It is
not necessarily the same as the Hamiltonian unit cell (HUC), spanned by the
Bravais lattice vectors $\bm{a}_{1}$ and $\bm{a}_{2}$, and which depends on
the gauge choice for the $u_{jk}$ variables.

We now discuss the physics of the vortex-free sector, for which $w_{p}=+1$ for
all plaquettes ($\rho=0$) [5], and of its dual, the vortex-full sector, where
$w_{p}=-1$ for all plaquettes of the honeycomb lattice ($\rho=1$) [9, 15].
They will constitute the backgrounds to our effective models. Their properties
are summarized in Table 1.

Vortex-free sector. – The standard gauge that realizes this configuration is
represented in Fig. 2. The unit cell is made of one hexagonal plaquette (with
zero flux) so that there are two bands. When $\kappa=0$, the two bands touch
at two Dirac points within the first Brillouin zone (same spectrum as
graphene). Introducing a TRS-breaking term such that $\kappa\ll J$ opens a gap
$\Delta=6\sqrt{3}\kappa$ [5, 6]. The resulting Chern number is $\nu=1$. As
$\kappa$ increases, the system remains gapped, and in the $\kappa\gg J$ limit,
$\Delta=2J$. As can be seen in Fig. 2, when $J=0$ and $\kappa>0$, the two
triangular sublattices of the honeycomb are decoupled and each exhibit an
alternating pattern of $\pm\pi/2$ flux per triangle. The resulting band
structure is gapless along nodal lines.

Translation sym. |  $P_{w}=-1$ |  $P_{w}=-1$ |  $P_{w}=+1$ |  $P_{w}=+1$
---|---|---|---|---
Unit cells |  HUC $=2$ GUC |  HUC $=2$ GUC |  HUC $=$ GUC |  HUC $=$ GUC
Gapped phases with Chern number $\nu$ |  Even $\nu$ |  Even $\nu$ |  Even or odd $\nu$ |  Even or odd $\nu$
Inversion symmetry |  $\mathcal{P}^{2}=-\mathbbm{1}$ |  $\mathcal{P}^{2}=+\mathbbm{1}$ |  $\mathcal{P}^{2}=-\mathbbm{1}$ |  $\mathcal{P}^{2}=+\mathbbm{1}$
Gapless phases |  Point-like |  Extended |  Point-like |  Extended
Vortex configurations and density |  Direct $\rho=1/n$ (Vortex-full $n\to 1$) |  Dual even $\rho=(n-1)/n$, $n$ even |  Dual odd $\rho=(n-1)/n$, $n$ odd (Vortex-free $n\to 1$) |  ?
Table 2: Four families built from translation ($P_{w}=\pm 1$) and inversion
($\mathcal{P}^{2}=\pm\mathbbm{1}$) symmetries. The geometric (GUC) and
Hamiltonian (HUC) unit cells are indicated. The direct and dual configurations
considered in this article occupy three families.

Vortex-full sector. – The simplest gauge choice that achieves this
configuration has a unit cell made of two hexagonal plaquettes (with $\pi$
flux each), which results in four bands. When $\kappa=0$, the energy spectrum
has four gapless Dirac cones [9]. Turning on a small $\kappa$ creates a gap
$\Delta\simeq 2\sqrt{3}\kappa$ and $\nu=2$. At $\kappa=J/2$, the gap closes.
Continuing to increase $\kappa$ re-opens a gap and $\nu=-2$. In the $\kappa\gg
J$ limit, the system remains gapped with $\Delta\simeq 2\sqrt{3}\kappa$. For
$J=0$ and $\kappa>0$, the two triangular sublattices are disconnected. Each
sublattice is similar to a Claro-Wannier model with uniform $\pi/2$ flux per
triangle [25]. This model has two bands of width $\sqrt{3}\kappa$ separated by
a gap $\Delta=2\sqrt{3}\kappa$ and a Chern number $\nu=-1$.

Throughout the rest of this article, unless otherwise stated, we shall define
our energy units such that $J=1$, with the assumption that $J$ is nonzero and
positive.

##  III Translation and inversion symmetries

As shown in Ref. [15], direct and dual vortex configurations can be separated
into three families based on translation and inversion symmetries. We review
here these three families and their properties, which will motivate the
introduction of the effective models.

Whether or not the gauge choice for the $u_{jk}$ variables is invariant under
translations of $\bm{A}_{1}$ and $\bm{A}_{2}$ is found to depend on the parity
of the vortex number $P_{w}$ inside the GUC [15, 14]. When $P_{w}=+1$, there
is an even number of vortices in the GUC. In this case, there exists a gauge
choice invariant under translations of $\bm{A}_{1}$ and $\bm{A}_{2}$ and the
HUC is the same as the GUC [8, 15]. By contrast, when $P_{w}=-1$, there is an
odd number of vortices in the GUC, and the HUC must be double that of the GUC.
Importantly, it has been shown that when $P_{w}=-1$, the system can only
exhibit even-valued Chern numbers, while if $P_{w}=+1$ both even and odd Chern
numbers are possible. The direct vortex configurations always exhibit
$P_{w}=-1$, while the dual configurations can either exhibit $P_{w}=+1$ or
$P_{w}=-1$ [15]. A vortex configuration with $P_{w}=+1$ is shown in Fig. 1(b),
while vortex configurations with $P_{w}=-1$ can be seen in Figs. 1(a) and (c).

We can also analyze how the $u_{jk}$ gauge variables are affected under
application of the inversion symmetry operator
$\mathcal{P}=\mathcal{G}\mathcal{I}$, which is the product of a pure spatial
inversion $\mathcal{I}$ and of a $\mathbb{Z}_{2}$ gauge transformation
$\mathcal{G}$. As discussed in Refs. [15, 14], the square of this operator
behaves as $\mathcal{P}^{2}=\pm\mathbbm{1}$. Extended gapless phases do not
exist for vortex configurations with $\mathcal{P}^{2}=-\mathbbm{1}$. Of the
vortex configurations discussed so far, only the dual even exhibit
$\mathcal{P}^{2}=+\mathbbm{1}$. The direct and the dual odd configurations,
however, both exhibit $\mathcal{P}^{2}=-\mathbbm{1}$.

In Table 2, we classify different families based on the signs of $P_{w}=\pm 1$
and $\mathcal{P}^{2}=\pm\mathbbm{1}$. There are neither direct nor dual
configurations for which $P_{w}=+1$ and $\mathcal{P}^{2}=+\mathbbm{1}$. It was
implicit in the analysis of Ref. [15] that for the vortex configurations
studied there, the gauge transformation can be chosen to have the periodicity
of the Hamiltonian (i.e. that the operator $\mathcal{P}$ commutes with the
translation operators defining the HUC). However, this might not be true for
other vortex configurations on the honeycomb lattice or for the Kitaev model
on other lattices, such as the square-octagon [23] or the triangle-dodecagon
[26, 27, 24, 28].

##  IV Effective low-energy models

Figure 3: Energy versus IDOS of $H$ plotted in blue for: (a) direct
($\rho=1/21$), (b) dual even ($\rho=15/16$), and (c) dual odd ($\rho=20/21$),
all with $\kappa=1$ (energy units are such that $J=1$). In (a), the thick red
line gives the IDOS of the vortex-free background, while for (b) and (c) it
gives that of the vortex-full background. The dotted grey lines indicate the
total bandwidth $\delta_{\rm v}$ of the vortex bands, while the dashed grey
lines indicate the energy gap $\Delta_{\rm b}$ of the background bands. The
second row shows a zoom on the vortex bands (there are two for the direct and
dual odd configurations, and four for the dual even configurations).

Previous works have shown that sufficiently sparse direct vortex
configurations can be described by an effective Majorana tight-binding model
that accurately reproduces their spectra and Chern numbers [12, 13]. The
intuition behind this model follows from analyzing the behavior of isolated
vortices: given a vortex-free background, exchanging $u_{jk}\to-u_{jk}$ for a
single link variable introduces two adjacent vortices into the system. When
these vortices are sufficiently separated from each other, these $w_{p}=-1$
plaquettes each host a single unpaired MZM [5]. The spectrum of a dilute
direct vortex lattice against a vortex-free background can then be understood
as arising from a set of high-energy bands coming from the background
(“background bands”) and a set of low-energy bands coming from the tunneling
of Majorana fermions between these zero modes (“vortex bands”). This can be
seen, for example, in Fig. 3(a) where we plot the IDOS of the Hamiltonian in
Eq. (1) exhibiting a direct vortex configuration with density $\rho=1/21$. The
total Chern number of the system in this case can then be decomposed as
$\nu=\nu_{\rm b}+\nu_{\rm v}$, where $\nu_{\rm b}=1$ and $\nu_{\rm v}$ are the
Chern numbers for the vortex-free background bands and the vortex bands
respectively.

In this section we expand upon the above analysis and derive low-energy
effective models for the two families of dual vortex configurations. These
cases are distinct from the direct configurations as an isolated $w_{p}=+1$
plaquette within a vortex-full background does not host an unpaired MZM.
Instead, in the dual case, the even Chern number of the vortex-full background
implies that an isolated “dual vortex” (i.e. a vortex-free plaquette
surrounded by vortex-full plaquettes) will host a pair of coupled Majorana
fermions equivalent to a single complex fermion mode with finite energy. The
tunneling between these isolated finite-energy modes in the dilute limit then
gives rise to a set of low-energy bands separated from the vortex-full
background (we will still refer to these bands as “vortex bands”). This can be
seen in Figs. 3(b) and (c), where we plot the IDOS for $\rho=15/16$ (dual
even) and $\rho=20/21$ (dual odd). The number of low-energy bands (either two
or four) depends on the symmetry family, as discussed below.

These effective models are valid as long as the total bandwidth $\delta_{\rm
v}$ of the vortex bands is much smaller than the energy gap of the background
bands $\Delta_{\rm b}$, as shown in Fig. 3. For justification and additional
details, see Appendix A.

###  IV.1 Warm-up: direct case

We begin with a review of the construction of the low-energy model introduced
by Lahtinen et al. for the direct vortex configurations [12, 13]. Defining the
HUC vectors as $\bm{a}_{1}=2\bm{A}_{1}$ and $\bm{a}_{2}=\bm{A}_{2}$, we are
then interested in providing an effective model describing Majorana fermions
tunneling between vortices. The simplest corresponding Hamiltonian is given by

$H_{1}=it_{1}\sum_{\langle
jk\rangle}s_{jk}\gamma_{j}\gamma_{k}+it_{2}\sum_{\langle\langle
jk\rangle\rangle}s_{jk}\gamma_{j}\gamma_{k},$ (3)

where $\gamma_{i}$ are Majorana operators defined on the $w_{p}=-1$ plaquettes
of the honeycomb lattice and $(t_{1},t_{2})\in\mathbb{R}^{2}$ represent the
hopping amplitudes between nearest and next-nearest neighbor vortices,
respectively. As explained in Ref. [12], one needs these two terms to
reproduce the Chern numbers observed. The $\mathbb{Z}_{2}$ gauge variables
$s_{jk}=\pm 1$ are chosen such that the total flux accumulated by a Majorana
fermion hopping counter-clockwise on a closed triangular circuit of $t_{1}$
links is $\phi_{1}=\pi/2$, the flux enclosed in a closed circuit of $t_{2}$
links is $\phi_{2}=-\pi/2$, and the flux enclosed in a closed circuit formed
by both $t_{1}$ and $t_{2}$ links is $\phi_{1,2}=\pi/2$ (see Fig 4). In
contrast with Ref. [12], we choose a periodic pattern of $s_{jk}$ variables
that both satisfies these requirements and gives rise to only two sites per
unit cell. Utilizing the gauge choice defined in Fig. 4, we Fourier transform
the Hamiltonian in Eq. (3) to obtain

Figure 4: Direct: Unit cell of the effective model defined in Eq. (3). Grey
and hatched hexagons indicate the two inequivalent sites in the unit cell of
the effective model. The directions of the blue dashed and red dotted arrows
represent the sign of the gauge variables for the nearest ($t_{1}$) and next-
nearest ($t_{2}$) neighbor hopping terms, respectively: If an arrow points
from vortex site $k$ to vortex site $j$, then $s_{jk}=-s_{kj}=+1$. Figure 5:
Direct: (a) Phase diagram and Chern numbers $\nu_{\rm v}$ of the effective
model given in Eq. (3), adapted from Ref. [12]. The phase boundaries are
$t_{2}=t_{1}$, $-t_{1}$, and $-t_{1}/2$. (b) Comparison of the IDOS of the
effective model and that of the vortex bands of $H$ (both IDOS are normalized
between 0 to 1). Here we have $\delta_{\rm v}/\Delta_{\rm b}\simeq 0.009$. The
parameters are $t_{1}=-0.0028$, and $t_{2}=-0.0002$. These parameter values
are indicated by the red dot in the phase diagram of (a), implying that
$\nu_{\rm v}=+1$, in agreement with the microscopic model.

$H_{1}=\sum_{\bm{k}}\Psi_{\bm{k}}^{\dagger}\mathcal{H}_{1}(\bm{k})\Psi_{\bm{k}}.$
(4)

Here $\Psi_{\bm{k}}=(\gamma_{\bm{k}},\overline{\gamma}_{\bm{k}})^{\rm T}$,
such that $\gamma_{i}$ and $\overline{\gamma}_{i}$ represent the Majorana
operators for the solid and striped vortices shown in Fig. 4, respectively,
and

$\displaystyle\begin{gathered}\mathcal{H}_{1}(\bm{k})=\begin{pmatrix}f(\bm{k})&g(\bm{k})\\\
g(\bm{k})^{*}&-f(\bm{k})\end{pmatrix},\\\
f(\bm{k})=-2t_{1}\sin(\bm{k}\cdot\bm{a}_{2})+2t_{2}\sin(\bm{k}\cdot\bm{b}_{1}),\\\
g(\bm{k})=2t_{1}[\sin(\bm{k}\cdot\bm{a}_{1}/2)-i\cos(\bm{k}\cdot\bm{a}_{3})]\\\
\hskip
85.35826pt-2t_{2}[\sin(\bm{k}\cdot\bm{b}_{2})+i\cos(\bm{k}\cdot\bm{b}_{3})].\end{gathered}$
(9)

Here we have defined the nearest neighbor directions as $\bm{a}_{1}/2$,
$\bm{a}_{2}$, and $\bm{a}_{3}=\bm{a}_{2}-\bm{a}_{1}/2$, as well as the next-
nearest neighbor directions as $\bm{b}_{1}=\bm{a}_{1}-\bm{a}_{2}$,
$\bm{b}_{2}=\bm{a}_{2}+\bm{a}_{3}$, and $\bm{b}_{3}=\bm{a}_{2}+\bm{a}_{1}/2$.
The above Majorana Hamiltonian in Eq. (9) can be transformed into a
Bogoliubov-de Gennes Hamiltonian by a $\bm{k}$-independent unitary
transformation, revealing that the system is equivalent to an anisotropic
$p_{x}+ip_{y}$ topological superconductor (see e.g. Sec. 16.3 in Ref. [29]).

In Fig. 5(a) we reproduce the phase diagram of this model from Ref. [12],
which features the possible vortex Chern numbers $\nu_{\rm v}=\pm 1,\pm 3$. We
call $E_{\rm min}(t_{1},t_{2})$ and $E_{\rm max}(t_{1},t_{2})$ the minimum and
maximum positive energy eigenvalues of Eq. (9), which are equal to

$\displaystyle\begin{gathered}E_{\rm
min}(t_{1},t_{2})=\min_{\bm{k}}\sqrt{f(\bm{k})^{2}+|g(\bm{k})|^{2}},\\\ E_{\rm
max}(t_{1},t_{2})=\max_{\bm{k}}\sqrt{f(\bm{k})^{2}+|g(\bm{k})|^{2}}.\end{gathered}$
(12)

Here, we used these extrema to optimize the fitting parameters $t_{1}$ and
$t_{2}$ by imposing the proper bound of the IDOS (they will also be used in
the following section). To give an example of the accuracy of this model, in
Fig. 5(b) we plot the IDOS of the low-energy bands of the Kitaev honeycomb
model in Eq. (1) along with that of the effective model given in Eq. (4).
Here, we fit the parameters $t_{1}$ and $t_{2}$ by matching the IDOS (and
Chern number) corresponding to Eq. (4) with that of Eq. (1). It is possible
that multiple choices of $t_{1}$ and $t_{2}$ may result in the same IDOS as
that of Eq. (1). In this case, as argued in Ref. [12], the relative magnitudes
of these parameters may be predicted by extracting $t_{1}$ and $t_{2}$ from
the energy splitting of MZMs obtained by a numerical calculation with two
spatially-separated vortices in a vortex-free background [12]. While this
method is less precise than exactly matching the IDOS, it has the advantage
that it provides a unique solution. Therefore, both approaches should be
combined.

###  IV.2 Dual even case

In this section we introduce a low-energy model for dual even vortex
configurations. The size of the HUC is double that of the GUC, and we set
$\bm{a}_{1}=2\bm{A}_{1}$, $\bm{a}_{2}=\bm{A}_{2}$, where the vectors
$\bm{A}_{1}$ and $\bm{A}_{2}$ now generate the dual vortex lattice. Each of
these isolated $w_{p}=+1$ plaquettes (i.e. a “dual vortex”) hosts two Majorana
modes at finite energy $\pm\epsilon$. The splitting $2\epsilon$ has been
calculated numerically and plotted in Fig. 6 of Ref. [15]. It vanishes in the
$J=0$ limit in which the two triangular sublattices of the honeycomb lattice
are decoupled (see the discussion on the vortex-full sector in Sec. II). This
means that the two Majorana modes in the same “dual vortex” are uncoupled and
therefore at zero energy. Inspired by this $J=0$ limit and by the effective
model in the direct case [see Eq. (3)], we construct a tight-binding model
that features two “layers” (or flavors) of Majorana modes. They tunnel between
the dual vortices according to the following effective Hamiltonian:

$\displaystyle H_{2}$ $\displaystyle=$ $\displaystyle it_{1}\sum_{\langle
jk\rangle}[s_{jk}^{(\alpha)}\alpha_{j}\alpha_{k}+s_{jk}^{(\beta)}\beta_{j}\beta_{k}]$
(13) $\displaystyle+it_{2}\sum_{\langle\langle
jk\rangle\rangle}[s_{jk}^{(\alpha)}\alpha_{j}\alpha_{k}+s_{jk}^{(\beta)}\beta_{j}\beta_{k}]+i\epsilon\sum_{j}\alpha_{j}\beta_{j}.\qquad$

Here, $\alpha_{i}$ and $\beta_{i}$ are Majorana operators defined on what we
refer to as the “upper” and “lower” layers, respectively. This is
schematically shown in Fig. 6, where the upper layer is colored white and the
lower layer is colored light blue.

Figure 6: Dual even: Unit cell of the effective model introduced in Eq. (14).
Each “dual vortex” hosts two Majorana modes, which we interpret as two
“layers” that we color white and light blue. Here, $\epsilon$ represents the
on-site coupling of two modes on the same site, and the blue dashed and red
dotted arrows represent the sign of the gauge variables for the nearest and
next-nearest neighbor hopping terms, respectively. The pattern of these gauge
variables for both layers is the same as that given in Fig. 4. Figure 7: Dual
even: (a) Phase diagram and Chern numbers $\nu_{\rm v}$ of the effective model
given in Eq. (14), where $t_{1}<0$. $E_{\rm min}(t_{1},t_{2})$ and $E_{\rm
max}(t_{1},t_{2})$ are defined in Eq. (12). (b) Comparison of the IDOS of the
effective model and that of the vortex bands of $H$ (both IDOS are normalized
between 0 to 1). Here, $\delta_{\rm v}/\Delta_{\rm b}\simeq 0.122$ and the
parameters are $\epsilon=0.370$, $t_{1}=-0.052$, and $t_{2}=0.004$. These
parameter values are indicated by the red dot in the phase diagram of (a),
implying that $\nu_{\rm v}=0$, in agreement with the microscopic model.

In this construction, each isolated “dual vortex” in Fig. 1(b) hosts two
Majorana modes: one belonging to $\alpha_{j}$, and one belonging to
$\beta_{j}$. The on-site coupling between them is given by $\epsilon$, and at
the current level we do not allow any inter-layer tunneling between
neighboring dual-vortex sites; nearest and next-nearest neighbor tunneling
occurs only between sites of the same layer. Including such terms would
violate the $\mathcal{P}^{2}=+\mathbbm{1}$ symmetry of the microscopic model.
As before, $s_{jk}^{(\alpha)}=\pm 1$ and $s_{jk}^{(\beta)}=\pm 1$ are
$\mathbb{Z}_{2}$ gauge variables that determine the total enclosed flux of
triangular circuits of dual vortices. From the $P_{w}=-1$ symmetry [which can
be seen in Fig. 1(b)], the total flux enclosed when hopping counter-clockwise
on a nearest neighbor circuit of $t_{1}$ links is $\phi_{1}=-\pi/2$, the flux
enclosed in a next-nearest neighbor circuit of $t_{2}$ links is
$\phi_{2}=\pi/2$, and the flux enclosed in a circuit composed of both $t_{1}$
and $t_{2}$ links is $\phi_{1,2}=-\pi/2$. We may observe that this is simply
the opposite of the flux pattern of the previous section. Noting that
exchanging $\kappa\to-\kappa$ reverses the flux, we may adopt the same pattern
of $s_{jk}$ variables defined in Fig. 4(b) for both layers of Eq. (13). That
is, we set $s_{jk}^{(\alpha)}=s_{jk}^{(\beta)}$. Doing so allows us to Fourier
transform the Hamiltonian of Eq. (13) to find

$H_{2}=\sum_{\bm{k}}\Psi_{\bm{k}}^{\dagger}\begin{pmatrix}\mathcal{H}_{1}(\bm{k})&i\epsilon\tau_{0}\\\
-i\epsilon\tau_{0}&\mathcal{H}_{1}(\bm{k})\end{pmatrix}\Psi_{\bm{k}}.$ (14)

Here,
$\Psi_{\bm{k}}=(\alpha_{\bm{k}},\overline{\alpha}_{\bm{k}},\beta_{\bm{k}},\overline{\beta}_{\bm{k}})^{\rm
T}$, such that $\alpha_{\bm{k}}$ and $\overline{\alpha}_{\bm{k}}$ refer to the
solid and striped white dual vortices, respectively (where solid and striped
indicate the two inequivalent sites within the unit cell), $\beta_{\bm{k}}$
and $\overline{\beta}_{\bm{k}}$ refer to the solid and striped light blue dual
vortices, respectively, $\tau_{0}$ is the identity matrix in the space of
solid and striped dual vortices, and $\mathcal{H}_{1}(\bm{k})$ is defined in
Eq. (9). As can clearly be seen, when $\epsilon=0$ the four energy bands of
Eq. (14) consist of two copies of the bands of the previous model. Due to the
lack of level repulsion resulting from the $\mathcal{P}^{2}=+\mathbbm{1}$
symmetry, these two pairs of bands begin to shift in opposite directions as
$|\epsilon|$ is increased from zero. Whenever $E_{\rm
min}(t_{1},t_{2})<|\epsilon|<E_{\rm max}(t_{1},t_{2})$, two of these bands
will overlap with each other at zero energy, leading to the appearance of a
Majorana-Fermi surface in the Brillouin zone [23, 24]. When $|\epsilon|<E_{\rm
min}(t_{1},t_{2})$, two copies of the valence band of Eq. (9) will have
negative energy, and the vortex Chern number of Eq. (14) will simply be twice
that of Eq. (9). When $|\epsilon|>E_{\rm max}(t_{1},t_{2})$, one valence band
and one conduction band of the model in Eq. (9) will be at negative energy.
These bands each have opposite Chern numbers, and thus the total vortex Chern
number of the system is zero. If we define $\nu^{(1)}$ as the vortex Chern
number of the model defined in Eq. (9), the Chern number of the present model
is thus given by

$\nu_{\rm v}=\begin{cases}2\nu^{(1)},&\text{as $|\epsilon|<E_{\rm
min}(t_{1},t_{2})$}\\\ \text{gapless},&\text{as $|\epsilon|\in[E_{\rm
min}(t_{1},t_{2}),E_{\rm max}(t_{1},t_{2})]$}\\\ 0,&\text{as
$|\epsilon|>E_{\rm max}(t_{1},t_{2})$}\end{cases}.$ (15)

In Fig. 7(a), we plot the phase diagram of Eq. (14) for $t_{1}<0$ ($t_{1}>0$
simply changes the sign of the Chern numbers and mirrors the plot). For all
values of $t_{1}$ and $t_{2}$ we find that we always have $E_{\rm
min}(t_{1},t_{2})<E_{\rm max}(t_{1},t_{2})$ [see Eq. (12)], i.e., there is
always a Majorana-Fermi surface for certain values of $\epsilon$. Because the
vortex-full background has the possible Chern numbers $\nu_{\rm b}=\pm 2$, the
above model thus predicts that the total Chern number can take the values
$\nu=0,\pm 2,\pm 4$, and $\pm 8$. The values $\nu=0,\pm 2,$ and $\pm 4$ have
indeed been observed, but not $\pm 8$ [15].

In Fig. 7(b) we plot the IDOS of the low-energy bands of $H$ along with that
of the effective Hamiltonian $H_{2}$. As before, we fit the parameters
$\epsilon$, $t_{1}$, and $t_{2}$ by matching the IDOS and Chern number with
Eq. (1).

###  IV.3 Dual odd case

Next we introduce a low-energy model for the dual odd configurations.
According to Table 2, in this case the $P_{w}=+1$ symmetry implies that there
is an even number of vortices within each GUC, as can be seen in the example
shown in Fig. 1(c) for $\rho=2/3$. Unlike the previous two families, for these
configurations the HUC and GUC are equivalent to one another, and we therefore
set $\bm{a}_{1}=\bm{A}_{1}$ and $\bm{a}_{2}=\bm{A}_{2}$.

We wish to construct a tight-binding model that features two layers of
Majorana modes that tunnel between dual vortices located. However, in
anticipation of the $\mathcal{P}^{2}=-\mathbbm{1}$ symmetry constraint, this
time we shall allow inter-layer tunneling between neighboring dual-vortex
sites. The Hamiltonian is then

$\displaystyle H_{3}=$ $\displaystyle it_{1}\sum_{\langle
jk\rangle}[s_{jk}^{(\alpha)}\alpha_{j}\alpha_{k}+s_{jk}^{(\beta)}\beta_{j}\beta_{k}]$
(16) $\displaystyle+$ $\displaystyle it_{2}\sum_{\langle\langle
jk\rangle\rangle}[s_{jk}^{(\alpha)}\alpha_{j}\alpha_{k}+s_{jk}^{(\beta)}\beta_{j}\beta_{k}]$
$\displaystyle+$ $\displaystyle it_{\perp 1}\sum_{\langle
jk\rangle}[s_{jk}^{(\perp)}\alpha_{j}\beta_{k}+\widetilde{s}_{jk}^{(\perp)}\beta_{j}\alpha_{k}]$
$\displaystyle+$ $\displaystyle it_{\perp 2}\sum_{\langle\langle
jk\rangle\rangle}[s_{jk}^{(\perp)}\alpha_{j}\beta_{k}+\widetilde{s}_{jk}^{(\perp)}\beta_{j}\alpha_{k}]$
$\displaystyle+$ $\displaystyle i\epsilon\sum_{j}\alpha_{j}\beta_{j}.$

Here, $(t_{1},t_{2})\in\mathbb{R}^{2}$ represent the hopping amplitudes of
modes between nearest and next-nearest neighbor sites in the same layer,
respectively. In contrast, the new terms $(t_{\perp 1},t_{\perp
2})\in\mathbb{R}^{2}$ represent the inter-layer hopping amplitudes between
nearest and next-nearest neighbors, respectively. This Hamiltonian then has
the four sets of gauge variables: $s^{(\alpha)}_{jk}$, $s^{(\beta)}_{jk}$,
$s^{(\perp)}_{jk}$, and $\widetilde{s}^{(\perp)}_{jk}$. To obey the
$\mathcal{P}^{2}=-\mathbbm{1}$ symmetry we must set
$s^{(\beta)}_{jk}=-s^{(\alpha)}_{jk}$, which leaves us with three available
sets of gauge variables to choose in order to analyze the model, along with
the signs of the hopping parameters and $\epsilon$.

Figure 8: Dual odd:  Unit cell of the effective model defined in Eq. (17). (a)
The blue dashed and red dotted arrows represent the sign of the gauge
variables for the nearest $(t_{1})$ and next-nearest neighbor $(t_{2})$ intra-
layer hopping terms, respectively. The upper layer corresponds to the gauge
variables $s^{(\alpha)}_{jk}$, while the lower layer corresponds to the gauge
variables $s^{(\beta)}_{jk}=-s^{(\alpha)}_{jk}$. (b) Schematic of the inter-
layer hopping amplitudes in Eq. (17) (top view). The terms $t_{\perp 1}$ and
$t_{\perp 2}$ are colored blue and red, respectively. Solid arrows refer to
the choice of signs for the $\widetilde{s}^{(\perp)}_{jk}$ gauge variables,
while dashed arrows refer to $s^{(\perp)}_{jk}$. Figure 9: Dual odd: Phase
diagram and Chern numbers $\nu_{\rm v}$ of the model given in Eq. (17). (a)
$t_{1}=0$ and $t_{2}=t_{\perp 2}$. The phase boundaries are $t_{\perp
1}=\epsilon/6,\epsilon,-\epsilon/2$, and $-\epsilon/3$. (b) $t_{2}=0$ and
$t_{1}=t_{\perp 2}$. The phase boundaries are $t_{\perp 1}=\epsilon/6$ and
$\pm\epsilon/2$.

The original honeycomb model in Eq. (1) shows that for this family the
accumulated flux in a triangular circuit of nearest-neighbor dual vortices is
found to be always either $\phi=0$ or $\phi=\pi$. Unlike in the previous two
sections, the gauge variables of Eq. (16) cannot be chosen in such a way that
recreates this pattern. Indeed, for Majorana fermions, triangular circuits
always carry a flux of $\pm\pi/2$. However, as schematically shown in Fig. 8,
the size of the HUC implies that all gauge choices of the $s^{(\alpha)}_{jk}$,
$s^{(\perp)}_{jk}$, and $\widetilde{s}^{(\perp)}_{jk}$ variables result in
patterns of alternating $\pm\pi/2$ flux within the nearest neighbor triangular
circuits, which is the closest a model of Majorana fermions can be to
obtaining a flux of $0$ or $\pi$. As such, we are therefore unrestricted in
our choice of gauge variables. Taking the gauge choices for
$s^{(\alpha)}_{jk}$, $s^{(\perp)}_{jk}$, and $\widetilde{s}^{(\perp)}_{jk}$
defined in Fig. 8, we may Fourier transform the Hamiltonian to find

$H_{3}=\sum_{\bm{k}}\Psi_{\bm{k}}^{\dagger}\mathcal{H}_{3}(\bm{k})\Psi_{\bm{k}},$
(17)

where here, $\Psi_{\bm{k}}=(\alpha_{\bm{k}},\beta_{\bm{k}})^{\rm T}$, and

$\displaystyle\begin{gathered}\mathcal{H}_{3}(\bm{k})=\begin{pmatrix}\widetilde{f}(\bm{k})&\widetilde{g}(\bm{k})\\\
\widetilde{g}(\bm{k})^{*}&-\widetilde{f}(\bm{k})\end{pmatrix},\\\
\widetilde{f}(\bm{k})=2t_{1}\sum_{j=1}^{3}\sin(\bm{k}\cdot\bm{a}_{j})+2t_{2}\sum_{j=1}^{3}\sin(\bm{k}\cdot\bm{b}_{j}),\\\
\widetilde{g}(\bm{k})=i\epsilon+i2t_{\perp
1}[\cos(\bm{k}\cdot\bm{a}_{1})-\cos(\bm{k}\cdot\bm{a}_{2})+\cos(\bm{k}\cdot\bm{a}_{3})]\\\
\hskip 28.45274pt-2t_{\perp
2}[\sin(\bm{k}\cdot\bm{b}_{1})-\sin(\bm{k}\cdot\bm{b}_{2})+\sin(\bm{k}\cdot\bm{b}_{3})].\end{gathered}$
(22)

The nearest neighbor directions are defined as $\bm{a}_{1}$, $\bm{a}_{2}$, and
$\bm{a}_{3}=\bm{a}_{2}-\bm{a}_{1}$, while the next-nearest neighbor directions
are defined as $\bm{b}_{1}=\bm{a}_{1}+\bm{a}_{2}$,
$\bm{b}_{2}=\bm{a}_{2}+\bm{a}_{3}$, and $\bm{b}_{3}=\bm{a}_{3}-\bm{a}_{1}$.
The above Majorana Hamiltonian in Eq. (22) can be transformed into a
Bogoliubov-de Gennes Hamiltonian by a $\bm{k}$-independent unitary
transformation, showing that the model is close to a $p_{x}+ip_{y}$
superconductor on a triangular lattice. This could have been anticipated given
that a “dual vortex” hosts a single complex fermionic mode, that the energy
spectrum must have particle-hole symmetry, and that the gapped bands carry
nonzero Chern numbers.

Figure 10: Dual odd: Comparison of the IDOS of the effective model given in
Eq. (17) and that of the vortex bands of $H$ (both IDOS are normalized between
0 to 1). Here, $\delta_{\rm v}/\Delta_{\rm b}\simeq 0.038$ and the parameters
are $\epsilon=0.395$, $t_{1}=0.095$, $t_{\perp 1}=-0.043$, and all the other
terms zero (note that these values do not correspond to one of the phase
diagrams given in Fig. 9). For an effective model given by these parameters we
calculate the vortex Chern number to be $\nu_{\rm v}=0$, in agreement with the
microscopic model.

In Figs. 9(a) and (b) we plot the phase diagrams of this model as $t_{1}=0$,
$t_{2}=t_{\perp 2}$ and $t_{2}=0$, $t_{1}=t_{\perp 2}$, respectively, which
collectively exhibit the Chern numbers $\nu_{\rm v}=0,\pm 1,\pm 2$, and $\pm
3$. Thus, knowing that the vortex-full background exhibits $\nu_{\rm b}=\pm
2$, the total Chern numbers arising from the above effective model are
$\nu=0,\pm 1,\pm 2,\pm 3,\pm 4$ and $\pm 5$. This is in agreement with what is
observed in Ref. [15], which, for this family, showcases all of these Chern
numbers except $\nu=\pm 4$. We also note that the above model is just one of
several possible choices for the gauge variables $s^{(\alpha)}_{jk}$,
$s^{(\perp)}_{jk}$, and $\widetilde{s}^{(\perp)}_{jk}$, and other choices may
give rise to additional Chern numbers. In total, the model in Eq. (22) is
determined by the five parameters $\epsilon,t_{1},t_{2},t_{\perp 1}$, and
$t_{\perp 2}$, along with the choice of four signs. These are the signs in
front of $\cos(\bm{k}\cdot\bm{a}_{2})$, $\cos(\bm{k}\cdot\bm{a}_{3})$,
$\sin(\bm{k}\cdot\bm{b}_{2})$ and $\sin(\bm{k}\cdot\bm{b}_{3})$ in Eq. (22).
All other signs can be absorbed into the five parameters.

In Fig. 10 we plot the IDOS of the low-energy bands of the Kitaev honeycomb
model in Eq. (1) along with that of the effective model in Eq. (17). The five
parameters $\epsilon$, $t_{1}$, $t_{2}$, $t_{\perp 1}$, and $t_{\perp 2}$ are
found by assuming a relevant gauge and matching the IDOS and Chern number with
Eq. (1). If the Chern number of the model is not known beforehand, then
multiple sets of these parameters may give rise to the same IDOS. In this
case, unlike the previous two symmetry families, there is no clear way to
predict the relative magnitudes of these five parameters by studying a two
“dual vortex” problem numerically.

In Appendix B, we show that not only the IDOS but also the dispersion relation
is well-reproduced by this effective model. As particularly challenging cases,
we consider gapless critical points separating gapped phases with different
Chern numbers.

##  V Discussion and perspectives

In this article, we extended the results of previous works by introducing low-
energy models that can describe both sparse triangular (direct) and dual
vortex configurations in Kitaev’s honeycomb model. In the direct case,
featuring a triangular pattern of vortices against a vortex-free background,
the behavior of the Kitaev honeycomb model may be interpreted as arising from
a combination of the background and of the tunneling of single MZMs between
vortex sites. In the dual case, a triangular pattern of vortex-absences is
instead superimposed on top of a vortex-full background. These isolated “dual
vortices” each contain a pair of coupled Majorana modes, and we have shown
that the properties of the Kitaev honeycomb model may be interpreted as
arising from a combination of the vortex-full background and of two “layers”
of finite-energy Majorana modes that tunnel to other dual-vortex sites in
either an intra-layer or inter-layer fashion.

In both the direct and dual cases the phase diagram of the Kitaev honeycomb
model is controlled by the TRS breaking term $\kappa$ that we allow to be
arbitrarily strong compared to the nearest neighbor interaction $J$. For
certain values of $\kappa/J$, the low-lying energy bands of the microscopic
model are well-separated from the high-energy bands, and in these cases we
have demonstrated that the low-energy bands may be accurately modeled with our
effective descriptions. In addition, the high-energy bands in the microscopic
model resemble the spectra of the vortex-free or vortex-full backgrounds.
Furthermore, for both the direct and dual configurations, we have found that
the total Chern number of the microscopic model $\nu$ is equal to the sum of
the background bands’ Chern number $\nu_{\rm b}$ and the low-energy bands’
Chern number $\nu_{\rm v}$, i.e., $\nu=\nu_{\rm b}+\nu_{\rm v}$. The Chern
numbers predicted by our effective models are in very good agreement with what
has been observed in Ref. [15].

We note that while the parameters of our effective models are found by
directly fitting the IDOS against that of the Kitaev honeycomb model, it would
be worthwhile to develop a deeper connection between the parameters of each
model. The dual odd family, for instance, has a large number of degrees of
freedom with five parameters and four signs of gauge variables to choose from.
Knowing how these parameters depend on the TRS breaking term $\kappa$, for
instance, may help us in better understanding the symmetries and topological
properties of the Kitaev honeycomb model.

Among perspectives, one could think of using the above effective models to
study, in detail, the anyonic excitations predicted by the sixteenfold way
[5]. Preliminary steps in this direction have been taken in Ref. [14]. One
could test the fusion and braiding properties of these excitations similarly
to what has been done in the vortex-free sector in Ref. [10].

In addition, there are several open questions raised by this work that should
be addressed. For example, is it possible to find a phase with Chern number
$\nu=\pm 7$ mod 16 by considering nontriangular vortex configurations [15]?
Are there models that realize the fourth family (see Table 2) of inversion-
symmetric and periodic lattices that is expected to feature extended gapless
phases and even or odd Chern numbers, or is there a general argument showing
that this is impossible? What happens in cases where inversion symmetry does
not commute with the translation operators defining the HUC? More generally,
can we use other crystalline symmetries than inversion to classify all type of
behavior (parity of Chern number, nature of gapless phases) in the Kitaev
model in different vortex sectors and on different lattices (the honeycomb but
also the square-octagon [23] or the triangle-dodecagon lattice [26, 27, 24,
28], for example)?

###### Acknowledgements.

We thank Surajit Basak for many fruitful discussions. D. J. A. acknowledges
LPS Orsay, where the writing of this manuscript was completed, for
hospitality. We acknowledge the financial support of the Quantum Information
Center at Sorbonne Université (QICS) via an Emergence SU grant for the TQCNAA
project and also that of the Domaine d’Intérêt Majeur (DIM) QuanTip of the
Île-de-France region.

##  Appendix A Validity of the effective models

In this Appendix we discuss the validity of our effective models. For a given
vortex configuration, either direct or dual, we may define $\delta_{\rm v}$ as
the total bandwidth of the low-energy vortex bands: two bands for direct or
dual odd, four bands for dual even. We may then define $\Delta_{\rm b}$ as the
energy gap of the next-highest bands. These parameters are defined explicitly
in Fig. 3, where the magnitudes of $\delta_{\rm v}$ and $\Delta_{\rm b}$ are
indicated by the dotted and dashed lines, respectively.

At small $\kappa$, the energy splitting of two isolated MZMs separated by a
distance $d$ is qualitatively given by $t(d,\kappa)\simeq\Delta_{\rm
b}e^{-d/\xi}\cos(k_{\rm F}d+\gamma)$, where $d=|\bm{A}_{1}|=|\bm{A}_{2}|$ is
the vortex separation, $\xi$ is the superconducting coherence length ($\xi\sim
W_{\rm b}/\Delta_{\rm b}$ with $W_{\rm b}$ the total bandwidth of the
background), $k_{F}$ is the distance from the center of the Brillouin zone to
the Dirac points and $\gamma$ is a phaseshift [30, 12]. These Friedel-like
oscillations are a hallmark of the two-vortex solution of isolated MZMs, and
naturally appear in other systems exhibiting Ising anyons, such as the Moore-
Read quantum Hall state, $p_{x}+ip_{y}$ superconductors, and topological
nanowires. The oscillations are related to the change in fermionic parity and
to the fact that MZMs behave as Ising anyons $\sigma$ that have the
possibility to fuse either in the vacuum $1$ or in a fermion $\psi$ according
to the well-known fusion rule $\sigma\times\sigma=1+\psi$ (see Fig. 1 in Ref.
[12]). We may estimate the total vortex bandwidth $\delta_{\rm v}$ for the
direct vortex configurations by $\delta_{\rm
v}\sim|t(d,\kappa)|\sim\Delta_{\rm b}e^{-d/\xi}$. Dilute vortex configurations
have MZMs well-separated from one another, i.e. $d\gg\xi$. This implies that
$\delta_{\rm v}\ll\Delta_{\rm b}$, which we take as a necessary condition for
the validity of the effective model for the direct configurations.

Figure 11: Ratio of the low-energy bandwidth $\delta_{\rm v}$ and the energy
gap $\Delta_{\rm b}$ as a function of $\kappa$ for a direct ($\rho=1/21$) and
dual ($\rho=20/21$) vortex configuration, where $\delta_{\rm v}$ and
$\Delta_{\rm b}$ are defined in Fig. 3. Low-energy tight-binding models are
valid only when $\delta_{\rm v}/\Delta_{\rm b}\ll 1$ as described in the text.
We arbitrarily draw the dotted line at $\delta_{\rm v}/\Delta_{\rm b}=1/2$ to
visually indicate the general regions where these tight-binding models are
most applicable (i.e., hatched regions indicate where they are not
applicable).

For the dual vortex configurations, the low-energy vortex bandwidth
$\delta_{\rm v}$ may be estimated as the energy level splitting $2\epsilon$ of
the two Majorana modes bound to an isolated white plaquette in a vortex-full
sea. Figure 6 of Ref. [15] shows $2\epsilon$ to always be less than the
vortex-full background’s energy gap $\Delta_{\rm b}$, except near $\kappa=0$
and $\kappa=1/2$ where $2\epsilon$ and $\Delta_{\rm b}$ vanish. Thus, the
$\delta_{\rm v}\ll\Delta_{\rm b}$ condition is essentially equivalent to
guaranteeing that $\kappa\neq 0$ and $\kappa\neq 1/2$. Dual vortex
configurations that are “dilute” actually mean that the vortex density
$\rho=(n-1)/n\to 1$.

In Fig. 11, we plot the ratio $\delta_{\rm v}/\Delta_{\rm b}$ for the direct
$\rho=1/21$ and dual odd $\rho=20/21$ vortex configurations as a function of
$\kappa$. We offer these two cases as representative examples, as other
choices of direct and dual configuration densities exhibit the same
qualitative features. As a visual guide, we draw a dotted line at $\delta_{\rm
v}/\Delta_{\rm b}=1/2$ to indicate the regions where $\delta_{\rm
v}/\Delta_{\rm b}\ll 1$. For the direct configuration in Fig. 11(a), we see
that the ratio $\delta_{\rm v}/\Delta_{\rm b}$ is equal to $1/2$ for the two
critical values $\kappa_{c_{1}}$ and $\kappa_{c_{2}}$, represented by the red
circle and triangle, respectively. In Table 1 we see that the energy gap of
the vortex-free ($\rho=0$) background (roughly equal to $\Delta_{\rm b}$)
vanishes as $\kappa=0$ and is of order $1$ as $\kappa\to\infty$. Thus, the
low-energy vortex bands are only well-separated from the background as
$\kappa_{c_{1}}<\kappa<\kappa_{c_{2}}$. For less dense triangular vortex
configurations in the $\rho\to 0$ limit, we find that this region of validity
increases, such that $\kappa_{c_{1}}\to 0$ and $\kappa_{c_{2}}\to\infty$.

For the dual configuration in Fig. 11(b), $\delta_{\rm v}/\Delta_{\rm b}$
equals $1/2$ for the three critical values $\kappa_{c_{i}}$ such that
$i=1,2,3$. These values are represented by the red circle, triangle, and
diamond, respectively. From Sec. II we know that the energy spectrum of the
vortex-full background ($\rho=1$) is gapless as $\kappa=0$ and $\kappa=1/2$.
This establishes two regions of validity for any low-energy effective theory
that we wish to build: a small region of low $\kappa$ values such that
$\kappa_{c_{1}}<\kappa<\kappa_{c_{2}}$, and a bigger region of large $\kappa$
values such that $\kappa>\kappa_{c_{3}}$. For denser dual configurations in
the $\rho\to 1$ limit, we find that these two regions of validity increase,
such that $\kappa_{c_{1}}\to 0$, $\kappa_{c_{2}}\to 1/2^{-}$, and
$\kappa_{c_{3}}\to 1/2^{+}$.

##  Appendix B Dispersion relations of vortex bands

In this Appendix, we analyze three particular cases of quantum phase
transitions occurring in the dual odd family of vortex configurations. These
transitions are identified by the number of contact points present within the
first Brillouin zone, exhibiting one, two, or six such points. The aim is here
to show that the effective model captures not only the integrated density of
states but really the $\bm{k}$-resolved dispersion relation (at least
qualitatively).

###  B.1 One contact point

Figure 12: (a) Dispersion of the low-lying vortex bands of the Kitaev
honeycomb model for the $\rho=12/13$ dual odd vortex configuration as
$\kappa_{c}=0.219461$. At this point, the vortex Chern number changes from
$\nu_{\rm v}=-1$ to $\nu_{\rm v}=0$. (b) Dispersion of the effective model of
Eq. (23). The parameter values are given by $\epsilon=0.124$, $t_{1}=0.0186$,
$t_{2}=0$, $t_{\perp 2}=-0.0124$, and $t_{\perp 1}=\epsilon/6$. The Brillouin
zone for both plots is given by
$\bm{k}=k_{1}\bm{a}_{1}^{*}+k_{2}\bm{a}_{2}^{*}$, where
$(\bm{a}_{1}^{*},\bm{a}_{2}^{*})$ are the reciprocal lattice vectors of
$(\bm{a}_{1},\bm{a}_{2})$ respectively, and $|k_{1,2}|\leq 1/2$.

This type of transition occurs when the vortex Chern number changes from
$\nu_{\rm v}=-1$ to $\nu_{\rm v}=0$ for the dual odd vortex configurations of
density $\rho=(n-1)/n$ such that $n$ is not a multiple of three. That is, for
densities such as $\rho=12/13$, $18/19$, $24/25$, and etc. These transitions
generally occur when $\kappa_{c}\sim 0.1-0.2$. Before the transition
($\kappa<\kappa_{c}$), the system features two gapped Dirac points of opposite
chirality and is a Chern insulator with $\nu_{\rm v}=-1$. As $\kappa$
increases, the Dirac points move and at a critical $\kappa_{c}$ they meet at a
time-reversal invariant momentum (TRIM) where they annihilate and the gap
closes. Exactly at the transition, there is a single contact point that does
not carry a topological charge. After this merging transition, the gap reopens
and the system becomes a trivial insulator with $\nu_{\rm v}=0$. This merging
transition is similar to the one occurring in the Kitaev honeycomb model when
evolving from the isotropic point $J_{x}=J_{y}=J_{z}$ to a very anisotropic
case [5], or also in uniaxially-strained graphene [31].

The Hamiltonian has the same form as that given in Eq. (22), except that the
function $\widetilde{g}(\bm{k})$ is modified. Focusing on the $\rho=12/13$
vortex density as a example, the “one contact point” Hamiltonian may be
expressed as:

$\displaystyle\widetilde{g}_{1}(\bm{k})$ $\displaystyle=$
$\displaystyle-i2t_{\perp
1}[\cos(\bm{k}\cdot\bm{a}_{1})+\cos(\bm{k}\cdot\bm{a}_{2})+\cos(\bm{k}\cdot\bm{a}_{3})]$
(23) $\displaystyle-$ $\displaystyle 2t_{\perp
2}[\sin(\bm{k}\cdot\bm{b}_{1})+\sin(\bm{k}\cdot\bm{b}_{2})+\sin(\bm{k}\cdot\bm{b}_{3})]$
$\displaystyle+$ $\displaystyle i\epsilon.$

We note that this Hamiltonian results from a different gauge choice for the
$s_{jk}^{(\alpha)}$, $s_{jk}^{(\perp)}$, and $\widetilde{s}_{jk}^{(\perp)}$
variables than what is shown in Fig. 8. In Fig. 12 we plot the dispersion of
both this effective model and the low-lying vortex bands of the Kitaev
honeycomb model of the $\nu_{\rm v}=-1\to 0$ transition in the $\rho=12/13$
dual odd vortex configuration as $\kappa_{c}=0.219461$ [15]. The proper
parameter values for the effective model are given in the caption. The
critical point of the effective model appears when $t_{\perp 1}=\epsilon/6$.
As $t_{\perp 1}<\epsilon/6$ the Chern number of the effective model is
$\nu_{\rm v}=-1$, and when $t_{\perp 1}>\epsilon/6$ the Chern number is
$\nu_{\rm v}=0$. The sign of $t_{\perp 2}$ controls the sign of $\nu_{\rm v}$,
and the choice of signs of the cosine terms of $\widetilde{g}_{1}(\bm{k})$
controls the position of the contact point within the Brillouin zone. If we
were to instead analyze a different vortex density exhibiting a one contact
point transition, one may have to alter the signs of the cosine terms in order
to properly match the location of the contact point within the Brillouin zone
(this is true for the two and six contact points models as well). Note that
because the value of $\delta_{\rm v}/\Delta_{\rm b}\simeq 0.38$ is relatively
high (i.e. the vortex bands are not so well separated from the background
bands), the bandwidth of the effective model overestimates that of the low-
lying vortex bands of the Kitaev honeycomb model.

###  B.2 Two contact points

Figure 13: (a) Dispersion of the low-lying vortex bands of the Kitaev
honeycomb model for the $\rho=20/21$ dual odd vortex configuration as
$\kappa_{c}=0.106486$. At this point, the vortex Chern number changes from
$\nu_{\rm v}=-2$ to $\nu_{\rm v}=0$. (b) Dispersion of the effective model of
Eq. (24). The parameter values are given by $\epsilon=0.1$, $t_{1}=0$,
$t_{2}=t_{\perp 2}=0.0009$, and $t_{\perp 1}=-\epsilon/3$. The Brillouin zone
for both plots is given by $\bm{k}=k_{1}\bm{a}_{1}^{*}+k_{2}\bm{a}_{2}^{*}$
with $|k_{1,2}|\leq 1/2$.

This type of transition results in a change of Chern number from $\nu_{\rm
v}=-2$ to $\nu_{\rm v}=0$, and occurs when $\kappa_{c}\sim 0.1-0.2$ for the
dual odd vortex configurations which have a multiple of three in the
denominator of their vortex density. That is, for $\rho=8/9$, $20/21$,
$26/27$, and etc. Before the transition, there are four gapped Dirac cones
(two of positive and two of negative chirality) and the system is a Chern
insulator with $\nu_{\rm v}=-2$. At the transition, Dirac points of opposite
chirality meet pairwise, merge, and the gap closes. Exactly at the transition,
there are two contact points that are not protected by a topological charge
and which do not occur at a TRIM. This is possible since this transition may
be described as two merging transitions occurring simultaneously. After the
transition, a merging gap opens and the system becomes a trivial insulator
with $\nu_{\rm v}=0$.

The Hamiltonian is similar to that of Eq. (22) except for a modified
$\widetilde{g}(\bm{k})$ term. Focusing on the $\rho=20/21$ case as an example,
the “two contact points ” Hamiltonian is given by:

$\displaystyle\widetilde{g}_{2}(\bm{k})$ $\displaystyle=$ $\displaystyle
i2t_{\perp
1}[\cos(\bm{k}\cdot\bm{a}_{1})+\cos(\bm{k}\cdot\bm{a}_{2})-\cos(\bm{k}\cdot\bm{a}_{3})]$
(24) $\displaystyle-$ $\displaystyle 2t_{\perp
2}[\sin(\bm{k}\cdot\bm{b}_{1})-\sin(\bm{k}\cdot\bm{b}_{2})+\sin(\bm{k}\cdot\bm{b}_{3})]$
$\displaystyle+$ $\displaystyle i\epsilon.$

Figure 14: (a) Dispersion of the low-lying vortex bands of the Kitaev
honeycomb model for the $\rho=6/7$ dual odd vortex configuration as
$\kappa_{c}=1.80057$. At this point, the vortex Chern number changes from
$\nu_{\rm v}=-3$ to $\nu_{\rm v}=3$. (b) Dispersion of the effective model of
Eq. (25). The parameter values are given by $\epsilon=t_{\perp 1}=0.22$,
$t_{1}=0.11$, and $t_{\perp 2}=t_{2}=0$. The Brillouin zone for both plots is
given by $\bm{k}=k_{1}\bm{a}_{1}^{*}+k_{2}\bm{a}_{2}^{*}$ with $|k_{1,2}|\leq
1/2$.

In Fig. 13 we plot the dispersion of both this effective model and the low-
lying vortex bands of the Kitaev honeycomb model of the $\nu_{\rm v}=-2\to 0$
transition in the $\rho=20/21$ dual odd vortex configuration, as
$\kappa_{c}=0.106486$. The proper parameter values for the effective model are
given in the caption. The critical points of the effective model appear when
$t_{\perp 1}=-\epsilon/3$. As $t_{\perp 1}>-\epsilon/3$ the Chern number of
the effective model is $\nu_{\rm v}=0$, and when $t_{\perp 1}<-\epsilon/3$ the
Chern number is $\nu_{\rm v}=-2$. The sign of $t_{2}$ controls the sign of
$\nu_{\rm v}$. Note that $\delta_{\rm v}/\Delta_{\rm b}\simeq 0.62$, and hence
the bandwidth of the effective model is slightly larger than that of the low-
lying vortex bands of the Kitaev honeycomb model.

We note that in this system, the energy spectrum is gapped and topologically
nontrivial only if the sign of one of the sine terms in
$\widetilde{g}_{2}(\bm{k})$ is dissimilar to the sign of the other two. If
every sine term were positive, for example, then the $2\times 2$ matrix
Hamiltonian may be expressed as the linear combination of only two Pauli
matrices and would either be gapless before the transition (i.e. a topological
semi-metal with gapless Dirac points) or gapped after the transition with a
trivial Chern number of $\nu_{\rm v}=0$, due to a merging transition.

###  B.3 Six contact points

This last type of transition results in a change of Chern number from
$\nu_{\rm v}=-3$ to $\nu_{\rm v}=3$, and only occurs for the $\rho=6/7$ dual
odd vortex configuration as $\kappa_{c}\simeq 1.8$. In contrast to the two
previous scenarios, it is not a merging transition of Dirac points, but six
Dirac points that each have a simultaneous change in sign of their mass (or
gap). The Hamiltonian is similar in form to the previous two except for the
$\widetilde{g}(\bm{k})$ term:

$\displaystyle\widetilde{g}_{6}(\bm{k})$ $\displaystyle=$ $\displaystyle
i2t_{\perp
1}[\cos(\bm{k}\cdot\bm{a}_{1})-\cos(\bm{k}\cdot\bm{a}_{2})+\cos(\bm{k}\cdot\bm{a}_{3})]$
(25) $\displaystyle-$ $\displaystyle 2t_{\perp
2}[\sin(\bm{k}\cdot\bm{b}_{1})+\sin(\bm{k}\cdot\bm{b}_{2})+\sin(\bm{k}\cdot\bm{b}_{3})]$
$\displaystyle+$ $\displaystyle i\epsilon.$

In Fig. 14 we plot the dispersion of both this effective model and the low-
lying vortex bands of the Kitaev honeycomb model of the $\nu_{\rm v}=-3\to 3$
transition in the $\rho=6/7$ dual odd vortex configuration, as
$\kappa_{c}=1.80057$. The proper parameter values for the effective model are
given in the caption. The critical points of the effective model appear when
$t_{\perp 2}=0$. As $t_{\perp 2}>0$ the Chern number of the effective model is
$\nu_{\rm v}=-3$, and when $t_{\perp 2}<0$ the Chern number is $\nu_{\rm
v}=3$. The ratio $\delta_{\rm v}/\Delta_{\rm b}\simeq 0.38$ indicates that the
bandwidth of the effective model is slightly larger than that of the low-lying
vortex bands of the Kitaev honeycomb model. Note that because $t_{1}\neq
t_{2}=0$, we do not have the same issue of the sine’s signs as discussed in
the two contact points section, and the $2\times 2$ matrix Hamiltonian is
still composed of three different Pauli matrices.

## References

  * Schnyder  _et al._ [2008] A. P. Schnyder, S. Ryu, A. Furusaki, and A. W. W. Ludwig, Classification of topological insulators and superconductors in three spatial dimensions, Phys. Rev. B 78, 195125 (2008).
  * Kitaev [2009] A. Kitaev, Periodic table for topological insulators and superconductors, AIP Conf. Proc. 1134, 22 (2009).
  * Ryu  _et al._ [2010] S. Ryu, A. P. Schnyder, A. Furusaki, and A. W. W. Ludwig, Topological insulators and superconductors: tenfold way and dimensional hierarchy, New J. Phys. 12, 065010 (2010).
  * Chiu  _et al._ [2016] C.-K. Chiu, J. C. Y. Teo, A. P. Schnyder, and S. Ryu, Classification of topological quantum matter with symmetries, Rev. Mod. Phys. 88, 035005 (2016).
  * Kitaev [2006] A. Kitaev, Anyons in an exactly solved model and beyond, Ann. Phys. (N. Y.) 321, 2 (2006).
  * Haldane [1988] F. D. M. Haldane, Model for a Quantum Hall Effect without Landau Levels: Condensed-Matter Realization of the “Parity Anomaly”, Phys. Rev. Lett. 61, 2015 (1988).
  * Lahtinen  _et al._ [2008] V. Lahtinen, G. Kells, A. Carollo, T. Stitt, J. Vala, and J. K. Pachos, Spectrum of the non-abelian phase in Kitaev’s honeycomb lattice model, Ann. Phys. (N. Y.) 323, 2286 (2008).
  * Kamfor  _et al._ [2011] M. Kamfor, S. Dusuel, K. P. Schmidt, and J. Vidal, Fate of Dirac points in a vortex superlattice, Phys. Rev. B 84, 153404 (2011).
  * Lahtinen and Pachos [2010] V. Lahtinen and J. K. Pachos, Topological phase transitions driven by gauge fields in an exactly solvable model, Phys. Rev. B 81, 245132 (2010).
  * Lahtinen [2011] V. Lahtinen, Interacting non-Abelian anyons as Majorana fermions in the honeycomb lattice model, New J. Phys. 13, 075009 (2011).
  * Kells  _et al._ [2011] G. Kells, J. Kailasvuori, J. K. Slingerland, and J. Vala, Kaleidoscope of topological phases with multiple Majorana species, New J. Phys. 13, 095014 (2011).
  * Lahtinen  _et al._ [2012] V. Lahtinen, A. W. W. Ludwig, J. K. Pachos, and S. Trebst, Topological liquid nucleation induced by vortex-vortex interactions in Kitaev’s honeycomb model, Phys. Rev. B 86, 075115 (2012).
  * Lahtinen  _et al._ [2014] V. Lahtinen, A. W. W. Ludwig, and S. Trebst, Perturbed vortex lattices and the stability of nucleated topological phases, Phys. Rev. B 89, 085121 (2014).
  * Zhang  _et al._ [2020] S.-S. Zhang, C. D. Batista, and G. B. Halász, Toward Kitaev’s sixteenfold way in a honeycomb lattice model, Phys. Rev. Research 2, 023334 (2020).
  * Fuchs  _et al._ [2020] J.-N. Fuchs, S. Patil, and J. Vidal, Parity of Chern numbers in the Kitaev honeycomb model and the sixteenfold way, Phys. Rev. B 102, 115130 (2020).
  * Grushin and Repellin [2023] A. G. Grushin and C. Repellin, Amorphous and Polycrystalline Routes toward a Chiral Spin Liquid, Phys. Rev. Lett. 130, 186702 (2023).
  * Cassella  _et al._ [2023] G. Cassella, P. D’Ornellas, T. Hodson, W. M. H. Natori, and J. Knolle, An Exact Chiral Amorphous Spin Liquid, Nat. Commun. 14, 6663 (2023).
  * [18] B. A. Bernevig and T. Neupert, in Topological Aspects of Condensed Matter Physics: Lecture Notes of Les Houches 2014 (Oxford University Press, Oxford, 2017), arXiv:1506.05805.
  * Ma and Feldman [2019] K. K. W. Ma and D. E. Feldman, The sixteenfold way and the quantum Hall effect at half-integer filling factors, Phys. Rev. B 100, 035302 (2019).
  * Trebst and Hickey [2022] S. Trebst and C. Hickey, Kitaev materials, Phys. Rep. 950, 1 (2022).
  * Volovik [1999] G. E. Volovik, Fermion zero modes on vortices in chiral superconductors, JETP Lett. 70, 609 (1999).
  * Caroli  _et al._ [1964] C. Caroli, P.-G. de Gennes, and J. Matricon, Bound Fermion states on a vortex line in a type II superconductor, Physics Letters 9, 307 (1964).
  * [23] G. Baskaran, G. Santhosh, and R. Shankar, Exact quantum spin liquids with Fermi surfaces in spin-half models, arXiv:0908.1614v3.
  * Tikhonov and Feigel’man [2010] K. S. Tikhonov and M. V. Feigel’man, Quantum Spin Metal State on a Decorated Honeycomb Lattice, Phys. Rev. Lett. 105, 067207 (2010).
  * Claro and Wannier [1979] F. H. Claro and G. H. Wannier, Magnetic subband structure of electrons in hexagonal lattices, Phys. Rev. B 19, 6068 (1979).
  * Yao and Kivelson [2007] H. Yao and S. A. Kivelson, Exact Chiral Spin Liquid with Non-Abelian Anyons, Phys. Rev. Lett. 99, 247203 (2007).
  * Dusuel  _et al._ [2008] S. Dusuel, K. P. Schmidt, J. Vidal, and R. L. Zaffino, Perturbative study of the Kitaev model with spontaneous time-reversal symmetry breaking, Phys. Rev. B 78, 125102 (2008).
  * Yao and Lee [2011] H. Yao and D.-H. Lee, Fermionic Magnons, Non-Abelian Spinons, and the Spin Quantum Hall Effect from an Exactly Solvable Spin-1/2 Kitaev Model with SU(2) Symmetry, Phys. Rev. Lett. 107, 087205 (2011).
  * [29] Topological Insulators and Topological Superconductors, B. A. Bernevig and T. L. Hughes, 1st ed. (Princeton University Press, 2013).
  * Cheng  _et al._ [2009] M. Cheng, R. M. Lutchyn, V. Galitski, and S. D. Sarma, Splitting of Majorana-Fermion Modes due to Intervortex Tunneling in a ${p}_{x}+i{p}_{y}$ Superconductor, Phys. Rev. Lett. 103, 107001 (2009).
  * G. Montambaux, F. Piéchon, J.-N. Fuchs and M.O. Goerbig [2009] G. Montambaux, F. Piéchon, J.-N. Fuchs and M.O. Goerbig, Merging of Dirac points in a two-dimensional crystal, Phys. Rev. B 80, 153412 (2009).