# Fractional charge bound to a vortex in two-dimensional topological
crystalline insulators

Eunwoo Lee  Department of Physics and Astronomy, Seoul National University,
Seoul 08826, Korea  Center for Correlated Electron Systems, Institute for
Basic Science (IBS), Seoul 08826, Korea  Center for Theoretical Physics (CTP),
Seoul National University, Seoul 08826, Korea  Akira Furusaki  RIKEN Center
for Emergent Matter Science, Wako, Saitama, 351-0198, Japan  Condensed Matter
Theory Laboratory, RIKEN, Wako, Saitama, 351-0198, Japan  Bohm-Jung Yang
<EMAIL_ADDRESS>Department of Physics and Astronomy, Seoul National
University, Seoul 08826, Korea  Center for Correlated Electron Systems,
Institute for Basic Science (IBS), Seoul 08826, Korea  Center for Theoretical
Physics (CTP), Seoul National University, Seoul 08826, Korea

###### Abstract

We establish the correspondence between the fractional charge bound to a
vortex in a textured lattice and the relevant bulk band topology in two-
dimensional (2D) topological crystalline insulators. As a representative
example, we consider the Kekule textured graphene whose bulk band topology is
characterized by a 2D $\mathbb{Z}_{2}$ topological invariant $\nu_{\rm 2D}$
protected by inversion symmetry. The fractional charge localized at a vortex
in the Kekule texture is shown to be related to the change in the bulk
topological invariant $\nu_{\rm 2D}$ around the vortex, as in the case of the
Su-Schriefer-Heeger model in which the fractional charge localized at a domain
wall is related to the change in the bulk charge polarization between
degenerate ground states. We show that the effective three-dimensional (3D)
Hamiltonian, where the angle $\theta$ around a vortex in Kekule-textured
graphene is a third coordinate, describes a 3D axion insulator with a
quantized magnetoelectric polarization. The spectral flow during the adiabatic
variation of $\theta$ corresponds to the chiral hinge modes of an axion
insulator and determines the accumulated charge localized at the vortex, which
is half-quantized when chiral symmetry exists. When chiral symmetry is absent,
electric charge localized at the vortex is no longer quantized, but the vortex
always carries a half-quantized Wannier charge as long as inversion symmetry
exists. For the cases when magnetoelectric polarization is quantized due to
the presence of symmetry that reverses the space-time orientation, we classify
all possible topological crystalline insulators whose vortex defect carries a
fractional charge.

##  I Introduction

Fractional charge localized at topological defects is closely related to the
topology of bulk electronic states Su et al. (1980); Jackiw and Rebbi (1976);
Jackiw and Rossi (1981); Goldstone and Wilczek (1981); Read and Green (2000);
Teo and Kane (2010a, b). For instance, the one-dimensional (1D) model proposed
by Su, Schrieffer, and Heeger (SSH) is a representative system where a half
electric charge is localized at a domain wall Su et al. (1980). Spontaneous
formation of lattice dimerization gives rise to two degenerate ground states,
and the zero-mode carrying the half electric charge is localized at the domain
wall interpolating between the degenerate ground states. The inherent
relationship between the zero-mode charge and the topology of the bulk
electronic states can be seen from the quantized charge polarization $P_{1}$
of the two degenerate ground states, given by $P_{1}=0$ and $P_{1}=1/2$,
respectively. The fractional charge accumulated at the domain wall is
determined by the difference of the bulk charge polarization $P_{1}$ of the
two ground states.

A remarkable idea realizing fractional charge in two dimensions is proposed by
Hou, Chamon, and Mudry (HCM) in Ref. [Hou et al., 2007]. They considered a
graphenelike system, whose low-energy excitations are described by massless
Dirac fermions, and showed that a vortex in the order parameter for the
Kekule-type dimerization accommodates zero-energy bound states with fractional
charge. Spontaneous formation of the Kekule texture in graphene leads to the
degenerate ground states with broken lattice symmetries Hou et al. (2007);
Chamon et al. (2008a, b); Ryu et al. (2009). In view of the connection between
the charge fractionalization and the topological properties (quantized
polarization) of the degenerate ground states in the SSH model, the charge
fractionalization found in the HCM model naturally leads to the following
question: Is it possible to understand the fractional charge at the vortex
from topological properties of the bulk electronic states in two dimensions?

In this paper, we show that the ground states of Kekule textured graphene
(KTG) are characterized by a $\mathbb{Z}_{2}$ topological invariant $\nu_{\rm
2D}$, which is quantized in two-dimensional (2D) systems with inversion
symmetry $P$ Hughes et al. (2011); Morimoto and Furusaki (2014); Qi et al.
(2008). Two insulating ground states with distinct quantized values of
$\nu_{\rm 2D}$ can be changed by switching the strong and weak bonds in the
Kekule texture, as in the case of the SSH model where two gapped phases with
distinct $P_{1}$ can be interchanged in an analogous way. We show that KTG
with $\nu_{\rm 2D}=1$ is a 2D inversion-symmetric second-order topological
insulator (SOTI) with zero-energy corner states. The nontrivial bulk band
topology of KTG is confirmed by calculating the parity eigenvalues, the Wilson
loop spectra, and the corner charge distribution. We note that the higher-
order band topology of KTG related with sixfold rotation symmetry was
discussed in Refs. Noh et al. (2018); Benalcazar et al. (2018). Similar to the
case of SSH model, the charge accumulation at a vortex in the order parameter
in the HCM model is related to the difference of the bulk topological quantity
$\nu_{\rm 2D}$ between the two degenerate ground states with $\nu_{\rm 2D}=0$
and $\nu_{\rm 2D}=1$, respectively. Furthermore, we demonstrate that the
change in $\nu_{\rm 2D}$ around a vortex is related to the quantized
magnetoelectric polarizability $P_{3}$ of an axion insulator, which guarantees
the presence of a zero-mode state localized at the vortex core.

Finally, considering that $P_{3}$ is quantized in the presence of the symmetry
reversing the space-time orientation, we classify all possible topological
crystalline insulators where an order parameter vortex can support fractional
vortex charge. It is worth noting that the correspondence between the bulk
band topology and the vortex bound state in KTG is beyond the ten-fold
classification scheme of defect states proposed by Teo and Kane Teo and Kane
(2010b), in which crystalline symmetries are not considered. We note that the
electric charge localized at the vortex is half-quantized only when chiral
symmetry exists, similar to the cases considered by Teo and Kane Teo and Kane
(2010b). However, due to the quantized $P_{3}$ protected by inversion
symmetry, the vortex always carries a quantized Wannier charge even if the
electric charge accumulated at the vortex is not quantized when chiral
symmetry is broken.

Figure 1: (a), (b) Finite-size structure of uniformly textured graphene with
$2\times 2$ unit cells. Here blue and yellow dots denote two sublattices. When
$\theta=0$, the intercell hopping is bigger than the intracell hopping, and
vice versa when $\theta=\pi$. (c), (d) Parity distribution of occupied states
at TRIMs for $\theta=0$ (c) and $\theta=\pi$ (d). The parity distributions at
three $M$ points are identical. (e), (f) Charge distribution in a finite-size
KTG with $28\times 28$ unit cells. Corner states appear (do not appear) when
$\nu_{\rm 2D}=1$ ($\nu_{\rm 2D}=0$).

##  II Higher-order band topology of KTG

The Hamiltonian for KTG proposed in Ref. Hou et al. (2007) is given by

$\displaystyle H=\sum_{\textbf{r}\in\Lambda_{A}}\sum_{j=1}^{3}\left(t+\delta
t_{j}(\textbf{r})\right)c_{A}^{\dagger}(\textbf{r})c_{B}(\textbf{r}+\textbf{d}_{j})+\textrm{h.c.},$
(1)

where $c_{A,B}(\textbf{r})$ indicates the (spinless) electron annihilation
operator on the sublattice $A$, $B$ at the position r, $\textbf{d}_{1,2,3}$
denotes the vectors connecting an $A$ sublattice site to its neighboring $B$
sublattice sites, and $\Lambda_{A}$ denotes the set of the sites belonging to
the A sublattice. $t$ denotes the hopping amplitude between nearest neighbor
sites, and $\delta t$ describes the modulation of the hopping amplitudes. When
$\delta t=0$, the Hamiltonian describes the low-energy band structure of
graphene having two Dirac points at the Brillouin zone corners with the
momentum $\bm{k}=\pm\bm{K}_{+}$, while $\delta
t_{j}(\textbf{r})=\Delta(\textbf{r})e^{i\textbf{K}_{+}\cdot\textbf{d}_{j}}e^{i\textbf{G}\cdot\textbf{r}}+\mathrm{c.c}$
couples the two Dirac fermions at the opposite valleys, where
$\textbf{G}=2\textbf{K}_{+}$ Hou et al. (2007). When $\Delta(\textbf{r})$ is a
nonzero constant, the unit cell becomes three times larger than that of
pristine graphene, and KTG is formed.

Let us examine the band topology of the gapped graphene with a uniform Kekule
order parameter. To this end, we define a parametrized Hamiltonian
$H_{\theta}$ describing graphene with a uniform Kekule texture
$\Delta(\textbf{r})=\Delta_{1}+i\Delta_{2}=\Delta_{0}e^{i\theta}$, where
$\Delta_{i=0,1,2}$ and $\theta$ are real constant. As shown in Figs. 1(a) and
1(b), we use the convention that the interunit-cell hopping is larger
(smaller) than the intraunit-cell hopping when $\theta=0$ ($\theta=\pi$). When
$\theta=0$ or $\pi$, the system is invariant under inversion $P$, the sixfold
rotation about the $z$ axis $C_{6z}$, two mirrors $M_{x}$ and $M_{y}$
[$M_{x}:(x,y)\rightarrow(-x,y)$, $M_{y}:(x,y)\rightarrow(x,-y)$], and time
reversal $T$ symmetries. On the other hand, when $\theta\neq 0,\pi$, inversion
symmetry is broken, and $H_{\theta}$ and $H_{-\theta}$ are related by
inversion $P$. Also, the chiral (or sublattice) symmetry exists when only the
hopping between nearest-neighbor sites is considered.

In general, a 2D $P$-symmetric spinless fermion system carries three $Z_{2}$
topological invariants, $P_{1x}$, $P_{1y}$, and $\nu_{\rm 2D}$ Hughes et al.
(2011). Here $P_{1a}$ ($a=x,~{}y$) indicates the quantized charge polarization
along the $a$ direction Zak (1989); Xiao et al. (2010), while $\nu_{\rm 2D}$
is a 2D topological invariant characterizing higher-order topological
insulator (HOTI) of $P$-symmetric fermion systems Ahn et al. (2018a); Ahn and
Yang (2018); Ahn et al. (2018b); Wieder and Bernevig (2018). $\nu_{\rm 2D}$
can be determined by using the parity eigenvalues at time-reversal invariant
momenta (TRIM) from the following relation

$\displaystyle(-1)^{\nu_{\rm
2D}}=\prod_{i=1}^{4}(-1)^{[N_{\text{occ}}^{-}(\Gamma_{i})/2]},$ (2)

where $N_{\text{occ}}^{-}(\Gamma_{i})$ is the number of occupied bands with
odd parity at the TRIM $\Gamma_{i}$ and the square bracket $[\alpha]$
indicates the integer part of $\alpha$ Song et al. (2018); Fu and Kane (2007);
Kim et al. (2015); Hughes et al. (2011); Turner et al. (2012); Po et al.
(2017); see Figs. 1(c) and 1(d).

Interestingly, we find that KTG with $\theta=0$ ($\theta=\pi$) is a 2D
inversion symmetric HOTI (a trivial insulator) with $\nu_{\rm 2D}=1$
($\nu_{\rm 2D}=0$) Noh et al. (2018); Benalcazar et al. (2018); Geier et al.
(2018); Khalaf (2018); Ahn et al. (2018a), which is also confirmed by the
Wilson loop spectra shown in the Supplemental Material (SM) SM (2020);
Benalcazar et al. (2017b); Xie et al. (2018). As shown in Fig. 1(e), KTG with
$\theta=0$ exhibits a pair of zero-energy corner states related by $P$
symmetry manifesting its higher-order band topology.

##  III Kekule-textured vortex and an axion insulator

A vortex structure of the order parameter $\Delta(\textbf{r})$ with the
winding number $n$ can be introduced by taking
$\Delta(\textbf{r})=|\Delta(\textbf{r})|e^{in\theta(\textbf{r})}$, where the
polar angle $\theta(\textbf{r})$ varies from $0$ to $2\pi$ encircling the
vortex core. In particular, when $n=1$, one zero-mode wave function $\psi\sim
e^{-\int_{0}^{r}d{r^{\prime}}|\Delta(r^{\prime})|}$ is found, yielding a half
integral charge $e/2$ localized at the vortex core Hou et al. (2007); Chamon
et al. (2008a); Weinberg (1981).

The topological property of a vortex can be understood in terms of the charge
pumping process of $H_{\theta}$ during the adiabatic variation of $\theta$
between $0$ and $2\pi$, which corresponds to the phase change
$\theta(\textbf{r})$ of the order parameter around a vortex. To observe the
corresponding spectral flow, we compute the energy spectrum of the Hamiltonian
$H_{\theta}$ under the open boundary condition along the $x$ and $y$
directions [as in the lattices shown in Figs. 1(a) and 1(b)]. As shown in Fig.
2(a), one eigenstate travels from the valence (conduction) bands to the
conduction (valence) bands during the variation of $\theta$ from $-\pi$ to
$\pi$. Two zero-energy modes at $\theta=0$ account for the corner charges
shown in Fig. 1(e).

The topological origin of the nontrivial spectral flow can be understood from
the symmetry of $H_{\theta}=\sum_{\textbf{k}}H(\textbf{k};\theta)$ as follows.
Here $\textbf{k}=(k_{x},k_{y})$ denotes a momentum in the 2D Brillouin zone.
Under inversion $P$, the parametrized Hamiltonian $H(k_{x},k_{y};\theta)$
transforms as

$\displaystyle PH(k_{x},k_{y};\theta)P^{-1}=H(-k_{x},-k_{y};-\theta).$ (3)

If $\theta$ is taken as a third momentum $k_{z}$, the effective 3D Hamiltonian
$H_{\text{3D}}(\bm{k})\equiv H(k_{x},k_{y};\theta=k_{z})$ can be considered as
a Hamiltonian for a $P$ invariant 3D insulator. Let us note that
$H_{\text{3D}}(\bm{k})$ restricted in the 2D $P$-invariant momentum subspaces
with $k_{z}=0$ and $\pi$ is characterized by the 2D $\mathbb{Z}_{2}$ invariant
$\nu_{\rm 2D}(k_{z})$, as the 2D Hamiltonian $H(k_{x},k_{y};\theta)$ with
$\theta=0$ and $\theta=\pi$ are the $P$ invariant models with a Kekule
texture. For a 3D insulator with $P$ symmetry, it was recently shown Ahn and
Yang (2018); Wieder and Bernevig (2018) that the quantized magnetoelectric
polarizability $P_{3}$ is equivalent to $\Delta\nu_{\rm 2D}=\nu_{\rm
2D}(k_{z}=\pi)-\nu_{\rm 2D}(k_{z}=0)$, namely,

$\displaystyle\Delta\nu_{\rm 2D}=2P_{3}=\frac{1}{4\pi^{2}}\int_{T^{2}\times
S^{1}}\text{Tr}\\!\left[AdA-\frac{2i}{3}A^{3}\right]\mod 2,$ (4)

where $A_{ij}=i\langle u_{i}|du_{j}\rangle$ is the non-Abelian Berry’s
connection characterizing the valence band eigenstates $|u_{i}\rangle$,
$T^{2}$ indicates the 2D Brillouin zone, $S^{1}$ denotes the unit circle
parametrized by $\theta$. Therefore, when $\Delta\nu_{\rm 2D}=1$ (mod 2), the
Hamiltonian $H_{\text{3D}}(\bm{k})$ exhibits quantized $P_{3}=1/2$, that is,
$H_{\mathrm{3D}}(\bm{k})$ describes an axion insulator Qi et al. (2008);
Turner et al. (2012). Due to the nontrivial band topology of
$H_{\text{3D}}(\bm{k})$, the energy spectrum of the Hamiltonian under open
boundary conditions in the $x$ and $y$ directions exhibits a topological
spectral flow as a function of $k_{z}$. Namely, one electron should be pumped
from the valence bands to the conduction bands and vice versa during the
variation of $k_{z}\in[-\pi,\pi]$, which corresponds to two chiral hinge modes
of the axion insulator described by $H_{\mathrm{3D}}(\bm{k})$.

The topological spectral flow of $H_{\mathrm{3D}}(\bm{k})$ explains the charge
accumulation at a vortex as follows. Let us consider a textured lattice with a
vortex-antivortex pair shown in Fig. 2(b), where the parameter $\theta$ is
fixed to be $\theta=\pi$ far away from the vortices. Since $H_{\theta}=\pi$ is
topologically trivial, the region with $\theta=\pi$ can be considered as a
vacuum. Then, in the region with $\theta=0$ between the vortex and the
antivortex, the centers of the vortex and the antivortex can be considered as
two corners of $H_{\theta=0}$ under open boundary condition; See Fig. 1(a).
Therefore, the vortex bound state can be considered as the corner state of the
HOTI.

The spectral flow in Fig. 2(a) implies that an in-gap state travels from the
valence (conduction) to the conduction (valence) bands when $\theta$ changes
from $0$ to $2\pi$ along a circle enclosing a vortex (antivortex). This means
that the in-gap state should cross the Fermi level in the middle of the
variation of $\theta\in[-\pi,\pi]$, leading to a zero-energy bound state at
the core of the vortex (antivortex). This shows that the topological charge
pumping process underlies the relationship between the existence of in-gap
states and the bulk band topology of KTG. We provide a field theoretical
explanation for the connection between topological spectral flow and charge
accumulation at a vortex in SM SM (2020).

Figure 2: (a) Evolution of the energy spectrum of $H_{\theta}$, with open
boundaries in the $x$ and $y$ directions, composed of $28\times 28$ unit
cells. Two zero-energy modes at $\theta=0$ correspond to the corner charges
shown in Fig. 1(e). During one cycle, one state travels from the valence
(conduction) bands to the conduction (valence) bands crossing the Fermi level
at $\theta=0$. (b) Schematic figure describing the lattice structure around a
vortex-antivortex pair in KTG. The blue and red circles encircling the vortex
and the antivortex, respectively, correspond to two chiral modes (the blue and
red lines) in (a). Far away from the vortex-antivortex pair, the lattice has a
uniform texture with $\theta=\pi$.

##  IV Generalization

A similar mechanism for charge accumulation at a vortex in a textured lattice
due to the topological spectral flow can be applied to any parametrized 2D
Hamiltonian $H(k_{x},k_{y};\theta)$ when its corresponding 3D Hamiltonian
$H_{\text{3D}}(\bm{k})$ exhibits quantized $P_{3}$. In the following
discussion, the presence or absence of spin-orbit coupling makes no difference
unless noted otherwise. Let us note that $P_{3}$ is quantized when the 3D
insulator described by $H_{\text{3D}}(\bm{k})$ has a space-time orientation
reversing symmetry such as $\mathcal{T}$, $\mathcal{P}$,
$\mathcal{C}_{n}\mathcal{T}$, and $\mathcal{C}_{n}\mathcal{P}$ ($n=2,3,4,6$)
Varjas et al. (2015); Qi et al. (2008); Schindler et al. (2018); Ahn and Yang
(2018); Ryu et al. (2010), where we have used calligraphic fonts for
symmetries in the 3D space to distinguish them from the symmetries in the 2D
space denoted in italic fonts. Among these symmetries, we can neglect
$\mathcal{C}_{3z}\mathcal{T}$ in that a
$\mathcal{C}_{3z}\mathcal{T}$-symmetric axion insulator is expected to carry
the same topological properties as the one with $\mathcal{T}$ symmetry only,
because $(\mathcal{C}_{3z}\mathcal{T})^{3}=\mathcal{T}$ and
$(\mathcal{C}_{3z}\mathcal{T})^{2}=\mathcal{C}_{3z}^{-1}$, and
$\mathcal{C}_{3z}$ by itself cannot quantize $P_{3}$. Likewise, we do not have
to consider $\mathcal{C}_{6z}\mathcal{T}$ since
$(\mathcal{C}_{6z}\mathcal{T})^{3}=\mathcal{C}_{2z}\mathcal{T}$. For
$\mathcal{C}_{nz}\mathcal{P}$ symmetries,
$(\mathcal{C}_{3z}\mathcal{P})^{3}=\mathcal{P}$ and
$(\mathcal{C}_{6z}\mathcal{P})^{3}=\mathcal{C}_{2z}\mathcal{P}=\mathcal{M}_{z}$.
Thus, we only have to consider $H_{\text{3D}}(\bm{k})$ with $\mathcal{T}$,
$\mathcal{P}$, $\mathcal{C}_{2z}\mathcal{P}$, $\mathcal{C}_{2z}\mathcal{T}$,
$\mathcal{C}_{4z}\mathcal{P}$, or $\mathcal{C}_{4z}\mathcal{T}$ symmetries. It
is worth noting that in each axion insulator described by
$H_{\text{3D}}(\bm{k})$, the 2D subspaces with $k_{z}=0$ and $k_{z}=\pi$,
respectively, support distinct 2D topological invariants. Namely, the
quantized $P_{3}$ of $H_{\text{3D}}(\bm{k})$ can be obtained from the
difference between the 2D topological invariants on the two symmetry invariant
planes with $k_{z}=0$ and $k_{z}=\pi$, respectively Fu and Kane (2006);
Bernevig et al. (2006); Teo et al. (2008); Van Miert and Ortix (2018);
Schindler et al. (2018). This ensures that the corresponding vortex structure
described by $H(k_{x},k_{y},\theta)$, which connects two distinct topological
phases at $\theta=0$ and $\pi$, hosts a fractional vortex charge.

Table 1:  The correspondence between the symmetry of $H_{\text{3D}}(\bm{k})$ exhibiting quantized $P_{3}$ and that of the 2D Hamiltonian $H(k_{x},k_{y};\theta=0,\pi)$. For $T$ and $C_{4}T$ symmetry, spin orbit coupling (SOC) is necessary for having $P_{3}=1/2$, while SOC is not required for other symmetries. (See the SM). The nature of the corresponding 2D topological insulator (TI) is shown in the fourth column. Here “QSHI” indicates a quantum spin Hall insulator. 3D Symmetry | 2D Symmetry | SOC | 2D TI
---|---|---|---
$\mathcal{P}=\mathcal{C}_{2z}\mathcal{M}_{z}$ |  $P$, $C_{2z}$, $C_{6z}$ |  | HOTI
$\mathcal{C}_{4z}\mathcal{P}=(\mathcal{C}_{4z}\mathcal{M}_{z})^{3}$ | $C_{4z}$ |  | HOTI
$\mathcal{T}$ | $T$ | ✓ | QSHI
$\mathcal{M}_{z}$ | $M_{z}$ |  | Mirror TI
$\mathcal{C}_{2z}\mathcal{T}$ |  $PT$, $C_{2z}T$, $C_{6z}T$ |  | HOTI
$\mathcal{C}_{4z}\mathcal{T}$ | $C_{4z}T$ | ✓ | HOTI

Now let us explain how the symmetry of the 3D Hamiltonian
$H_{\text{3D}}(\bm{k})$ can be related to that of the physical 2D parametrized
Hamiltonian $H(k_{x},k_{y};\theta)$. For instance, as in the case of KTG, if
the system with a uniform order parameter is invariant under $n$-fold rotation
$C_{n}$ about the $z$ axis at $\theta=0$ and $\pi$, but not at other $\theta$
values, then the $C_{n}$ symmetry can be implemented as

$\displaystyle
C_{n}H(k_{x},k_{y};\theta)C_{n}^{-1}=H(k^{\prime}_{x},k^{\prime}_{y};-\theta),$
(5)

where $k^{\prime}_{x}$, $k^{\prime}_{y}$ are the rotated momenta after $C_{n}$
operation. Due to the sign change of $\theta$ under $C_{n}$, the 2D system is
$C_{n}$ invariant only at $\theta=0,~{}\pi$. Then the $C_{n}$ symmetry of the
2D system can be implemented in $H_{\text{3D}}(\bm{k})$ as the
$\mathcal{C}_{n}\mathcal{M}_{z}$ symmetry. Now we ask whether
$\mathcal{C}_{n}\mathcal{M}_{z}$ symmetry can quantize $P_{3}$ and also
whether $C_{n}$ symmetry supports a 2D topological invariant $\nu_{\rm 2D}$
such that its difference $\Delta\nu_{\rm 2D}=\nu_{\rm 2D}(\theta=\pi)-\nu_{\rm
2D}(\theta=0)$ is identical to $2P_{3}$ (modulo two). Similar ideas can be
applied to all the other symmetries that can give quantized $P_{3}$ SM (2020).

Table I summarizes the correspondence between the symmetry of the 2D insulator
with a uniform gap-opening order parameter like the Kekule texture, described
by the Hamiltonian $H(k_{x},k_{y};\theta=0,\pi)$, and the symmetry of the
relevant 3D Hamiltonian $H_{\text{3D}}(\bm{k})$ exhibiting quantized $P_{3}$.
When the two Hamiltonian $H(k_{x},k_{y};\theta=0)$ and
$H(k_{x},k_{y};\theta=\pi)$ support distinct bulk topological properties, a
fractional charge can be localized at a vortex of the order parameter
parametrized by $\theta$ due to the topological spectral flow associated with
the quantized $P_{3}$.

##  V Chiral symmetry and charge quantization

In fact, the electric charge accumulated at a vortex is half-quantized only
when chiral symmetry exists. When chiral symmetry is broken, the midgap states
can be shifted from the zero energy and the vortex charge takes an arbitrary
value. Yet, the topological invariant $P_{3}=\nu_{\rm 2D}(\theta=\pi)-\nu_{\rm
2D}(\theta=0)$ is quantized irrespective of the presence or absence of chiral
symmetry, and this guarantees a quantized half-integral Wannier charge
localized at a vortex as long as the relevant crystalline symmetry is
preserved.

Figure 3: (a) The location of Wannier centers when $\theta=0$ corresponding to
the 2D HOTI with $\nu_{\rm 2D}=1$. (b) A similar figure when $\theta=\pi$
corresponding to a trivial insulator with $\nu_{\rm 2D}=0$. (c) Distribution
of Wannier centers around a vortex. The Wannier functions centered at the unit
cell boundary when $\theta=0$ contributes the quantized half-integer Wannier
charge in the shaded region including the vortex core.

To define a Wannier charge, we assume that an electron is sharply localized at
the Wannier center for the counting of electron numbers, and the Wannier
centers are located at the Wyckoff positions that are invariant under the
crystalline symmetres up to a lattice constant vector. Then the electron
number sharply localized at a Wannier center is the Wannier charge.

The quantized Wannier charge can be seen from the distribution of Wannier
center positions described in Fig. 3. When $\theta=0$ ($\theta=\pi$) with $P$
symmetry, the locations of Wannier centers are fixed at the boundary (center)
of the unit cell. Thus, these two states cannot be adiabatically connected
when $P$ exists. On the other hand, for other $\theta$ with broken $P$, the
Wannier centers are located at generic points within a unit cell, and always
contribute an integer number of electrons. Thus, a group of connected unit
cells including the vortex core, such as the shaded region in Fig. 3 (c),
always contains a half-quantized Wannier charge due to the Wannier centers at
the unit cell boundary when $\theta=0$. This is similar to the situation of
the $P$ symmetric SSH chain having quantized charge polarization, which gives
a half-integral Wannier charge at the domain wall irrespective of chiral
symmetry SM (2020).

This Wannier charge description can generally be applied to the cases with
other symmetries SM (2020). Let us note that the topological properties of the
2D HOTI lacking chiral symmetry, especially its relation with the third-order
TI and the related filling anomly, have been extensively studied in recent
papers Călugăru et al. (2019); Okuma et al. (2019); Song et al. (2017); Lee et
al. (2020); Hwang et al. (2019).

##  VI Discussion

In real materials, the $U(1)$ symmetry of the order parameter may be reduced
to a discrete $\mathbb{Z}_{n=2,3,4,6}$ symmetry by lattice potentials.
However, fractional charges can still be localized at a vortex with discrete
symmetry as long as the phase $\theta$ of the order parameter winds $2\pi$
around the vortex. For instance, in the case of KTG, fractional charge can be
localized at a junction where three domains with $\theta=0$, $\pi/3$,
$5\pi/3$, respectively, meet. Alternatively, one can construct
heterostructures composed of 2D topological and normal insulators to observe
fractional bound charges SM (2020); Qi et al. (2008b).

One interesting direction for future research is to extend the idea of
symmetry protected topological vortices to defects of various codimensions.
For example, it is shown in Ref. Teo and Kane (2010b) that a point defect
carrying fractional charges in 3D systems can be described by using the
Hamiltonian for five-dimensional insulators whose topological invariant has
the Chern-Simons $5$-form. Systematic classification of such defect structures
taking into account space group symmetries is desirable to complete the
classification table for defect Hamiltonian beyond the tenfold classification
scheme proposed before Ref. Teo and Kane (2010b).

###### Acknowledgements.

E.L. was supported by IBS-R009-D1. A.F. was supported by JSPS KAKENHI (Grant
No. 15K05141 and 19K03680) and JST CREST (Grant No. JPMJCR19T2). B.-J.Y. was
supported by the Institute for Basic Science in Korea (Grant No. IBS-R009-D1)
and Basic Science Research Program through the National Research Foundation of
Korea (NRF) (Grant No. 0426-20180011), and the POSCO Science Fellowship of
POSCO TJ Park Foundation (No. 0426-20180002). This work was supported in part
by the U.S. Army Research Office under Grant No. W911NF-18-1-0137. We
appreciate the helpful discussions with J. Ahn and Y. Hwang.

###### Contents

  1. I Introduction
  2. II Higher-order band topology of KTG
  3. III Kekule-textured vortex and an axion insulator
  4. IV Generalization
  5. V Chiral symmetry and charge quantization
  6. VI Discussion
  7. S1 Two different ways of describing charge fractionalization in SSH model
  8. S2 Nested Wilson loop
  9. S3 Influence of core electrons to the higher order band topology
  10. S4 Second- and third-order topological insulators
  11. S5 Classification of Axion insulators by crystalline symmetries.
  12. S6 Field-theoretical description for a topological vortex
  13. S7 Wannier center and quantized Wannier charge in topological vortex

##  S1 Two different ways of describing charge fractionalization in SSH model

Here let us review the physics of charge fractionalization described by the
Su-Schriefer-Heeger (SSH) model Su et al. (1980). The Hamiltonian can be
written as

$\displaystyle H_{\textrm{SSH}}=\sum(t+\delta
t_{i})a^{\dagger}_{i}a_{i+1}+H.c,$ (S1)

where $t$ is the hopping amplitude and its modulation is described by $\delta
t_{i}=(-1)^{i}\Delta$. Also, we assume that the system is half-filled. If
$\Delta$ is 0, the energy dispersion is gapless and the unit cell contains one
atom. However, when $\Delta$ is non-zero constant, the unit cell becomes
doubled so that the size of the Brillouin zone is reduced by half and gapped
band structure appears. Spontaneous formation of lattice dimerization gives
rise to two degenerate ground states distinguished by the opposite sign of the
mass term $\Delta$. Depending on the sign of $\Delta$, the system becomes
either trivial or topological. In our convention, $\Delta>0$ corresponds to
the topological phase and $\Delta<0$ corresponds to the trivial phase. At the
boundary between topological and trivial phases, there exists a localized
zero-mode wave function and the bound charge is quantized to 1/2. Let us
introduce two different ways describing charge fractionalization at the domain
wall.

In the first approach, we solve the low energy Dirac Hamiltonian given by

$H(x)=\left(\begin{array}[]{cc}-i\partial_{x}&\Delta(x)\\\
\Delta^{*}(x)&i\partial_{x}\end{array}\right).$ (S2)

Across the domain wall, $\Delta(x)$ changes from $-\Delta_{0}$ to
$\Delta_{0}$, and it is complex in between. The equation satisfied by the
zero-mode solution is $H(x)\psi(x)=0$ where $\psi(x)=(u(x),v(x))^{T}$. More
explicitly,

$\displaystyle-i\partial_{x}u(x)+\Delta(x)v(x)=0,$
$\displaystyle\Delta^{*}(x)u(x)+i\partial_{x}v(x)=0.$ (S3)

Under the solitonic mass background, one can find a normalizable zero energy
wave function $\psi(x)$ $\sim
e^{-\int_{0}^{x}|\Delta(x^{\prime})|dx^{\prime}}$, which gives a half electric
charge localized at the domain wall.

In the second approach, one can explain the existence of fractional charge in
terms of the bulk topological property. Let us consider a Hamiltonian
$H(k_{x},\theta)=\sin{k_{x}}\sigma_{z}+\Delta_{x}(k_{x},\theta)\cos\theta\sigma_{x}+\Delta_{y}(k_{x},\theta)\sin\theta\sigma_{y}$
that varies adiabatically along a circle parametrized by $\theta\in[0,2\pi]$
without closing a bulk gap, as in Ref. Teo and Kane, 2010b (See Fig. S1(a)).
We assume that Hamiltonian describes a trivial (topological) phase when
$\theta=0$ $(\pi)$ with the corresponding polarization $P_{1}=0$
($P_{1}=1/2$), where $P_{1}=\frac{1}{2\pi}\int_{T^{1}}\text{Tr}[A]$ is
quantized to be either $0$ or $1/2$ $\mod 1$ due to the inversion $P$
symmetry. On the other hand, $P_{1}$ is not quantized for the states with
$\theta\neq 0,~{}\pi$, since $P$ is broken due to the term proportional to
$\sigma_{y}$.

In the presence of $P$ symmetry, $P_{1}$ can also be computed by multiplying
the parity eigenvalues at the time reversal invariant momenta (TRIM) Hughes et
al. (2011); Turner et al. (2012) as

$\displaystyle(-1)^{2P_{1}}=(-1)^{N_{\text{occ}}^{-}(k=0)}(-1)^{N_{\text{occ}}^{-}(k=\pi)},$
(S4)

where $N_{\text{occ}}^{-}(\Gamma_{i})$ is the number of occupied states with
negative parity at $\Gamma_{i}$.

Figure S1: (Color online) (a) An adiabatic cycle interpolating two insulators
with the charge polarization $P_{1}=0$ and $P_{1}=1/2$, respectively. (b)
Schematic figures describing the trivial (topological) phase with the
parameter $\theta=0$ ($\theta=\pi$)

The bound charge between two insulators with $P_{1}=0$ ($\theta=0$) and
$P_{1}=1/2$ ($\theta=\pi$) is expressed by the integral of the polarization
gradient

$\displaystyle q_{\text{bound}}=\int-\partial_{\theta}P_{1}d\theta=1/2\textrm{
mod 1},$ (S5)

which indicates the half integral electric charge localized between two
phases.

While the parametrized Hamiltonian $H(k_{x},\theta)$ is $P$ invariant only at
$\theta=0$ and $\pi$, $P$ symmetry can be implemented as

$\displaystyle PH(k_{x},\theta)P^{-1}=H(-k_{x},-\theta).$ (S6)

Interestingly, when $\theta$ is replaced by $k_{y}$, the Hamiltonian
$H(k_{x},k_{y})$ describes a 2D Chern insulator with $P_{2D}$ symmetry,

$\displaystyle P_{2D}H(k_{x},k_{y})P_{2D}^{-1}=H(-k_{x},-k_{y}).$ (S7)

This is because the Chern number C can be related to the change of the charge
polarization through the parity formula for occupied bands as

$\displaystyle(-1)^{C}=(-1)^{\sum
N_{\text{occ}}^{-}(\Gamma_{i})}=(-1)^{2P_{1}(\theta=\pi)-2P_{1}(\theta=0)},$
(S8)

where $\Delta(2P_{1})=C\mod 2$. Thus, $\Delta(2P_{1})=1$ implies that there is
a nontrivial spectral flow from the conduction band to the valence band during
the adiabatic evolution of $\theta$ from $0$ to $2\pi$, which corresponds to
the chiral edge mode of the Chern insulator Qi et al. (2008). Since the chiral
edge mode connects the valence and conduction bands, even if chiral symmetry
is broken, the half charge should be localized at the domain wall between the
topological and trivial phases in $P$ symmetric SSH model.

As shown in the main text, the fractional charges bound to the vortex of the
order parameter can also be explained by the bulk topological response of the
effective Hamiltonian for axion insulators in three dimensions. Table I
compares the charge fractionalization at the domain wall of the SSH model and
at the order parameter vortex of the HCM model.

| Domain wall of SSH model | Vortex of HCM model
---|---|---
$H_{\rm D+1}$ | Chern insulator | Axion insulator
Invariant | $C_{1}=\Delta(2P_{1})$ | $2P_{3}=\Delta\nu_{2D}$
Table 2: Comparison of the charge fractionalization mechanism at a domain wall
in 1D systems described by the SSH model and that at an order parameter vortex
in 2D systems described by the HCM model. Fractional charges localized at a
point defect in $D$-dimensional systems can be explained by the topological
response of the $D+1$ dimensional effective Hamiltonian $H_{D+1}$.

##  S2 Nested Wilson loop

A Wilson loop is a gauge invariant observable whose eigenvalue spectrum
contains the information on the topological properties of the Hamiltonian. In
$2D$, a Wilson loop operator is defined by

$\displaystyle W_{(k_{1}+2\pi,k_{2})\leftarrow(k_{1},k_{2})}$
$\displaystyle\equiv W_{x,\textbf{k}}$
$\displaystyle=\lim_{N\rightarrow\infty}F_{N-1}F_{N-2}\dotsb F_{1}F_{0}$
$\displaystyle=Pe^{-i\oint_{C}A_{k}dk},$ (S9)

where $[F_{i}]_{nm}=\langle
u_{m}(k_{i+1},k_{2})|u_{n}(k_{i},k_{2})\rangle,k_{i}=\frac{2\pi}{N}i,$ and
$m,n=1,\dotsb,N_{\rm occ},i=1,\dotsb,N$. Since a Wilson loop operator is
unitary, the eigenvalue equation is given by
$W_{x,\textbf{k}}|\nu_{x,\textbf{k}}^{\ell}\rangle=e^{i\nu_{x}^{\ell}(k_{y})}|\nu_{x,\textbf{k}}^{\ell}\rangle$,
where $\nu_{x}^{\ell}(k_{y})$ corresponds to the $x$-component of the Wannier
center of $\ell$th Wannier functions. It follows that the electron charge
polarization is expressed as

$\displaystyle p_{x}$ $\displaystyle=\frac{1}{2\pi
N_{y}}\sum_{k_{y}}\sum_{\ell=1}^{N_{occ}}\nu_{x}^{\ell}(k_{y})$
$\displaystyle=-\frac{i}{2\pi
N_{y}}\sum_{k_{y}}\log\det[W_{x,\textbf{k}}]=-\frac{1}{2\pi}\oint\text{Tr}[A_{x,k}]d^{2}k.$
(S10)

In the presence of $P$ symmetry, the set of eigenvalues satisfy
$\\{\nu_{x}(k_{y})\\}\equiv\\{-\nu_{x}(-k_{y})\\}$ mod $2\pi$, so that the
polarization is quantized to either $0$ or $1/2$ modulo 1.

Recently, the nested Wilson loop method was developed to study higher order
topological properties benalcazar2017quantized; Benalcazar et al. (2017b). The
procedure for computing the nested Wilson loop is as follows. First, we
calculate a Wilson loop operator along a reciprocal vector $\bf{G_{1}}$, where
its eigenvalues are $\\{e^{i\nu_{1}(k_{2})}\\}$. Next, we choose a certain
subset of Wilson loop eigenvalues $\\{\nu_{1}\\}$ and find the corresponding
Wannier eigenfunctions$\\{|\nu_{1}(k_{2})\rangle\\}$. With the eigenfunctions,
we calculate a Wilson loop operator along the other reciprocal vector
$\bf{G_{2}}$, where we obtain a nested Wilson loop operator $\widetilde{W}$ as

$\displaystyle\widetilde{W}_{k_{2}+2\pi\leftarrow
k_{2}}=\widetilde{W}_{2}=\lim_{N\rightarrow\infty}\widetilde{F}_{N-1}\widetilde{F}_{N-2}\dotsb\widetilde{F}_{1}\widetilde{F}_{0},$
(S11)

where
$[\widetilde{F}_{i}]_{nm}=\langle\nu_{m}(k_{i+1})|\nu_{n}(k_{i})\rangle,k_{i}=\frac{2\pi}{N}i$
and $m,n=1,\dotsb,N_{\text{sub}}$. Here $N_{\text{sub}}$ denotes the number of
the subbands of a gapped Wilson loop spectrum. By using this method, it is
possible to detect the electric multipole moments of the system
benalcazar2017quantized; Benalcazar et al. (2017b); Xie et al. (2018).

Figure S2: (Color online) (a) The determinant of nested Wilson loops for the
model Hamiltonian describing the Kekule textured graphene, which is plotted
for $0\leq\theta\leq\pi$. At $\theta=\theta_{c}$, $\det(\widetilde{W})$
changes abruptly. (b) The Wilson loop spectrum for $\theta=0$. The spectrum
has a Dirac-like crossing point at $k=0$, and the determinant of the nested
Wilson loop is $-1$. (c) The Wilson loop spectrum for $\theta=\pi$. The Wilson
loop spectrum of a trivial insulator has no band crossing and the determinant
of the nested Wilson loop is $1$. (d) The Wilson loop spectrum for
$\theta=\theta_{c}$, where the gap is closed, which implies the phase
transition between $\det(\widetilde{W})=-1$ and $+1$.

The Kekule textured graphane has inversion $P$ symmetry, where the Wilson loop
spectrum satisfies $\\{\nu_{1}(k_{2})\\}\equiv\\{-\nu_{1}(-k_{2})\\}\mod
2\pi$. Also, it has spinless $T$ symmetry, where the Wilson loop spectrum
satisfies $\\{\nu_{1}(k_{2})\\}\equiv\\{\nu_{1}(-k_{2})\\}\mod 2\pi$. Put
together, $\\{\nu_{1}(k_{2})\\}\equiv\\{-\nu_{1}(k_{2})\\}\mod 2\pi$. Let us
note that $T$ is not crucial in quantizing $\det{\widetilde{W}}$, but we have
introduced the symmetry for the sake of simplicity. Thus, a Wilson band must
be $\nu_{1}(k_{2})=0$, $\pi$ or exist as a pair
$\\{+\nu_{1}(k_{2}),-\nu_{1}(k_{2})\\}$, so that the spectrum can be divided
into two subsets that are centered at either $\nu_{1}=0$ or $\nu_{1}=\pi$ Ahn
et al. (2018b). We choose two symmetric Wannier bands that are above and below
$\nu_{1}=\pi$, and calculate the relevant nested Wilson loop $\widetilde{W}$.
The associated Wannier sector polarization is
$p_{y}^{\nu_{1}}=-\frac{i}{2\pi}\log\det[\widetilde{W}]$. Under $PT$ symmetry,
$p_{2}^{\nu_{1}}$ is quantized to $1/2$ or $0$ since $p_{2}^{\nu_{1}}\equiv-
p_{2}^{-\nu_{1}}$ mod 1 benalcazar2017quantized; Benalcazar et al. (2017b).
Also, it is known that $\det{\widetilde{W}}$ is equal to $(-1)^{\nu_{2D}}$ Ahn
et al. (2018b). Thus, $\nu_{2D}=2p_{2}\mod 1$.

For the Kekule textured graphene, $P$ exists when $\theta=0$ or $\pi$. In Fig.
S2, the Wilson loop spectrum of Kekule textured graphene is shown, where
$\det(\widetilde{W})$ is quantized as $-1$ for $\theta=0$ while
$\det(\widetilde{W})$ is quantized as $1$ for $\theta=0$. As $\theta$ varies,
we observe that the Wilson loop spectrum undergoes a phase transition between
$\theta=0$ and $\theta=\pi$. $\det(\widetilde{W})$ changes abruptly and the
gap of the Wilson spectrum is closed at a critical point.

##  S3 Influence of core electrons to the higher order band topology

Figure S3: (Color online) (a) $sp^{2}$ orbitals on the honeycomb lattice. When
only nearest neighbor hopping is considered, electrons are localized at
$\sigma$ bonds. (b) The Wilson loop spectrum for core levels composed of
$sp^{2}$ orbital exhibiting $\nu_{2D}=1$. There are nine Wilson bands. Seven
of them are at $\nu=0$ and two of them are at $\nu=\pi$.

To construct the tight-binding Hamlitonian for the Kekule textured graphene,
only $p_{z}$ orbitals are taken into account since they are mainly responsible
for the low energy electronic property near the Fermi level. Core energy
levels that are far away from the Fermi level are assumed to have negligible
effect. However, when the bulk topological invariant is concerned, core levels
far below the Fermi level can also play a crucial role. For instance, in
graphene, $sp^{2}$ orbitals form core levels below the Fermi level since the
$\sigma$ bonding between $sp^{2}$ orbitals is stronger than the $\pi$ bonding
between $p_{z}$ orbitals. To see the effect of core energy level in higher
order topological physics, let us construct a simple tight-binding Hamiltonian
for $sp^{2}$ orbitals, which only takes into account nearest neighbor hopping,

$\displaystyle H=\sum_{n=1}^{9}t_{\sigma}a_{2i-1}^{\dagger}a_{2i}+H.c,$ (S12)

where $t_{\sigma}$ is hopping parameter for a $\sigma$ bonding and $a_{i}$ is
the annihilation operator of the $i$-th $sp^{2}$ orbital. There are 9 $sp^{2}$
orbitals filled in the unit cell as shown in Fig. S3(a). The energy
eigenvalues are degenerate and the electronic bands are dispersionless because
electrons are localized between neighboring two atoms. From the Wilson loop
spectrum or the parity eigenvalues at TRIM points, the $2D$ topological
invariant $\nu_{2D}$ turns out to be nontrivial. Thus, in Kekule textured
graphene, when the lower energy levels are included, the Kekule texture with
$\theta=0$ becomes trivial ($\nu_{2D}=0$) while that with $\theta=\pi$ becomes
nontrivial ($\nu_{2D}=1$) (See Fig. S3(b)). However, the change of $\nu_{2D}$
remains the same, and thus fractional charges bound to vortices are also
robust against the adding additional bands below the Fermi level.

##  S4 Second- and third-order topological insulators

The zero-mode corner states of a higher-order topological insulator (HOTI) in
2D systems are protected only if chiral symmetry exists. When chiral symmetry
is broken, the energy of the corner states can be shifted from the zero
energy, and even merged into the bulk spectrum. However, since the nontrivial
value of $\nu_{\rm 2D}$ is determined only by the parity eigenvalues, the
insulator without chiral symmetry cannot be adiabatically connected to a
trivial insulator, which implies that the insulator without chiral symmetry
still possesses nontrivial topology that is distinct from a trivial insulator.
In fact, the 2D topological system with/without chiral symmetry is called the
second-order/third-order topological insulator (SOTI/TOTI) Călugăru et al.
(2019); Okuma et al. (2019).

The nontrivial band topology of the inversion symmetric HOTIs is manifested in
filling anomaly, which is described as follows Song et al. (2017);
benalcazar2019quantization. In the chiral symmetric system, the number of
occupied and unoccupied states are the same: $N_{\rm occ}=N_{\rm unocc}$. For
the inversion symmetric SOTI, there are two zero-mode corner states connected
by the inversion symmetry. However, the global geometry of the SOTI cannot
satisfy the inversion symmetry when the system is half-filled, since only one
of the corner states should be occupied (see Fig. 1(e), for example). This
phenomenon is called filling anomaly, which still holds even if chiral
symmetry is broken, since the energy of the in-gapped states connected by the
inversion symmetry are the same Song et al. (2017);
benalcazar2019quantization.

The low energy properties of 2D TOTI and its relation with 2D SOTI are given
as follows. The low energy Dirac Hamiltonian for the boundary states of a 2D
first-order TI can be written as
$H_{edge}=-iR^{-1}\tilde{\sigma}_{z}\partial_{\theta}$, which describes the
gapless chiral edge mode at the boundary. We use $\tilde{\sigma}$ to describe
the symmetry representation for the surface states. For simplicity, we assume
a disk-shaped geometry of the system with the radius $R$ and the corresponding
polar angle $\theta$. For the inversion symmetric SOTIs, on the other hand,
one can introduce mass term $m_{1}(\theta)\tilde{\sigma}_{x}$ which satisfies
the constraint $m_{1}(\theta+\pi)=-m_{1}(\theta)$ where inversion
$P=\tilde{\sigma}_{z}$. Thus, the edge Hamiltonian becomes

$\displaystyle
H_{edge}=-iR^{-1}\tilde{\sigma}_{z}\partial_{\theta}+\tilde{m}_{1}(\theta)\tilde{\sigma}_{x},$
(S13)

which is obviously chiral symmetric. Due to the constraint above, the mass gap
is closed at two points at least, which indicates two corner states related by
the inversion symmetry.

In the case of a 2D TOTI, chiral symmetry is broken so that two mass terms
$m_{1}(\theta)\tilde{\sigma}_{x}$ and $m_{2}(\theta)\tilde{\sigma}_{y}$ are
allowed, where $m_{1,2}(\theta+\pi)=-m_{1,2}(\theta)$. Thus, the surface
spectrum is gapped in general. Nevertheless, it is still topological in the
sense that mass winding number $w$ defined below is non-zero.

$\displaystyle w=\int
d\theta\frac{1}{2\pi}\partial_{\theta}\tan^{-1}\left(\frac{m_{2}(\theta)}{m_{1}(\theta)}\right)$
(S14)

This non-zero winding number is proportional to the charge accumulation at the
boundary, which directly reflects the filling anomaly Wieder and Bernevig
(2018); Hwang et al. (2019).

##  S5 Classification of Axion insulators by crystalline symmetries.

In the main text, we mentioned that space-time orientation reversing
symmetries quantize magnetoelectric polarizability $P_{3}$, and
$H_{\text{3D}}(\bm{k})$ describes a 3D axion insulator. The relevant
symmetries of $H(k_{x},k_{y},\theta=0,~{}\pi)$ are $P$, $T$, $M_{z}$,
$C_{2z}T$, $C_{4z}P$ and $C_{4z}T$. The quantized $P_{3}$ of each 3D axion
insulator can be described by a pumping process of 2D invariants that are
defined at two different symmetry invariant planes with $k_{z}=0$ and $\pi$.
Among the symmetries, $T$ and $M_{z}$ symmetric axion insulators are described
by a pumping process of topological phases with the first order band topology
(QSHI and mirror Chern insulators) and trivial phases. For the other
symmetries, the fractional charge at the order parameter vortex is described
by a pumping process of a higher order topological insulator (HOTI) and a
trivial insulator. Here, we explicitly describe the low energy Hamiltonian and
its physical property for each symmetry class.

Figure S4: (Color online) (a) A heterostructure where fractional charges are
bound at the junction between a QSHI ($\theta=0$) and two ferromagnetic
insulators (FM, $\theta=\pm\alpha$) with opposite in-plane magnetization.
Since two FMs introduce surface mass terms with opposite sign, a zero-mode
state is localized at the junction. (b) Similar heterostructure for a mirror
Chern insulator. (c) A $\mathbb{Z}_{3}$ vortex of Kekule textured graphene
where three domains with lattice structure corresponding to $\theta=0$,
$\pi/3$, $-\pi/3$ meet at a junction. Due to the Wannier function centered at
the unit cell boundary at $\theta=0$, the vortex geometry contains half-
integer electrons. (d) A $\mathbb{Z}_{4}$ vortex relevant to $C_{4z}$
symmetric insulators with a quantized quadrupole moment where four domains
with lattice structure corresponding to $\theta=0$, $\pi/2$, $\pi$, and
$3\pi/2$ meet at a junction. Due to the Wannier functions centered at the unit
cell corner at $\theta=0$, half-integer electrons are bound to the vortex
center.

$\mathcal{T}$ symmetry.— 3D time reversal invariant $Z_{2}$ TI is described by
a pumping process between a QSHI and a trivial insulator defined at the
$k_{z}=0$ and $\pi$ planes Fu and Kane (2006). Let us note that SOC is crucial
since 2D AI class is trivial according to the tenfold classification Ryu et
al. (2010). The low energy Hamiltonian for a 2D QSHI is

$\displaystyle
H(\textbf{k})=k_{x}\tau_{x}\sigma_{x}+k_{y}\tau_{x}\sigma_{y}+M\tau_{z}\sigma_{0},$
(S15)

where Pauli matrices $\tau$ and $\sigma$ represent orbital and spin degrees of
freedom and time reversal symmetry is $T=\sigma_{y}K$ satisfying
$TH(\textbf{k})T^{-1}=H(-\textbf{k})$. The corresponding gapped insulator
should be either a QSHI ($M>0$) or a trivial insulator ($M<0$) depending on
the sign of the mass $M$ Bernevig et al. (2006). The Hamiltonian with a vortex
can be expressed as

$\displaystyle H(k_{x},k_{y},\theta)$
$\displaystyle=k_{x}\tau_{x}\sigma_{x}+k_{y}\tau_{x}\sigma_{y}$
$\displaystyle+\Delta\cos\theta\tau_{z}\sigma_{0}+\Delta\sin\theta\tau_{x}\sigma_{z},$
(S16)

where $\tau_{x}\sigma_{z}$ anticommutes with other matrices and it is odd
under time reversal symmetry, so that the vortex structure breaks time-
reversal symmetry except at $\theta=0$ (a trivial insulator) and at
$\theta=\pi$ (a QSHI).

As shown in Fig. S4 (a), a vortex structure can also be realized at the
junction where a QSHI ($\theta=\pi$) and two ferromagnetic insulators
($\theta=\pm\alpha$) meet. The Hamiltonian describing the gapless boundary of
the QSHI is

$\displaystyle H_{0}$ $\displaystyle=v_{F}\int
dx(\psi_{R}^{\dagger}i\partial_{x}\psi_{R}-\psi_{L}^{\dagger}i\partial_{x}\psi_{L}),$
$\displaystyle=v_{F}\int dx\Psi^{\dagger}i\tilde{\sigma}_{z}\partial_{x}\Psi,$
(S17)

where $\psi_{R/L}$ is the right/left moving states and
$\Psi^{\dagger}=(\psi_{R}^{\dagger},\psi_{L}^{\dagger}).$ The ferromagnetic
insulators introduce time-reversal breaking mass terms that open the band gap
of the edge state. The resulting the surface Hamiltonian is $H=H_{0}+\int
dx\Psi^{\dagger}(m_{1}\tilde{\sigma}_{x}+m_{2}\tilde{\sigma}_{y})\Psi$.
Parametrizing $m_{1}=m\cos\phi$, $m_{2}=m\sin\phi$, the charge density at the
magnetic domain is expressed as $\rho=\frac{1}{2\pi}\partial_{x}\phi$
Goldstone and Wilczek (1981). When the spin direction is reversed, the sign of
the mass terms is also reversed ($\phi\rightarrow\phi+\pi$). Thus, the mass
gap is closed at the junction, indicating a zero-mode state and fractional
charges bound at the vortex center $Q=\int\rho dx=1/2$ Qi et al. (2008b).

$\mathcal{M}_{z}$ symmetry.— For $\mathcal{M}_{z}$ symmetry, a mirror Chern
number $C_{m}$ can be defined on the $k_{z}=0$ and $k_{z}=\pi$ planes. $C_{m}$
is given by the difference of the Chern number $C$ of two sectors with
different mirror eigenvalues. The difference of the mirror Chern numbers on
these planes corresponds to $2P_{3}$, that is,
$C_{m}(k_{z}=0)-C_{m}(k_{z}=\pi)\equiv 2P_{3}\mod 2$ Ahn and Yang (2018). Let
us consider the low energy Hamiltonian describing a 2D mirror Chern insulator
(MCI).

$\displaystyle H(\textbf{k})=k_{x}\Gamma_{1}+k_{y}\Gamma_{2}+M\Gamma_{3},$
(S18)

where $\Gamma_{1}=\tau_{x}\sigma_{x}$, $\Gamma_{2}=\tau_{x}\sigma_{y}$,
$\Gamma_{3}=\tau_{z}$. The mirror operator
$M_{z}=-i\Gamma_{1}\Gamma_{2}\Gamma_{3}=\tau_{z}\sigma_{z}$ satisfies
$M_{z}H(\textbf{k})M_{z}^{-1}=H(\textbf{k})$. Considering the projection
operator $P_{\pm}=\frac{1}{2}(1\pm\tau_{z}\sigma_{z})$, the Hamiltonian can be
expressed as a direct product of two Hamiltonians that have opposite mirror
eigenvalues: $H=P_{+}H\oplus P_{-}H$. Depending on the mass sign, the
Hamiltonian describes either a MCI ($M>0,H=H_{C=+1}\oplus H_{C=-1}$) or a
trivial insulator ($M<0$) Teo et al. (2008). The Hamiltonian with a vortex is
given by

$\displaystyle H(\textbf{k},\theta)$
$\displaystyle=k_{x}\tau_{x}\sigma_{x}+k_{y}\tau_{x}\sigma_{y}$
$\displaystyle+\Delta_{1}\cos\theta\tau_{z}\sigma_{0}+\Delta_{2}\sin\theta\tau_{x}\sigma_{z},$
(S19)

where $\tau_{x}\sigma_{z}$ is odd under mirror symmetry, so that the vortex
structure breaks mirror symmetry except at $\theta=0$ and $\pi$. Similar to
the case with $\mathcal{T}$ symmetry, a vortex structure is realizable at a
junction where a MCI and two mirror symmetry broken insulators meet. As shown
in the Fig. S4 (b), in spin orbit coupled systems, two ferromagnetic
insulators with anti-parallel in-plane magnetization break mirror symmetry,
introducing surface mass terms with opposite sign, so that fractional charges
can be localized at the vortex core. In spinless fermion systems, the
ferromagnetic insulators can be replaced by two ferroelectric insulators with
anti-parallel out-of-plane charge polarization.

$\mathcal{P}$ symmetry.— Except for the cases with $\mathcal{T}$ and
$\mathcal{M}_{z}$ symmetries, the quantized magnetoelectric polarizability
$P_{3}$ of the axion insulators associated with other symmetries can be
described by the pumping process between a higher order topological insulator
and a trivial insulator. The inversion symmetric axion insulator can be
described by the pumping process between a $P$ protected higher order TI and a
trivial insulator. The corresponding low energy Dirac Hamiltonian can be
written as

$\displaystyle
H(k_{x},k_{y},\theta)=k_{x}\Gamma_{1}+k_{y}\Gamma_{2}+M\cos\theta\Gamma_{3}+M\sin\theta\Gamma_{4},$
(S20)

where we take $P=\Gamma_{3}=\tau_{z}$, $\Gamma_{1}=\tau_{x}$,
$\Gamma_{2}=\tau_{y}\sigma_{z}$, and $\Gamma_{4}=\tau_{y}\sigma_{y}$. One can
clearly see that $H(k_{x},k_{y},\theta)$ is $P$ symmetric only at
$\theta=0,~{}\pi$. Since the coefficient of the $\Gamma_{3}$ term has the
opposite sign at $\theta=0$ and $\theta=\pi$, the number of the occupied bands
with the negative $P$ eigenvalue is also different by two. Considering that
the $Z_{2}$ invariant $\nu_{2D}$ of a $P$ symmetric 2D insulator is expressed
by the product of $P$ eigenvalues at time reversal invariant momenta (TRIM)
$\Gamma_{i}$ as

$\displaystyle(-1)^{\nu_{2D}}=\prod_{i=1}^{4}(-1)^{[N_{\text{occ}}^{-}(\Gamma_{i})/2]},$
(S21)

where $N_{\text{occ}}^{-}(\Gamma_{i})$ is the number of the occupied states
with the negative $P$ eigenvalues at the momentum $\Gamma_{i}$. For a 2D
higher order TI, $\nu_{2D}\equiv 1\mod 2$ and the corresponding topological
invariant for the axion insulator is given by
$\nu_{2D}(k_{z}=0)-\nu_{2D}(k_{z}=\pi)\equiv 2P_{3}\mod 2$ Wieder and Bernevig
(2018). The higher order topology of the $\nu_{2D}=1$ phase can be understood
as follows. Using the symmetry representation of $\Gamma$ matrices, the above
Dirac Hamiltonian at $\theta=0,~{}\pi$ can be viewed as two copies of Chern
insulators with the opposite Chern numbers. After suitable regularization by
addition terms quadratic in momentum, one can make $H(\theta=0)$ to have a
pair of counter-propagating edge modes whereas $H(\theta=\pi)$ has no edge
mode. Then by adding $P$ symmetric mass terms to $H(\theta=0)$, one can find
corner charges at the location where the sign of the surface mass terms
changes Khalaf (2018).

$\mathcal{C}_{4z}\mathcal{P}$ symmetry.— The $\mathcal{C}_{4z}\mathcal{P}$
invariant axion insulator can be described by the pumping process between two
different $C_{4z}$ invariant insulators defined on the $k_{z}=0$ and
$k_{z}=\pi$ planes Van Miert and Ortix (2018). For a 2D $C_{4z}$ invariant
insulator, the corner charge $Q_{C}$ is quantized as
$Q_{C}=\frac{N_{C}}{4}~{}\mod 1$, where $N_{C}$ denotes the number of
electrons whose Wannier center is located at the Wyckoff position C. Thus, the
higher order band topology of $2D$ $C_{4z}$ invariant insulators is
characterized by a $\mathbb{Z}_{4}$ invariant. However, a $C_{4z}P$ symmetric
axion insulator is characterized by a $\mathbb{Z}_{2}$ invariant, which
detects the difference of the corner charges defined on the $k_{z}=0$ and
$k_{z}=\pi$ planes, respectively Van Miert and Ortix (2018). The
$\mathbb{Z}_{2}$ nature arises because $C_{2z}$ symmetry exists on every
$k_{z}$ plane so that the set of $C_{2z}$ eigenvalues on the $k_{z}=0$ and
$k_{z}=\pi$ planes must be the same in an insulating phase. Thus, the change
of $N_{C}$ must be an even integer such as

$\displaystyle N_{C}(k_{z}=\pi)-N_{C}(k_{z}=0)=0~{}\text{or}~{}2~{}\mod 4.$
(S22)

Now let us consider the following low energy Dirac Hamiltonian

$\displaystyle
H(\textbf{k},\theta)=k_{x}\Gamma_{1}+k_{y}\Gamma_{2}+M\cos\theta\Gamma_{3}+M\sin\theta\Gamma_{4},$
(S23)

where $\Gamma_{1}=\sigma_{x}$, $\Gamma_{2}=\sigma_{y}$,
$\Gamma_{3}=\tau_{z}\sigma_{z}$, $\Gamma_{4}=\tau_{y}\sigma_{z}$,
$\Gamma_{5}=\tau_{x}\sigma_{z}$, and
$C_{4z}=$diag$[1,i,-1,-i]=\tau_{z}\otimes\begin{pmatrix}1&0\\\
0&i\end{pmatrix}$. Then from $C_{4}\Gamma_{1}C_{4}^{-1}=\Gamma_{2}$ and
$C_{4}\Gamma_{2}C_{4}^{-1}=-\Gamma_{1}$, one can see the $C_{4z}$ invariance
of the Hamiltonian at $\theta=0,~{}\pi$,
$C_{4}H(\textbf{k},\theta=0,~{}\pi)C_{4}^{-1}=H(C_{4}\textbf{k},\theta=0,~{}\pi)$.

At $\theta=0,~{}\pi$, the wave functions of the occupied states are
$[0,1,0,0]^{T}$ and $[0,0,1,0]^{T}$ when $M>0$, with the corresponding
$C_{4z}$ eigenvalues, $i$ and $-1$, respectively. On the other hand, when
$M<0$, the wave functions of the occupied states are $[1,0,0,0]^{T}$ and
$[0,0,0,1]^{T}$ with the corresponding $C_{4z}$ eigenvalues, $1$ and $-i$,
respectively. One can clearly see that the $C_{2z}$ eigenvalues of the
occupied states remains the same independent of the sign of $M$. In this case,
$\Delta N_{C}=2\mod 4$ as shown in Ref. [Van Miert and Ortix, 2018], and the
sign reversal of $M$ describes the nontrivial $\mathbb{Z}_{2}$ invariant of
the axion insulator.

Figure S5: (Color online) (a) Schematic figure describing the lattice
structure around a vortex in a $C_{4z}$ invariant insulator. (b) Schematic
figure describing the $\mathcal{C}_{4z}\mathcal{P}$ symmetric axion insulator
in momentum space. Here the $k_{z}(\theta)=0$ plane corresponds to the $2D$
insulator with a quantized quadrupole moment while the $k_{z}(\theta)=\pi$
plane corresponds to a trivial insulator.

On the other hand, when $\theta\neq 0,~{}\pi$, $C_{4z}$ symmetry is broken
while $C_{2z}$ is preserved. This means that around a vortex core, the system
evolves continuously between two $C_{4z}$ invariant insulators while keeping
the $C_{2z}$ symmetry. Fig. S5 shows an example of the lattice structure
modulation around the vortex. In real materials, since the $U(1)$ symmetry of
the order parameter is reduced to a discrete $C_{n}$ symmetry due to the
lattice potential. An example of $\mathbb{Z}_{4}$ vortex structure is shown in
Fig. S4 (d) in which fractional charges are localized at the junction between
four domains with $\theta=0$, $\pi/4$, $\pi/2$, $3\pi/4$, respectively Chamon
et al. (2008b); Benalcazar et al. (2017b). Here the four domains are obtained
due to four different ways choosing the unit cell for given lattice structure.

$\mathcal{C}_{2z}\mathcal{T}$ symmetry.— Let us first consider spinless
fermion systems. The low energy Dirac Hamiltonian describing a $C_{2z}T$
invariant insulator can be written as

$\displaystyle
H_{D}(\theta)=-i\Gamma_{1}\partial_{1}-i\Gamma_{2}\partial_{2}+M\cos\theta\Gamma_{3}+M\sin\theta\Gamma_{4},$
(S24)

where $\Gamma_{1}=\tau_{x}$, $\Gamma_{2}=\sigma_{y}\tau_{y}$,
$\Gamma_{3}=\tau_{z}$, $\Gamma_{4}=\sigma_{z}\tau_{y}$,
$\Gamma_{5}=\sigma_{x}\tau_{y}$. Here we assume the following symmetry
representations $C_{2z}=\Gamma_{3}=\tau_{z}$, $T=\Gamma_{3}K$, $C_{2z}T=K$
where $K$ denotes complex conjugation. One can easily see that
$H_{D}(\theta=0,~{}\pi)$ is $C_{2z}T$ invariant while
$M\sin\theta\Gamma_{4}\neq 0$ breaks $C_{2z}$ and $C_{2z}T$ while keeping $T$
invariance, and thus it can connect two $C_{2z}T$ invariant insulators with
distinct topological properties.

The higher order nature of $C_{2z}T$ invariant insulators can be seen by
following the similar idea as in the case of $C_{2z}$ symmetric systems.
Namely, $H_{D}(\theta=0,~{}\pi)$ can be considered as two copies of quantum
Hall insulators with opposite Chern numbers. Then after a suitable
regularization, one can make $H_{D}(\theta=0)$ to have a pair of counter-
propagating chiral edge modes whereas $H_{D}(\theta=0,~{}\pi)$ has no edge
state. By adding $C_{2z}T$ symmetric mass terms, corner charges can be found
at the domain wall of the surface mass terms. Recently, it is shown that a
$C_{2z}T$ invariant axion insulator, dubbed a 3D strong Stiefel Whitney
insulator (SWI), can be described by using a pumping process between a 2D SWI
and a trivial insulator with the corresponding $Z_{2}$ invariant $w_{2}=1$ and
$w_{2}=0$, respectively. Also, it is shown that
$w_{2}(k_{z}=0)-w_{2}(k_{z}=\pi)\equiv 2P_{3}\mod 2$ Wieder and Bernevig
(2018); Ahn and Yang (2018).

In spinful fermion systems, one can use a different basis for symmetry
representation. For instance, we can choose $C_{2z}=i\sigma_{y}$,
$T=i\sigma_{y}K$, $C_{2z}T=K$, and also $\Gamma_{1}=\sigma_{x}$,
$\Gamma_{2}=\sigma_{z}$, $\Gamma_{3}=\sigma_{y}\tau_{y}$,
$\Gamma_{4}=\sigma_{y}\tau_{x}$, $\Gamma_{5}=\sigma_{y}\tau_{z}$ where
$\sigma$ denotes the spin degrees of freedom. One can see that $\Gamma_{4}$
and $\Gamma_{5}$ terms break $T$ symmetry while keeping $C_{2z}$. Hence the
vortex structure breaks $T$ when $\theta\neq 0,~{}\pi$ in spinful fermion
systems, whereas it breaks $C_{2z}$ while keeping $T$ in spinless fermion
systems.

$\mathcal{C}_{4z}\mathcal{T}$ symmetry.— The low energy Dirac Hamiltonian
describing a QSHI protected by $C_{4z}$ and $T$ is

$\displaystyle
H=-i\Gamma_{1}\partial_{1}-i\Gamma_{2}\partial_{2}+M\cos\theta\Gamma_{3}+M\sin\theta\Gamma_{4},$
(S25)

where $\Gamma_{1}=\sigma_{x}\tau_{x}$, $\Gamma_{2}=\sigma_{y}\tau_{x}$,
$\Gamma_{3}=\tau_{z}$, $\Gamma_{4}=\sigma_{z}\tau_{x}$, $\Gamma_{5}=\tau_{y}$,
$C_{4}=\cos(\pi/4)+\sin(\pi/4)\Gamma_{1}\Gamma_{2}$ and $T=i\sigma_{y}K$.
Using $C_{4z}\Gamma_{1}C_{4z}^{-1}=-\Gamma_{2}$,
$C_{4z}\Gamma_{2}C_{4z}^{-1}=\Gamma_{1}$, one can easily check the invariance
of the Hamiltonian under $C_{4z}$ and $T$. Also
$\Gamma_{4}=\sigma_{z}\tau_{x}$ and $\Gamma_{5}=\tau_{y}$ break $T$ but
satisfy $C_{4z}$. Hence a constant mass term $M\sin\theta\Gamma_{4}\neq 0$
describing a ferromagnetic ordering with out-of-plane magnetization, breaks
$T$ while keeping $C_{4z}$, and thus can connect two $C_{4z}T$ symmetric
insulators with distinct topological properties.

The higher order nature of the topological insulator protected by $C_{4z}T$
symmetry can be understood as follows. For convenience, let us consider a
disk-shaped finite-size system with a circular boundary. We assume that $M>0$
inside the insulator and $M<0$ outside the insulator. In polar coordinates,
the Hamiltonian is expressed as

$\displaystyle
H=-i\Gamma_{1}(\theta)\partial_{r}-i\Gamma_{2}(\theta)r^{-1}\partial_{\theta}+M(r)\Gamma_{3},$
(S26)

where $\Gamma_{1}(\theta)=\cos\theta\Gamma_{1}+\sin\theta\Gamma_{2}$ and
$\Gamma_{2}=-\sin\theta\Gamma_{1}+\cos\theta\Gamma_{2}$. Employing a
projection operator $P(\theta)=$
$\frac{1}{2}(1+i\Gamma_{1}(\theta)\Gamma_{3})$, we get a gapless surface
Hamiltonian $H|_{r=R}=-iR^{-1}\tilde{\sigma_{z}}\partial_{\theta}$, which is a
characteristic of a QSHI. If $T$ symmetry is broken but non-local $C_{4}T$
symmetry is protected, surface mass terms are allowed:
$H_{m}=m_{4}(\textbf{r})\Gamma_{4}+m_{5}(\textbf{r})\Gamma_{5}$. Projection to
the surface gives
$H_{m}|_{r=R}=m_{4}\tilde{\sigma_{x}}+m_{5}\tilde{\sigma_{y}}$. Under $C_{4}T$
symmetry, mass terms change sign: $m_{i}(\textbf{r})=-m_{i}(C_{4}\textbf{r})$.
Since two masses do not vanish simultaneously in general, the system does not
exhibit anomalous corner states. However, in the presence of chiral symmetry,
for instance, $\Gamma_{5}H\Gamma_{5}^{-1}=-H$, $m_{5}$ must vanish, so that
$C_{4z}$ symmetric corner states appear.

The corresponding lattice Hamiltonian describing a $C_{4}T$ invariant axion
insulator is proposed in Ref. [Schindler et al., 2018] as

$\displaystyle H(\textbf{k})$ $\displaystyle=(M+t\sum_{i}\cos
k_{i})\tau_{z}\sigma_{0}+\Delta_{1}\sum_{i}\sin k_{i}\tau_{x}\sigma_{i}$
$\displaystyle+\Delta_{2}(\cos k_{x}-\cos k_{y})\tau_{y}\sigma_{0}$ (S27)

where $-3<M/t<-1$. Between $k_{z}=0$ and $\pi$, mass term changes sign, which
means that one of them describes a higher order topological insulator and the
other is a trivial insulator. Thus, the axion insulator exhibits a pumping
process between a 2D chiral symmetric higher order topological insulator and a
trivial insulator. Similar to the case with $\mathcal{T}$ symmetry,
$(\mathcal{C}_{4z}\mathcal{T})^{4}=-1$ condition is crucial to define an axion
insulator phase. For instance, let us suppose that there is
$\mathcal{C}_{4z}\mathcal{T}$ symmetric chiral hinge modes for a finite size
system. On the top surface, there must be gapless modes that respect
$\mathcal{C}_{4z}\mathcal{T}$ symmetry. However, gapless modes are not
protected for spinless fermion systems. Whereas gapless modes are protected
for spinful electron systems since $+i$ and $-i$ eigenvalues of
$\mathcal{C}_{2z}$ symmetry are always paired due to
$\mathcal{C}_{4z}\mathcal{T}$ symmetry.

##  S6 Field-theoretical description for a topological vortex

Let us explain the connection between the topological spectral flow and charge
accumulation on a vortex by exploiting field theoretical approach. Since the
number of topological spectral flow lines is equivalent to the winding number
of the vortex, it suffices to show the relation between the winding number and
vortex charge. The low energy Hamiltonian describing a vortex can be written
as follows (we assume that the vortex is protected by inversion symmetry,
without loss of generality) Chamon et al. (2008a); Ryu et al. (2009).

$\displaystyle H=\Psi^{\dagger}[\alpha\cdot{\bf
p}+\beta(\Delta_{1}-\gamma_{5}\Delta_{2})+\mu]\Psi,$ (S28)

where
$p=-i(\partial_{x},\partial_{y}),\alpha=(\sigma_{x}\tau_{z},\sigma_{y}\tau_{z}),\beta=\tau_{x},\gamma_{5}=\tau_{z},\Delta(r)=\Delta_{1}(r)+\Delta_{2}(r)=\Delta_{0}(r)e^{i\theta(r)},\mu$
is the next-nearest interaction term that breaks chiral symmetry
$S=\sigma_{z}\tau_{z}$, and inversion symmetry is represented as $\beta$, such
that the Hamiltonian is invariant under inversion only at $\theta=0,\pi$.
$\theta=0$ corresponds to a higher-order topological insulator (HOTI) phase,
and $\theta=\pi$ corresponds to a trivial insulator phase as in the main text.

Under a unitary transformation, the Hamiltonian can be written as
$H=\begin{pmatrix}\mu&D\\\ D^{\dagger}&\mu\end{pmatrix}$, where
$D=i\sigma^{i}\partial_{i}+i\Delta_{1}+\sigma_{3}\Delta_{2}$, and
$D^{\dagger}=i\sigma^{i}\partial_{i}-i\Delta_{1}+\sigma_{3}\Delta_{2}$ Chamon
et al. (2008a). One can see that there is a bulk gap between $E=-m+\mu$ and
$E=m+\mu$ $(m=|\Delta|)$. If the vortex has a winding number $n$, there are
$n$ bound states between the gap. To prove this, we solve the Schrodinger
equation.

$\displaystyle H\Psi_{E}=E\Psi_{E},\Psi_{E}=(u_{E},v_{E})^{t}$ (S29)

First, let us find the solution where $E=\mu$, such that

$\displaystyle\begin{pmatrix}0&D\\\
D^{\dagger}&0\end{pmatrix}\begin{pmatrix}u_{E}\\\ v_{E}\end{pmatrix}=0.$ (S30)

The equation is equivalent to $Dv_{E}=0,u_{E}=0$ or
$D^{\dagger}u_{E}=0,v_{E}=0.$ According to the index theorem, the number of
solutions $Dv_{E}=0$ and $D^{\dagger}u_{E}=0$ are determined by the winding
number of vortex order parameter $\Delta_{0}(r)e^{i\theta}$: Weinberg (1981)

$\displaystyle n>0:\dim kerD^{\dagger}=0,\quad\dim kerD=n,$ $\displaystyle
n<0:\dim kerD^{\dagger}=|n|,\quad\dim kerD=0.$ (S31)

The index theorem tells us that the number of in-gapped bound states is equal
to the winding number of the vortex.

Second, let us find the solution where $E\neq\mu$. Let $E=\mu+\epsilon$, then
one can find that valence and conduction states are always paired
($E=\mu\pm\epsilon$) since $D^{\dagger}Du_{E}=\epsilon^{2}u_{E}$. Thus, the
number of occupied bound states is equal to the number of unccupied bound
states.

In conclusion, as vortex order parameter winds the vortex core once, a
nontrivial bound state is created. Especially in the presence of chiral
symmetry, the nontrivial state is fixed at the zero energy, which gives an
exact half electric charge localized at the vortex core.

##  S7 Wannier center and quantized Wannier charge in topological vortex

Here we use symmetric Wannier function description to show that fractional
Wannier charge bound to a topological vortex is quantized even without chiral
or particle-hole symmetry. Let us first consider a 1D domain wall system of
the SSH Hamiltonian. The value of polarization should be quantized to $0$ or
$1/2\mod 1$ in the presence of inversion symmetry. When $P_{1}=0$, Wannier
function is located at the center of the unit cell. In the case of
$P_{1}=1/2$, on the other hand, Wannier function is located at the unit cell
boundary. Fig. S6 (a) describes the domain wall structure of the SSH
Hamiltonian between trivial and topological phases. Here we assume that chiral
symmetry is broken due to, for instance, the next nearest neighbor hopping. On
the left-hand side($P_{1}=0$), the Wannier center is located at the unit cell
center. On the other side($P_{1}=1/2$), the Wannier center is located at the
unit cell boundary, where two adjacent cells share an electron: the electron
gives a half charge contribution to each unit cell. Therefore, the domain
structure composed of unit cells contains half-integer electrons due to the
electron located at the unit cell boundary on the right-hand side. Namely,
irrespective of chiral or particle-hole symmetry, a vortex always hosts a
quantized half-integral Wannier charge as long as the relevant symmetry is
preserved. Although the quantized Wannier charge is not identical to the
physical electric charge localized at a domain wall when chiral symmetry is
broken, the half-quantization of the Wannier charge at a domain wall indicates
the fact that the two insulators meeting at the domain carry distinct
quantized topological invariants and therefore cannot be smoothly connected to
each other in the presence of inversion symmetry.

Similarly, the existence of half-quantized charge bound to a topological
vortex can be explained by the Wannier function description. Except for $T$
and $M_{z}$, there are four symmetries $P$ (or $C_{2z}$), $C_{2z}T,C_{4z},$
and $C_{4z}T$ that are relevant to the vortices with higher order topology.
Since $T$ is local-in-space symmetry, it does not change the position of
Wannier center. Thus, to prove that half-integer charge is bound to a
topological vortex, we only have to consider $P$ (or $C_{2z}$) and $C_{4z}$
symmetries.

Figure S6: (Color online) (a) Domain wall structure between trivial
($P_{1}=0$) and topological($P_{1}=1/2$) insulator. Since Wannier function is
located at the unit cell boundary for $P_{1}=1/2$, the domain wall structure
composed of integer unit cells (shaded region) contains half integer
electrons. (b) Wyckoff positions of $P$ and $C_{4z}$ symmetric lattices. (c)
Topological vortex protected by inversion symmetry. For $\theta=0$ (inversion
symmetric HOTI), $N_{B}=N_{C}=N_{D}=1\mod 2$. Since Wannier function is
located at the unit cell boundary at $\theta=0$, the charge bound to the
vortex is quantized to half-integer. (d) Topological vortex protected by
$C_{4z}$ symmetry. Since $N_{C}(\theta=0)=2\mod 4$ ($C_{4z}$ symmetric HOTI),
the charge bound to the vortex is
$\frac{1}{4}(2N_{C}(\theta=0)+2N_{c}(\theta=\pi))=\frac{1}{2}\mod 1$.

$P$ (or $C_{2z}$).— Inversion symmetric insulator has four Wyckoff positions
A, B, C, and D which are invariant under inversion up to lattice translation
vector (Fig. S6 (a)). We define the number of symmetric Wannier functions
centered at each Wyckoff position as $N_{A},N_{B},N_{C},$ and $N_{D}$. For
inversion symmetric HOTI, polarization is zero:

$\displaystyle p_{x}=\sum_{i}X_{i}=\frac{1}{2}(N_{B}+N_{C})=0\mod 1,$ (S32)
$\displaystyle p_{y}=\sum_{i}Y_{i}=\frac{1}{2}(N_{C}+N_{D})=0\mod 1,$ (S33)

where $X_{i}$ and $Y_{i}$ denote $x$ and $y$ coordinates of $i$th Wannier
function, respectively. On the other hand, quadrupole moment is defined as

$\displaystyle
q_{xy}=\sum_{i}X_{i}Y_{i}=\frac{1}{4}N_{C}=\frac{1}{4}\mod\frac{1}{2}.$ (S34)

Let us note that $N_{C}$ is well defined modulo 2 since $N_{C}$ can be changed
only by multiples of two keeping inversion symmetry. Taken altogether,
$N_{B}=N_{C}=N_{D}=1\mod 2$. For a trivial insulator, on the other hand,
quadrupole moment $q_{xy}=\frac{1}{4}N_{C}=0\mod\frac{1}{2}$. Thus,
$N_{B}=N_{C}=N_{D}=0\mod 2$. Let us consider a vortex structure
$PH(k_{x},k_{y},\theta)P^{-1}=H(-k_{x},-k_{y},-\theta)$ where $\theta=0$
$(\pi)$ corresponds to inversion symmetric HOTI (trivial insulator) (Fig. S6
(c)). Since an odd number of Wannier functions are located at the unit-cell
boundary for $\theta=0$, the vortex geometry composed of integer unit-cells
contains a half-integer charge.

$C_{4z}$.— $C_{4z}$ symmetric insulator has two Wyckoff positions A and C
which are invariant under $C_{4z}$ operation and two Wyckoff positions B and D
which are invariant under $C_{2z}$ operation up to lattice translation vector
(Fig. S6 (b)). The vortex structure satisfies
$C_{4z}H(k_{x},k_{y},\theta)C_{4z}^{-1}=H(k_{x}^{\prime},k_{y}^{\prime},-\theta)$,
followed by $C_{2z}H(k_{x},k_{y},\theta)C_{2z}^{-1}=H(-k_{x},-k_{y},\theta)$.
Since we are not interested in inversion symmetry protected topological phase,
$N_{B}=N_{C}=N_{D}=0\mod 2$ for any $\theta$. As shown earlier, in the case of
$C_{4z}$ symmetric HOTI, $N_{C}=2\mod 4$. Also, $N_{C}=0\mod 4$ for trivial
insulator. Note that $N_{C}$ is well defined modulo 4 since $N_{C}$ can be
changed by multiples of four keeping $C_{4}$ symmetry. The charge bound to the
vortex is, therefore,
$\frac{1}{4}[N_{C}(\theta=0)+N_{c}(\theta=\pi)]=\frac{1}{2}\mod 1$, for the
electron located at the Wyckoff position $C$ gives $1/4$ charge contribution
to each unit cell (Fig. S6 (d)).

## References

  * Su et al. (1980) W.-P. Su, J. R. Schrieffer, and A. J. Heeger, Physical Review B 22, 2099 (1980).
  * Jackiw and Rebbi (1976) R. Jackiw and C. Rebbi, Physical Review D 13, 3398 (1976).
  * Jackiw and Rossi (1981) R. Jackiw and P. Rossi, Nuclear Physics B 190, 681 (1981).
  * Goldstone and Wilczek (1981) J. Goldstone and F. Wilczek, Physical Review Letters 47, 986 (1981).
  * Read and Green (2000) N. Read and D. Green, Physical Review B 61, 10267 (2000).
  * Teo and Kane (2010a) J. C. Y. Teo and C. L. Kane, Physical review letters 104, 046401 (2010a).
  * Teo and Kane (2010b) J. C. Y. Teo and C. L. Kane, Physical Review B 82, 115120 (2010b).
  * Hou et al. (2007) C.-Y. Hou, C. Chamon, and C. Mudry, Physical review letters 98, 186809 (2007).
  * Chamon et al. (2008a) C. Chamon, C.-Y. Hou, R. Jackiw, C. Mudry, S.-Y. Pi, and G. Semenoff, Physical Review B 77, 235431 (2008a).
  * Chamon et al. (2008b) C. Chamon, C.-Y. Hou, R. Jackiw, C. Mudry, S.-Y. Pi, and A. P. Schnyder, Physical review letters 100, 110405 (2008b).
  * Ryu et al. (2009) S. Ryu, C. Mudry, C.-Y. Hou, and C. Chamon, Physical Review B 80, 205319 (2009).
  * Hughes et al. (2011) T. L. Hughes, E. Prodan, and B. A. Bernevig, Physical Review B 83, 245132 (2011).
  * Morimoto and Furusaki (2014) T. Morimoto and A. Furusaki, Physical Review B 89, 235127 (2014).
  * Qi et al. (2008) X.-L. Qi, T. L. Hughes, and S.-C. Zhang, Physical Review B 78, 195424 (2008).
  * Noh et al. (2018) J. Noh, W. A. Benalcazar, S. Huang, M. J. Collins, K. P. Chen, T. L. Hughes, and M. C. Rechtsman, Nature Photonics 12, 408 (2018).
  * Benalcazar et al. (2018) W. A. Benalcazar, T. Li, and T. L. Hughes, Physical Review B 99, 245151 (2019).
  * Zak (1989) J. Zak, Physical review letters 62, 2747 (1989).
  * Xiao et al. (2010) D. Xiao, M.-C. Chang, and Q. Niu, Reviews of modern physics 82, 1959 (2010).
  * Ahn et al. (2018a) J. Ahn, D. Kim, Y. Kim, and B.-J. Yang, Physical review letters 121, 106403 (2018a).
  * Ahn and Yang (2018) J. Ahn and B.-J. Yang, Physical Review B 99, 235125 (2019).
  * Ahn et al. (2018b) J. Ahn, S. Park, and B.-J. Yang, Physical Review X 9, 021013 (2019).
  * Wieder and Bernevig (2018) B. J. Wieder and B. A. Bernevig, arXiv preprint arXiv:1810.02373 (2018).
  * Song et al. (2018) Z. Song, T. Zhang, and C. Fang, Physical Review X 8, 031069 (2018).
  * Fu and Kane (2007) L. Fu and C. L. Kane, Physical Review B 76, 045302 (2007).
  * Kim et al. (2015) Y. Kim, B. J. Wieder, C. L. Kane, and A. M. Rappe, Physical review letters 115, 036806 (2015).
  * Turner et al. (2012) A. M. Turner, Y. Zhang, R. S. K. Mong, and A. Vishwanath, Physical Review B 85, 165120 (2012).
  * Po et al. (2017) H. C. Po, A. Vishwanath, and H. Watanabe, Nature Communications 8, 50 (2017).
  * Geier et al. (2018) M. Geier, L. Trifunovic, M. Hoskam, and P. W. Brouwer, Physical Review B 97, 205135 (2018).
  * Khalaf (2018) E. Khalaf, Physical Review B 97, 205136 (2018).
  * SM (2020) See Supplemental Material at [URL] for more details on nested Wilson loop, higher order band topology, classification of axion insulators, field-theoretical description for a topological vortex, and Wannier charge.
  * Benalcazar et al. (2017b) W. A. Benalcazar, B. A. Bernevig, and T. L. Hughes, Physical Review B 96, 245115 (2017b).
  * Xie et al. (2018) B.-Y. Xie, H.-F. Wang, H.-X. Wang, X.-Y. Zhu, J.-H. Jiang, M.-H. Lu, and Y.-F. Chen, Physical Review B 98, 205147 (2018).
  * Weinberg (1981) E. J. Weinberg, Physical Review D 24, 2669 (1981).
  * Varjas et al. (2015) D. Varjas, F. de Juan, and Y.-M. Lu, Physical Review B 92, 195116 (2015).
  * Schindler et al. (2018) F. Schindler, Z. Wang, M. G. Vergniory, A. M. Cook, A. Murani, S. Sengupta, A. Y. Kasumov, R. Deblock, S. Jeon, I. Drozdov, et al., Nature physics 14, 918 (2018).
  * Ryu et al. (2010) S. Ryu, A. P. Schnyder, A. Furusaki, and A. W. Ludwig, New Journal of Physics 12, 065010 (2010).
  * Fu and Kane (2006) L. Fu and C. L. Kane, Physical Review B 74, 195312 (2006).
  * Bernevig et al. (2006) B. A. Bernevig, T. L. Hughes, and S.-C. Zhang, Science 314, 1757 (2006).
  * Teo et al. (2008) J. C. Y. Teo, L. Fu, and C. L. Kane, Physical Review B 78, 045426 (2008).
  * Van Miert and Ortix (2018) G. van Miert and C. Ortix, Physical Review B 98, 081110(R) (2018).
  * Schindler et al. (2018) F. Schindler, A. M. Cook, M. G. Vergniory, Z. Wang, S. S. Parkin, B. A. Bernevig, and T. Neupert, Science advances 4, eaat0346 (2018).
  * Călugăru et al. (2019) D. Călugăru, V. Juričić, and B. Roy, Physical Review B 99, 041301(R) (2019).
  * Okuma et al. (2019) N. Okuma, M. Sato, and K. Shiozaki, Physical Review B 99, 085127 (2019).
  * Song et al. (2017) Z. Song, Z. Fang, and C. Fang, Physical review letters 119, 246402 (2017).
  * Lee et al. (2020) E. Lee, R. Kim, J. Ahn, and B.-J. Yang, npj Quantum Materials 5, 1 (2020).
  * Hwang et al. (2019) Y. Hwang, J. Ahn, and B.-J. Yang, Physical Review B 100, 205126 (2019).
  * Qi et al. (2008b) X.-L. Qi, T. L. Hughes, and S.-C. Zhang, Nature Physics 4, 273 (2008b).