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

TL;DR

This paper links the fractional charge bound to vortices in 2D topological crystalline insulators, like Kekule textured graphene, to bulk topological invariants and explores their relation to 3D axion insulators and symmetry protections.

Contribution

It establishes a theoretical framework connecting vortex-bound fractional charge in 2D topological crystalline insulators to bulk topology and 3D axion insulator models, including symmetry considerations.

Findings

Fractional charge at vortices relates to bulk topological invariants.

Effective 3D Hamiltonian describes axion insulator with quantized magnetoelectric polarization.

Vortex charge is half-quantized under chiral symmetry, otherwise it is not quantized but remains half-quantized with inversion symmetry.

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.

Figures

Tables

Table 1. Table 1: The correspondence between the symmetry of H 3D ​ ( 𝒌 ) subscript 𝐻 3D 𝒌 H_{\text{3D}}(\bm{k}) exhibiting quantized P 3 subscript 𝑃 3 P_{3} and that of the 2D Hamiltonian H ( k x , k y ; θ = 0 , π ) H(k_{x},k_{y};\theta=0,\pi) . For T 𝑇 T and C 4 ​ T subscript 𝐶 4 𝑇 C_{4}T symmetry, spin orbit coupling (SOC) is necessary for having P 3 = 1 / 2 subscript 𝑃 3 1 2 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
$𝒫={𝒞}_{2z}{ℳ}_z$	$P$, $C_{2z}$, $C_{6z}$		HOTI
${𝒞}_{4z}𝒫={\left({𝒞}_{4z}{ℳ}_z\right)}^3$	$C_{4z}$		HOTI
$𝒯$	$T$	✓	QSHI
${ℳ}_z$	$M_z$		Mirror TI
${𝒞}_{2z}𝒯$	$PT$, $C_{2z}T$, $C_{6z}T$		HOTI
${𝒞}_{4z}𝒯$	$C_{4z}T$	✓	HOTI

Table 2. 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 𝐷 D -dimensional systems can be explained by the topological response of the D + 1 𝐷 1 D+1 dimensional effective Hamiltonian H D + 1 subscript 𝐻 𝐷 1 H_{D+1} .

	Domain wall of SSH model	Vortex of HCM model
$H_{D+1}$	Chern insulator	Axion insulator
Invariant	$C_1=\Delta \left(2P_1\right)$	$2P_3=\Delta {\nu }_{2D}$

Equations

$$H=\sum _{\text{r}\in {\Lambda }_A}\sum _{j=1}^3\left(t+\delta t_j(\text{r})\right)c_A^{†}(\text{r})c_B(\text{r}+{\text{d}}_j)+\text{h.c.},$$

$$(-1{)}^{{\nu }_{2D}}=\prod _{i=1}^4(-1{)}^{[N_{\text{occ}}^{-}({\Gamma }_i)/2]},$$

$$PH(k_x,k_y;\theta )P^{-1}=H(-k_x,-k_y;-\theta ).$$

$$\Delta {\nu }_{2D}=2P_3=\frac{1}{4{\pi }^2}{\int }_{T^2\times S^1}\text{Tr}\left[AdA-\frac{2i}{3}{A}^3\right]\mathrm{mod}2,$$

$$C_{n}H(k_x,k_y;\theta )C_n^{-1}=H(k_x^{\prime },k_y^{\prime };-\theta ),$$

$$H_{\text{SSH}}=\sum (t+\delta t_i)a_i^{†}{a}_{i+1}+H.c,$$

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

$$-i{\partial }_{x}u(x)+\Delta (x)v(x)=0,$$

$${\Delta }^{\ast }(x)u(x)+i{\partial }_{x}v(x)=0.$$

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

$$q_{\text{bound}}=\int -{\partial }_{\theta }{P}_{1}d\theta =1/2\text{ mod 1},$$

$$PH(k_x,\theta )P^{-1}=H(-k_x,-\theta ).$$

$$P_{2D}H(k_x,k_y)P_{2D}^{-1}=H(-k_x,-k_y).$$

$$(-1{)}^C=(-1{)}^{\sum N_{\text{occ}}^{-}({\Gamma }_i)}=(-1{)}^{2P_1(\theta =\pi )-2P_1(\theta =0)},$$

$$W_{(k_1+2\pi ,k_2)\leftarrow (k_1,k_2)}$$

$$=\underset{N\to \infty }{\lim }{F}_{N-1}{F}_{N-2}\cdots F_1{F}_0$$

$$=P{e}^{-i{\oint }_C{A}_{k}dk},$$

$$p_x$$

$$=-\frac{i}{2\pi N_y}\sum _{k_y}\log \det [W_{x,\text{k}}]=-\frac{1}{2\pi }\oint \text{Tr}[A_{x,k}]d^{2}k.$$

$${\tilde{W}}_{k_2+2\pi \leftarrow k_2}={\tilde{W}}_2=\underset{N\to \infty }{\lim }{\tilde{F}}_{N-1}{\tilde{F}}_{N-2}\cdots {\tilde{F}}_1{\tilde{F}}_0,$$

$$H=\sum _{n=1}^9{t}_{\sigma }{a}_{2i-1}^{†}{a}_{2i}+H.c,$$

$$H_{edge}=-i{R}^{-1}{\tilde{\sigma }}_z{\partial }_{\theta }+{\tilde{m}}_1(\theta ){\tilde{\sigma }}_x,$$

$$w=\int d\theta \frac{1}{2\pi }{\partial }_{\theta }{\tan }^{-1}\left(\frac{m_2(\theta )}{m_1(\theta )}\right)$$

$$H(\text{k})=k_x{\tau }_x{\sigma }_x+k_y{\tau }_x{\sigma }_y+M{\tau }_z{\sigma }_0,$$

$$H(k_x,k_y,\theta )$$

$$+\Delta \cos \theta {\tau }_z{\sigma }_0+\Delta \sin \theta {\tau }_x{\sigma }_z,$$

$$H_0$$

$$=v_F\int dx{\Psi }^{†}i{\tilde{\sigma }}_z{\partial }_x\Psi ,$$

$$H(\text{k})=k_x{\Gamma }_1+k_y{\Gamma }_2+M{\Gamma }_3,$$

$$H(\text{k},\theta )$$

$$+{\Delta }_1\cos \theta {\tau }_z{\sigma }_0+{\Delta }_2\sin \theta {\tau }_x{\sigma }_z,$$

$$H(k_x,k_y,\theta )=k_x{\Gamma }_1+k_y{\Gamma }_2+M\cos \theta {\Gamma }_3+M\sin \theta {\Gamma }_4,$$

$$(-1{)}^{{\nu }_{2D}}=\prod _{i=1}^4(-1{)}^{[N_{\text{occ}}^{-}({\Gamma }_i)/2]},$$

$$N_C(k_z=\pi )-N_C(k_z=0)=0\text{or}2\mathrm{mod}4.$$

Full text

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.

II Higher-order band topology of KTG

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

[TABLE]

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

[TABLE]

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 [math] 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 [math] 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

[TABLE]

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,

[TABLE]

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 [math] 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).

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.

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

[TABLE]

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.

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

[TABLE]

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

[TABLE]

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,

[TABLE]

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 [math] 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

[TABLE]

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

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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 [math] 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.

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

[TABLE]

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

[TABLE]

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 [math] 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

[TABLE]

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).

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 [math] 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

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,

[TABLE]

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

[TABLE]

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.

[TABLE]

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.

$\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

[TABLE]

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

[TABLE]

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

[TABLE]

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).

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

Now let us consider the following low energy Dirac Hamiltonian

[TABLE]

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.

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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).

[TABLE]

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.

[TABLE]

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

[TABLE]

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)

[TABLE]

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 [math] 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.

$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:

[TABLE]

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

[TABLE]

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)).