Fragile Topology Protected by Inversion Symmetry: Diagnosis, Bulk-Boundary Correspondence, and Wilson Loop

TL;DR

This paper investigates fragile topological insulators protected by inversion symmetry, proposing a new diagnosis method, analyzing bulk-boundary correspondence, and exploring Wilson loop spectra, revealing their prevalence and complex boundary phenomena.

Contribution

It introduces an efficient symmetry-based diagnosis for fragile topology, extends the concept of higher-order TIs, and links adiabatic pumping to fragile topological phases.

Findings

At least 17% of 2D inversion-symmetric insulators have fragile topology.

Fragile TIs can be viewed as higher-order TIs with boundary mass terms carrying odd winding numbers.

Wilson loop spectra determine the minimal bands needed to resolve Wannier obstructions.

Abstract

We study the bulk and boundary properties of fragile topological insulators (TIs) protected by inversion symmetry, mostly focusing on the class A of the Altland-Zirnbauer classification. First, we propose an efficient method for diagnosing fragile band topology by using the symmetry data in momentum space. Using this method, we show that among all the possible parity configurations of inversion-symmetric insulators, at least 17 percent of them have fragile topology in 2D while fragile TIs are less than 3 percent in 3D. Second, we study the bulk-boundary correspondence of fragile TIs protected by inversion symmetry. In particular, we generalize the notion of $d$-dimensional ($d$D) $k$th-order TIs, which is normally defined for $0<k\le d$, to the cases with $k>d$, and show that they all have fragile topology. In terms of the Dirac Hamiltonian, a $d$D $k$th-order TI has $(k-1)$ boundary mass terms. We show that a minimal fragile TI with the filling anomaly can be considered as the $d$D $(d+1)$th-order TI, and all the other $d$D $k$th-order TIs with $k>(d+1)$ can be constructed by stacking $d$D $(d+1)$th-order TIs. Although $d$D $(d+1)$th-order TIs have no in-gap states, the boundary mass terms carry an odd winding number along the boundary, which induces localized charges on the boundary at the positions where the boundary mass terms change abruptly. In the cases with $k>(d+1)$, we show that the net parity of the system with boundaries can distinguish topological insulators and trivial insulators. Also, by studying the (nested) Wilson loop spectra, we determine the minimal number of bands to resolve the Wannier obstruction, which is consistent with the prediction from our diagnosis method of fragile topology. Finally, we show that a $(d+1)$D $(k-1)$th-order TI can be obtained by an adiabatic pumping of $d$D $k$th-order TI, which generalizes the previous study of the 2D 3rd-order TI.

Tables

Table 1. Table 1: The parity of the occupied state at inversion-invariant momenta of inversion-symmetric atomic insulators. s 𝑠 s and p 𝑝 p denote the even- and odd-parity states, respectively. A 𝐴 A and B 𝐵 B are Wyckoff positions of 1 1 1 D inversion-symmetric lattices shown in Fig.  1 (a).

Orbital	$k=0$	$k=\pi$
$s\left(A\right)$	$+$	$+$
$p\left(A\right)$	$-$	$-$
$s\left(B\right)$	$+$	$-$
$p\left(B\right)$	$-$	$+$

Table 2. Table 2: The relation between N + subscript 𝑁 N_{+} , N − subscript 𝑁 N_{-} , N cell subscript 𝑁 cell N_{\rm cell} , and ℐ ℐ \mathcal{I} for doubled SSH models of finite-size, which corresponds to the 1 1 1 D third and fourth-order topological insulator described by ℋ ( 1 , 3 ) subscript ℋ 1 3 \mathcal{H}_{(1,3)} and ℋ ( 1 , 4 ) subscript ℋ 1 4 \mathcal{H}_{(1,4)} .

$N_{\mathrm{cell}}$	Phase	$N_{+}$	$N_{-}$	$ℐ$
Even	Trivial	$N_{\mathrm{cell}}$	$N_{\mathrm{cell}}$	$0$
	Topological	$N_{\mathrm{cell}}-1$	$N_{\mathrm{cell}}+1$	$2$
Odd	Trivial	$N_{\mathrm{cell}}+1$	$N_{\mathrm{cell}}-1$	$-2$
	Topological	$N_{\mathrm{cell}}$	$N_{\mathrm{cell}}$	$0$

Table 3. Table 3: The relation between N + subscript 𝑁 N_{+} , N − subscript 𝑁 N_{-} , N ± subscript 𝑁 plus-or-minus N_{\pm} , and ℐ ℐ \mathcal{I} for the 1 1 1 D second-order TI described by ℋ ( 1 , 2 ) subscript ℋ 1 2 \mathcal{H}_{(1,2)} when the extra charge δ 𝛿 \delta is added due to the filling anomaly.

$N_{\mathrm{cell}}$	Phase	$N_{+}$	$N_{-}$	$ℐ$
Even	Trivial	$\frac{1}{2}{N}_{\mathrm{cell}}$	$\frac{1}{2}{N}_{\mathrm{cell}}$	$0$
	Topological	$\frac{1}{2}\left(N_{\mathrm{cell}}+\delta -1\right)$	$\frac{1}{2}\left(N_{\mathrm{cell}}+\delta +1\right)$	$1$
Odd	Trivial	$\frac{1}{2}\left(N_{\mathrm{cell}}+1\right)$	$\frac{1}{2}\left(N_{\mathrm{cell}}-1\right)$	$-1$
	Topological	$\frac{1}{2}\left(N_{\mathrm{cell}}+\delta \right)$	$\frac{1}{2}\left(N_{\mathrm{cell}}+\delta \right)$	$0$

Table 4. Table 4: Results of the parity configuration diagnosis in 2D. When n occ subscript 𝑛 occ n_{\rm occ} is sufficiently large, N ​ ( A ) 𝑁 𝐴 N(A) , N ​ ( F ) 𝑁 𝐹 N(F) , and N ​ ( S ​ t ) 𝑁 𝑆 𝑡 N(St) approach to 1 3 ​ N tot 1 3 subscript 𝑁 tot \frac{1}{3}N_{\rm tot} , 1 6 ​ N tot 1 6 subscript 𝑁 tot \frac{1}{6}N_{\rm tot} , and 1 2 ​ N tot 1 2 subscript 𝑁 tot \frac{1}{2}N_{\rm tot} . At least, two-third of parity configurations are topologically nontrivial.

$n_{\mathrm{occ}}$	$N_{\mathrm{tot}}$	$N\left(A\right)$	$N\left(F\right)$	$N\left(St\right)$
$1$	$16$	$8$	$0$	$8$
$2$	$81$	$33$	$8$	$40$
$3$	$256$	$96$	$32$	$128$
$4$	$625$	$225$	$88$	$312$
$5$	$1296$	$456$	$192$	$648$
$6$	$2401$	$833$	$368$	$1200$
$7$	$4096$	$1408$	$640$	$2048$
$8$	$6561$	$2241$	$1040$	$3280$
$9$	$10000$	$3400$	$1600$	$5000$
$10$	$14641$	$4961$	$2360$	$7320$
$n_{\mathrm{occ}}\to \infty$	${\left(n_{\mathrm{occ}}+1\right)}^4$	$\frac{1}{3}{N}_{\mathrm{tot}}$	$\frac{1}{6}{N}_{\mathrm{tot}}$	$\frac{1}{2}{N}_{\mathrm{tot}}$

Table 5. Table 5: Results of the parity configuration diagnosis in 3D. When n occ subscript 𝑛 occ n_{\rm occ} is sufficiently large, N ​ ( A ) 𝑁 𝐴 N(A) is 1 315 ​ N tot 1 315 subscript 𝑁 tot \frac{1}{315}N_{\rm tot} . As for N ​ ( F ) 𝑁 𝐹 N(F) and N ​ ( S ​ t ) 𝑁 𝑆 𝑡 N(St) , although their exact numbers in the n occ → ∞ → subscript 𝑛 occ n_{\rm occ}\rightarrow\infty limit are not known, considering their n occ subscript 𝑛 occ n_{\rm occ} dependence, we conjecture that N ​ ( F ) = 283 10080 ​ N tot 𝑁 𝐹 283 10080 subscript 𝑁 tot N(F)=\frac{283}{10080}N_{\rm tot} and N ​ ( S ​ t ) = 31 32 ​ N tot 𝑁 𝑆 𝑡 31 32 subscript 𝑁 tot N(St)=\frac{31}{32}N_{\rm tot} in the large n occ subscript 𝑛 occ n_{\rm occ} limit. In contrast to 2 2 2 D, the majority of parity configurations are stable topological phases. Here, z 4 = − 4 ​ ν W 1 subscript 𝑧 4 4 subscript 𝜈 subscript 𝑊 1 z_{4}=-4\nu_{W_{1}} denotes the ℤ 4 subscript ℤ 4 \mathbb{Z}_{4} symmetry indicator Po  et al. ( 2017 ); Khalaf  et al. ( 2018 ); Song  et al. ( 2018a ) for inversion symmetry. The relation between symmetry indicators and ν W subscript 𝜈 𝑊 \nu_{W} indices is discussed in Appendix  C . Let us also note that stable topological phases are Weyl semimetal when z 4 subscript 𝑧 4 z_{4} is odd  Turner  et al. ( 2012 ); Hughes  et al. ( 2011 ) .

$n_{\mathrm{occ}}$	$N_{\mathrm{tot}}$	$N\left(A\right)$	$N\left(F\right)$	$N\left(St\right)$
				$z_4=0$	$z_4=1$	$z_4=2$	$z_4=3$
$1$	$256$	$16$	$0$	$56$	$64$	$56$	$64$
$2$	$6561$	$129$	$112$	$1400$	$1640$	$1640$	$1640$
$3$	$65536$	$704$	$1344$	$14336$	$16384$	$16384$	$16384$
$4$	$390625$	$2945$	$9536$	$85176$	$97656$	$97656$	$97656$
$5$	$1679616$	$10128$	$42368$	$367416$	$419904$	$419896$	$419904$
$6$	$5764801$	$29953$	$151248$	$1260000$	$1441200$	$1441200$	$1441200$
$7$	$16777216$	$78592$	$445696$	$3670016$	$4194304$	$4194304$	$4194304$
$8$	$43046721$	$187137$	$1160944$	$9413600$	$10761680$	$10761680$	$10761680$
$9$	$100000000$	$411280$	$2713728$	$21875000$	$25000000$	$24999992$	$25000000$
$10$	$214358881$	$845185$	$5859936$	$46884600$	$53589720$	$53589720$	$53589720$

Table 6. Table 6: The extra charge accumulation in the 2 2 2 D third-order TI. Q i = 1 , 2 , 3 , 4 num . subscript superscript 𝑄 num 𝑖 1 2 3 4 Q^{\rm num.}_{i=1,2,3,4} ( Q i = 1 , 2 , 3 , 4 subscript 𝑄 𝑖 1 2 3 4 Q_{i=1,2,3,4} ) indicates the extra charge in the i 𝑖 i th quadrant computed from the lattice model (from the induced charge of the continuum model). Note that Q 1 num . = Q 3 num . subscript superscript 𝑄 num 1 subscript superscript 𝑄 num 3 Q^{\rm num.}_{1}=Q^{\rm num.}_{3} and Q 2 num . = Q 4 num . subscript superscript 𝑄 num 2 subscript superscript 𝑄 num 4 Q^{\rm num.}_{2}=Q^{\rm num.}_{4} due to inversion symmetry, and similar relations hold for Q i subscript 𝑄 𝑖 Q_{i} .

$Q_1^{\mathrm{num}.}$	$Q_1$	$Q_2^{\mathrm{num}.}$	$Q_2$
$-0.282329$	$-0.284825$	$-0.217611$	$-0.215175$

Table 7. Table 7: The extra charge accumulation in the 3 3 3 D fourth-order TI. Q i = 1 , … , 8 num . subscript superscript 𝑄 num 𝑖 1 … 8 Q^{\rm num.}_{i=1,\dots,8} ( Q i = 1 , … , 8 subscript 𝑄 𝑖 1 … 8 Q_{i=1,\dots,8} ) indicates the extra charge in the i 𝑖 i th octant computed from the lattice model (from the induced charge of the continuum model). Note that Q 1 num . = Q 7 num . subscript superscript 𝑄 num 1 subscript superscript 𝑄 num 7 Q^{\rm num.}_{1}=Q^{\rm num.}_{7} , Q 2 num . = Q 8 num . subscript superscript 𝑄 num 2 subscript superscript 𝑄 num 8 Q^{\rm num.}_{2}=Q^{\rm num.}_{8} , Q 3 num . = Q 5 num . subscript superscript 𝑄 num 3 subscript superscript 𝑄 num 5 Q^{\rm num.}_{3}=Q^{\rm num.}_{5} , and Q 4 num . = Q 6 num . subscript superscript 𝑄 num 4 subscript superscript 𝑄 num 6 Q^{\rm num.}_{4}=Q^{\rm num.}_{6} due to inversion symmetry, and similar relations hold for Q i subscript 𝑄 𝑖 Q_{i} .

$Q_1^{\mathrm{num}.}$	$Q_1$	$Q_2^{\mathrm{num}.}$	$Q_2$
$0.116481$	$0.113613$	$0.203660$	$0.223673$
$Q_3^{\mathrm{num}.}$	$Q_3$	$Q_4^{\mathrm{num}.}$	$Q_4$
$0.089725$	$0.081586$	$0.090134$	$0.081127$

Table 8. Table 8: Change of parity configuration and the Wilson loop spectrum after adding trivial orbitals.

Band rep.		$\left\{n_{\theta }^W\left(k_2=0\right)\right\}$				$\left\{n_{\theta }^W\left(k_2=\pi \right)\right\}$				$N_W$
		$n_0^W\left(0\right)$	$n_{\pi }^W\left(0\right)$	$n_c^W\left(0\right)$		$n_0^W\left(\pi \right)$	$n_{\pi }^W\left(\pi \right)$	$n_c^W\left(\pi \right)$
${ℬ}_{\gamma }(d,2^{d-1})\oplus p\left(A\right)$		$1$	$2^{d-1}$	$0$		$2^{d-1}-1$	$0$	$2$		$2^{d-2}-1$
${ℬ}_{\gamma }(d,2^{d-1})\oplus p\left(B\right)$		$0$	$2^{d-1}+1$	$0$		$2^{d-1}$	$1$	$0$		$2^{d-2}$
${ℬ}_{\gamma }(d,2^{d-1})\oplus p\left(C\right)$		$0$	$2^{d-1}+1$	$0$		$2^{d-1}$	$1$	$0$		$2^{d-2}$
${ℬ}_{\gamma }(d,2^{d-1})\oplus p\left(D\right)$		$1$	$2^{d-1}$	$0$		$2^{d-1}+1$	$0$	$0$		$2^{d-2}$

Table 9. Table 9: Classification of generic higher-order TIs/TSCs protected by inversion symmetry  Khalaf ( 2018 ) . Among the ten AZ classes, class A and AIII belong to complex AZ classes while other eight classes belong to real AZ classes. s 𝑠 s labels the AZ classes: s = 0 , 1 𝑠 0 1 s=0,1 for complex AZ classes and s = 0 , 1 , … , 7 𝑠 0 1 … 7 s=0,1,\dots,7 for real AZ classes. T 𝑇 T , P 𝑃 P , and S 𝑆 S denote time-reversal, particle-hole and chiral symmetries. T 2 superscript 𝑇 2 T^{2} are denoted as 0 0 when the symmetry is absent. The same notation is adopted for P 𝑃 P and S 𝑆 S . Inversion symmetry operator I 𝐼 I satisfies [ I , T ] = { I , P } = 0 𝐼 𝑇 𝐼 𝑃 0 [I,T]=\{I,P\}=0 . δ 𝛿 \delta is given by dimension d 𝑑 d and order k 𝑘 k , δ = d − k + 1 𝛿 𝑑 𝑘 1 \delta=d-k+1 (mod 8 8 8 ). Here, k 𝑘 k is an arbitrary positive integer. As for ℤ 2 subscript ℤ 2 \mathbb{Z}_{2} and ` ​ ` ​ ℤ 2 ​ " ` ` subscript ℤ 2 " ``\mathbb{Z}_{2}" classifications, ℤ 2 subscript ℤ 2 \mathbb{Z}_{2} and ` ​ ` ​ ℤ 2 ​ " ` ` subscript ℤ 2 " ``\mathbb{Z}_{2}" nontrivial classes exhibit gapless boundary states when k ≤ d 𝑘 𝑑 k\leq d . When k = d + 1 𝑘 𝑑 1 k=d+1 , only ` ​ ` ​ ℤ 2 ​ " ` ` subscript ℤ 2 " ``\mathbb{Z}_{2}" nontrivial classes exhibit the filling anomaly.

AZ class	$s$	$T^2$	$P^2$	$S^2$	$\delta =0$	$\delta =1$	$\delta =2$	$\delta =3$	$\delta =4$	$\delta =5$	$\delta =6$	$\delta =7$
A	$0$	$0$	$0$	$0$	$``{\mathbb{Z}}_2"$	$0$	${\mathbb{Z}}_2$	$0$	${\mathbb{Z}}_2$	$0$	${\mathbb{Z}}_2$	$0$
AIII	$1$	$0$	$0$	$1$	$0$	${\mathbb{Z}}_2$	$0$	${\mathbb{Z}}_2$	$0$	${\mathbb{Z}}_2$	$0$	${\mathbb{Z}}_2$
AI	$0$	$1$	$0$	$0$	$``{\mathbb{Z}}_2"$	$0$	$0$	$0$	${\mathbb{Z}}_2$	$0$	${\mathbb{Z}}_2$	${\mathbb{Z}}_2$
BDI	$1$	$1$	$1$	$1$	${\mathbb{Z}}_2$	${\mathbb{Z}}_2$	$0$	$0$	$0$	${\mathbb{Z}}_2$	$0$	${\mathbb{Z}}_2$
D	$2$	$0$	$1$	$0$	${\mathbb{Z}}_2$	${\mathbb{Z}}_2$	${\mathbb{Z}}_2$	$0$	$0$	$0$	${\mathbb{Z}}_2$	$0$
DIII	$3$	$-1$	$1$	$1$	$0$	${\mathbb{Z}}_2$	${\mathbb{Z}}_2$	${\mathbb{Z}}_2$	$0$	$0$	$0$	${\mathbb{Z}}_2$
AII	$4$	$-1$	$0$	$0$	$``{\mathbb{Z}}_2"$	$0$	${\mathbb{Z}}_2$	${\mathbb{Z}}_2$	${\mathbb{Z}}_2$	$0$	$0$	$0$
CII	$5$	$-1$	$-1$	$1$	$0$	${\mathbb{Z}}_2$	$0$	${\mathbb{Z}}_2$	${\mathbb{Z}}_2$	${\mathbb{Z}}_2$	$0$	$0$
C	$6$	$0$	$-1$	$0$	$0$	$0$	${\mathbb{Z}}_2$	$0$	${\mathbb{Z}}_2$	${\mathbb{Z}}_2$	${\mathbb{Z}}_2$	$0$
CI	$7$	$1$	$-1$	$1$	$0$	$0$	$0$	${\mathbb{Z}}_2$	$0$	${\mathbb{Z}}_2$	${\mathbb{Z}}_2$	${\mathbb{Z}}_2$

Equations

$$H_{\mathrm{SSH}}(k)=\sin k{\sigma }_x+(t+\cos k){\sigma }_z,$$

$$I{H}_{\mathrm{SSH}}(k)I^{-1}$$

$$S{H}_{\mathrm{SSH}}(k)S^{-1}$$

$${\mathcal{H}}_{(1,1)}=k{\sigma }_x+{\phi }_1(x){\sigma }_z,$$

$${\mathcal{H}}_{(1,2)}=k{\sigma }_x+{\phi }_1(x){\sigma }_z+{\phi }_2(x){\sigma }_y,$$

$$\rho (x)$$

$$\displaystyle=\frac{1}{2\pi}\partial_{x}\tan^{-1}\Big{(}\frac{\varphi_{2}}{\varphi_{1}}\Big{)}.$$

$$H(k)$$

$$=\sin k{\tau }_x{\sigma }_x+(t+\cos k){\tau }_x{\sigma }_z.$$

$$\mathcal{H}=k{\tau }_x{\sigma }_x+{\phi }_1(x){\tau }_x{\sigma }_z,$$

$${\mathcal{H}}_{(1,3)}=$$

$$+$$

$${\mathcal{H}}_{(1,4)}=$$

$$+{\phi }_3(x){\tau }_y{\sigma }_0+{\phi }_4(x){\tau }_z{\sigma }_0,$$

$$\mathcal{I}\equiv N_{-}-N_{+}.$$

$${\mathcal{H}}_0(\mathbf{k})=\sum _{i=1}^d{k}_i{\gamma }_i+\lambda \mathcal{M},$$

$$\mathcal{H}(\mathbf{k},\mathbf{r})=\sum _{i=1}^d{k}_i{\gamma }_i+\lambda (\mathbf{r})\mathcal{M}+\sum _{a=1}^{2n-d}{m}_{a,\mathbf{r}}{M}_a,$$

$$I\mathcal{H}(\mathbf{k},\mathbf{r})I^{-1}=\mathcal{H}(-\mathbf{k},-\mathbf{r}),$$

$$\mathcal{H}(\mathbf{k},\mathbf{r})=\sum _{i=1}^d{k}_i{\gamma }_i+\lambda (\mathbf{r})\mathcal{M}+\sum _{a=1}^d{m}_{a,\mathbf{r}}{M}_a.$$

$$h(\mathbf{k},\mathbf{r})=P_{+}(\mathbf{r})\mathcal{H}(\mathbf{k},\mathbf{r})P_{+}(\mathbf{r}),$$

$$h(\mathbf{k},\mathbf{r})={\mathbf{k}}_S\cdot \tilde{\mathbf{\gamma }}+\sum _{a=1}^d{m}_{a,\mathbf{r}}{\tilde{M}}_a,$$

$$\mathbf{\phi }(\mathbf{r})=\frac{{\vec{m}}_{\mathbf{r}}}{|{\vec{m}}_{\mathbf{r}}|},$$

$$w_1$$

$$w_2$$

$$\mathbf{\phi }(\mathbf{r})=-\mathbf{\phi }(-\mathbf{r}).$$

$${\mathcal{H}}_{(d,d+1)}=-\sum _{i=1}^{d}i{\gamma }_i{\partial }_i+\sum _{a=1}^{d+1}{m}_{a,\mathbf{r}}{M}_a,$$

$$J_d^{\mu }=$$

$$\times {\Phi }_{a_1}{\partial }_{{\mu }_1}{\Phi }_{a_2}{\partial }_{{\mu }_2}\dots {\partial }_{{\mu }_d}{\Phi }_{a_{d+1}},$$

$$j_{d-1}^{\mu }=$$

$$\simeq$$

$$\times {\phi }_{a_1}{\partial }_{{\mu }_1}{\phi }_{a_2}{\partial }_{{\mu }_2}\dots {\partial }_{{\mu }_{d-1}}{\phi }_{a_d},$$

$$\delta ={\int }_{S^{d-1}}{j}_{d-1}^0=-(-1{)}^{\frac{d(d+1)}{2}}{w}_{d-1}.$$

$$j_0^0$$

$$j_1^0$$

$$j_2^0$$

Full text

I Introduction

Understanding the role of crystalline symmetries in the topological properties of materials is a central topic in recent studies of condensed matter. Soon after the discovery of topological insulators protected by time-reversal symmetry Kane and Mele (2005); Hasan and Kane (2010), topological insulators (TIs) and topological superconductors (TSCs) have been systematically classified into 10-fold Altland-Zirnbauer (AZ) classes Altland and Zirnbauer (1997), based on non-spatial symmetries including time-reversal, particle-hole, and chiral symmetries Schnyder et al. (2008); Kitaev (2009); Ryu et al. (2010). In this classification, TIs are characterized by nontrivial topological invariants and have a Wannier obstruction that is stable against adding trivial bands. Also $d$-dimensional ($d$D) TIs possess stable gapless states at the $(d-1)$D boundaries, which is called the bulk-boundary correspondence.

On the other hand, the topological properties of crystalline solids can also be protected by space-group symmetries, which are realized in topological crystalline insulators (TCIs) Fu (2011); Hsieh et al. (2012). Compared to the TIs and TSCs classified based on 10-fold AZ classes, TCIs have extremely rich structures because of the large variety of crystalline symmetries forming $230$ space groups and $1651$ magnetic space groups Kruthoff et al. (2017); Bradlyn et al. (2017); Po et al. (2017); Watanabe et al. (2018); Song et al. (2018a); Thorngren and Else (2018); Shiozaki et al. (2018a). For instance, recent studies of TCIs have revealed two intriguing properties that are absent in the 10-fold classification of TIs and TSCs, that is, the fragile topology and the higher-order bulk-boundary correspondence.

Fragile topology is one interesting characteristics of some TCIs discovered in recent theoretical studies Po et al. (2018); Cano et al. (2018); Bradlyn et al. (2019); Bouhon et al. (2018); Wieder and Bernevig (2018); Po et al. (2019); Else et al. (2019); Ahn et al. (2019); Liu et al. (2019a); Song et al. (2019); Kooi et al. (2019); Wieder et al. (2019). In general, TCIs are not adiabatically deformable to atomic insulators without breaking symmetries, that is, they have an obstruction to constructing exponentially localized symmetric Wannier states Thonhauser and Vanderbilt (2006); Soluyanov and Vanderbilt (2011); Marzari et al. (2012). However, the Wannier obstruction of fragile TIs can be resolved by adding appropriately chosen trivial bands. This property of fragile TIs is clearly distinct from the properties of TCIs with stable topology whose band topology is not affected by additional trivial bands. Because of this, TCIs with stable topology can be classified by $K$-theory Freed and Moore (2013); Shiozaki and Sato (2014); Shiozaki et al. (2016, 2017); Kruthoff et al. (2017); Stehouwer et al. (2018), where the equivalence class of insulators is examined when adding trivial bands are allowed. This, at the same time, implies that $K$-theory is not an appropriate tool to study fragile TIs. In general, two different approaches can be considered to examine the Wannier obstruction of fragile TIs. One is to investigate the symmetry representation of exponentially localized Wannier states Po et al. (2018, 2019); Else et al. (2019); Ahn et al. (2019); Liu et al. (2019a); Kooi et al. (2019); Wieder et al. (2019), and the other is to study the winding pattern of the Wilson loop spectrum Cano et al. (2018); Bradlyn et al. (2019); Bouhon et al. (2018); Wieder and Bernevig (2018); Ahn et al. (2019); Kooi et al. (2019); Wieder et al. (2019). While both approaches have been applied successfully to characterize the fragile topology of some specific models, the general characterization of fragile TIs based on them has not been achieved yet.

Another intriguing property of some TCIs is that they do not follow the conventional bulk-boundary correspondence. Recently, it has been found that some TCIs have gapless states only in some subspace of the boundary such as hinges or corners Benalcazar et al. (2017a, b); Schindler et al. (2018a); Langbehn et al. (2017); Song et al. (2017a); Fang and Fu (2017); Ezawa (2018); Khalaf et al. (2018); Khalaf (2018); van Miert and Ortix (2018); Geier et al. (2018); Matsugatani and Watanabe (2018); Wang et al. (2018); Schindler et al. (2018b); Ahn and Yang (2019); Călugăru et al. (2019); Trifunovic and Brouwer (2019); Liu and Wakabayashi (2017); Liu et al. (2019b); Cornfeld and Chapman (2019); Kooi et al. (2019); Song et al. (2018b); Okuma et al. (2019); Shiozaki et al. (2018b); Shiozaki (2019); Benalcazar et al. (2019); Schindler et al. (2019). In general, $d$D topological phases with gapless states in $(d-k)$D boundaries are called $d$D $k$th-order topological phases. In particular, when $1<k\leq d$, they are called higher-order topological phases and their bulk-boundary correspondence has recently been well-established Khalaf et al. (2018); Khalaf (2018); Geier et al. (2018); Trifunovic and Brouwer (2019).

On the other hand, the nature of boundary states of fragile TIs is not clearly understood yet. It is generally expected that the bulk-boundary correspondence of fragile TIs belongs to a rather trivial category. This is because fragile TIs cannot have a peculiar spectral feature at the boundary distinct from that of atomic insulators, since the boundary energy spectrum is insensitive to adding trivial bands below the bulk gap. This is also consistent with the fact that topological phases without gapless boundary states are either fragile topological or atomic Liu et al. (2019a); Trifunovic and Brouwer (2019). Nevertheless, some characteristic boundary features of fragile TIs have been reported Wieder and Bernevig (2018); Wang et al. (2018). For instance, it was shown that, in a two-dimensional (2D) fragile TI protected by inversion symmetry, fractional corner charges appear at the boundary when the geometry of the system preserves inversion symmetry Wieder and Bernevig (2018); Wang et al. (2018). Moreover, many atomic insulators protected by rotation symmetry were shown to host fractional corner charges in $2$D, which is similar to the case of $2$D second-order TIs Benalcazar et al. (2017a, b); Song et al. (2017a); Benalcazar et al. (2019); van Miert and Ortix (2018); Schindler et al. (2019). The origin of corner charges in fragile topological and atomic insulators has been attributed to the filling anomaly Benalcazar et al. (2019), which denotes an obstruction to fulfill the electron filling for charge neutrality when the relevant crystalline symmetry is preserved Wieder and Bernevig (2018); Song et al. (2017a); Wang et al. (2018); Benalcazar et al. (2019); Schindler et al. (2019). To understand the generic boundary properties of fragile TIs and atomic insulators, it is highly desirable to establish their bulk-boundary correspondence beyond that of stable topological insulators.

In this work, we fully characterize the bulk and boundary properties of fragile TIs protected by inversion symmetry. First, we establish the bulk-boundary correspondence of fragile TIs protected by inversion symmetry. For this purpose, we generalize the notion of the $d$D $k$th-order TIs to the cases with $k>d$, and show that the $d$D $(d+1)$th-order TIs have fragile topology and exhibit the filling anomaly with the associated boundary charge accumulation. On the other hand, $d$D $k$th-order TIs with $k>d+1$ have a featureless boundary even though they have fragile topology. Second, we propose an efficient scheme to diagnose fragile band topology by using the symmetry data in momentum space. Using this method, we show that among all the possible parity configurations of inversion-symmetric insulators, at least 17$\%$ of them have fragile topology in $2$D while only 3$\%$ are fragile TIs in $3$D. Finally, we show that the fragile topology of $d$D $k$th-order TI with $k>d$ is manifested in the relative winding of the Wilson loop and nested Wilson loop spectra. Here, we mainly focus on the class A of AZ classification that does not have time-reversal, particle-hole, chiral symmetries. However, we anticipate that the idea we propose here can be generalized to other AZ classes. In particular, the Dirac Hamiltonians of generalized higher-order topological phases are discussed in Appendix E.

To generalize the notion of the higher-order topology to the cases with $k>d$, we use the Dirac Hamiltonian approach that has been developed recently to understand the higher-order topology with $1<k\leq d$ Khalaf (2018). In Ref. Khalaf, 2018, it is shown that when a Dirac Hamiltonian has $(k-1)$ mass terms odd under inversion in addition to the mass term even under inversion, the Dirac Hamiltonian describes a $k$th-order topological phase. The mass term even under inversion is the bulk mass term $\mathcal{M}$ while the other mass terms odd under inversion are boundary mass terms $m_{a,\bm{r}}$ ($a=1,\dots,k-1$). Then, the energy gap at the boundary $E_{g}=2\sqrt{\sum_{a=1}^{k-1}m^{2}_{a,\bm{r}}}$ vanishes in the $(d-k)$D subspace of the boundary, because the boundary mass terms vanish there due to the constraint $m_{a,\bm{r}}=-m_{a,-\bm{r}}$ imposed by inversion symmetry.

Using the Dirac Hamiltonian formalism, it is straightforward to extend the notion of higher-order topology to the cases with $k>d$ Călugăru et al. (2019); Okuma et al. (2019); Shiozaki (2019). For instance, $d$D $(d+1)$th-order TIs have $d$ boundary mass terms, so the boundary is fully gapped in general. In this case, even though the boundary is gapped, inversion symmetry requires that the $(d+1)$ mass terms, including $d$ boundary mass terms and the bulk mass term, form a skyrmionlike structure in real space, and have an odd winding number along the boundary. In Refs. Okuma et al. (2019) and Shiozaki (2019), it was shown that 3D fourth-order Dirac Hamiltonian in $\mathbb{R}^{3}$ is equivalent to a bound state at the origin. However, the relevant bulk-boundary correspondence was not discussed. In this paper, we show that the winding number of boundary mass terms gives rise to the filling anomaly using the Goldstone-Wilczek formula Goldstone and Wilczek (1981). A related idea has already been proposed for $d=2$ in Ref. Wieder and Bernevig, 2018 but our theory is generally valid in arbitrary spatial dimensions. We also clarify the spatial location of extra charge accumulation or depletion. According to the Goldstone-Wilczek formula, since the extra charge density is determined by the mass winding density, the extra charges are localized in the region where the mass terms change rapidly. Since the bulk mass term changes rapidly at the boundary, the induced charges are naturally located at the boundary. The actual location of the extra charge is the position on the boundary where the boundary mass terms vary abruptly. Because of this reason, when the geometry of the system has sharp corners, the extra charges appear there. On the other hand, when the geometry of the system has only smooth boundary, the extra charges are spread over the whole boundary, as confirmed by numerical calculations.

Under the periodic boundary condition in which the fragile band topology is defined, $d$D $(d+1)$th-order TIs can be constructed in the following way. We start from an insulator with $2^{d-1}$ occupied bands in which all the occupied bands have positive parities at all inversion-invariant momenta. Let us call such an insulator as a trivial-parity insulator. Then, replacing all the positive-parity states at an odd number of inversion-invariant momenta by negative-parity states gives a $d$D $(d+1)$th-order TIs. In particular, when a $d$D $(d+1)$th-order TI has negative-parity states only at a single inversion-invariant momentum, the transition between the trivial-parity insulator and the $d$D $(d+1)$th-order TI can naturally be described by a Dirac Hamiltonian. In this construction, $d$D $k$th-order TIs with $k>d$ have the fragile band topology. In particular, a minimal fragile TI with the filling anomaly can be considered as a $d$D $(d+1)$th-order TI. One can construct a $d$D $k$th-order TI with $k>(d+1)$ by stacking $d$D $(d+1)$th-order TIs. In general, a $d$D $k$th-order TI with $k>(d+1)$ does not have any nontrivial feature at the boundary although it has fragile topology. Nevertheless, we show that the trivial-parity insulator and $d$D $k$th-order TI can be distinguished by the net parity $\mathcal{I}=N_{-}-N_{+}$ Thorngren and Else (2018); Song et al. (2017b); Huang et al. (2017) of a finite-size system with open boundary, defined by the difference between the number of occupied states with odd parity $N_{-}$ and the number of occupied states with even parity $N_{+}$. In fact, the net parity $\mathcal{I}$ defined for the given geometry of the finite-size system is identical to one of bulk topological invariants for inversion-symmetric fragile TIs and atomic insulators. This is consistent with the layer construction proposed in Refs. Trifunovic and Brouwer, 2019; Song et al., 2017b; Huang et al., 2017; Shiozaki et al., 2018b; Song et al., 2018b.

The properties of fragile TIs discussed above are rigorously confirmed by analyzing the symmetry data, that is, the inversion parities at inversion-invariant momenta in the Brillouin zone. In fact, in recent studies, it is found that parities can be organized in the form of symmetry indicators that distinguish stable topological phases from fragile topological or atomic insulators Po et al. (2017). However, this method cannot analyze fragile TIs. Here, we establish a simple criterion that can be used to distinguish fragile TIs from atomic insulators efficiently. Using this method, we estimate the number of the parity configurations corresponding to inversion-symmetric stable TIs, fragile TIs, and atomic insulators, respectively, in both $2$D and $3$D.

To investigate the fragile topology further, we systematically examine how many trivial bands are required to turn a fragile TI to an atomic insulator. Using our criterion for fragile topology, we show that for the $d$D $(d+1)$th-order TI which has $2^{d-1}$ occupied bands, $(2^{d-1}-1)$ bands are needed to trivialize the fragile topology. We explain why $(2^{d-1}-1)$ bands are necessary for trivialization by using the Wilson loop method. The idea is that all the spectral winding of Wilson loops Yu et al. (2011); Alexandradinata et al. (2014, 2016) and (generalized) nested Wilson loops Benalcazar et al. (2017a, b) have to be unwound to completely resolve the Wannier obstruction in fragile topological phases. Explicitly, we show that the total number of bands needed to trivialize the $d$D $(d+1)$th-order TIs can be decomposed into $2^{d-1}-1=2^{d-2}+2^{d-3}+\cdots+1$, where $2^{d-2-l}$ is the number of bands needed to remove the spectral winding of the $l$th nested Wilson loop ($l=0,1,2,\dots$). Here, the zeroth nested Wilson loop denotes the conventional Wilson loop. We also generalize this idea to $d$D $k$th-order TIs with $k>(d+1)$ which also have fragile topology.

The rest of this paper is organized as follows. We start Sec. II with the Su-Schrieffer-Heeger (SSH) model in order to introduce the concepts of $d$D $(d+1)$th-order TIs and its higher-order generalizations. Although there is no fragile TI in one dimension because there is no Wannier obstruction in one dimension Marzari et al. (2012); Alexandradinata et al. (2014); Read (2017); Alexandradinata et al. (2019), all the ideas related with the bulk-boundary correspondence can be generalized straightforwardly to higher dimensions. In Sec. III, we introduce the Dirac Hamiltonians of $d$D $k$th-order TIs for arbitrary $d$ and $k$. Then, we show that the $d$ boundary mass terms of $d$D $(d+1)$th-order TI have a nontrivial winding number on the boundary in the presence of inversion symmetry. The winding number of boundary mass terms is directly connected to the induced current which dictates the nontrivial charge accumulation at the boundary. In Sec. IV, we construct a mapping between a parity configuration in the Brillouin zone and the inversion-symmetric Wannier states. With the mapping, we obtain $(2^{d}+1)$ invariants for the class A inversion-symmetric insulators in $d$ dimensions. From this, we develop an efficient method for diagnosing the atomic, fragile topological, and stable topological phases and classify all the possible parity configurations in two and three dimensions. Also, we introduce the net parity of a finite-size system with open boundary and the condition for the presence of the filling anomaly. In Sec. V, we demonstrate our theoretical prediction through numerical calculations of tight-binding models with open boundary condition. In Sec. VI, we show that fragile topological phases protected by inversion symmetry are characterized by the relative winding Alexandradinata et al. (2014) in the (nested) Wilson loop spectrum, which is fragile against adding trivial bands. Using the $3$D fourth-order TI as a concrete example, we count the number of trivial bands needed to unwind the relative winding in the Wilson loop and the nested Wilson loop spectra and show that this matches with the result expected from our diagnosis method. Then, we generalize this counting to the case of $d$D $k$th-order TIs. In Sec. VII, we study the general pumping process of inversion-symmetric insulators. Specifically, we construct $(d+1)$D $(k-1)$th-order TI by a pumping of $d$D $k$th-order TI, which generalizes the result in Ref. Wieder and Bernevig, 2018 for $d=2$ and $k=3$. Finally, we conclude in Sec. VIII.

II Symmetry data and generalized higher-order topology in one-dimension

Here we introduce the notion of the generalized higher-order topology in one-dimensional (1D) systems by considering the Su-Schrieffer-Heeger (SSH) model Su et al. (1979) and its multiple copies. Let us note that the meaning of “topological phases” is a bit subtle in 1D, since only atomic and obstructed atomic insulators Bradlyn et al. (2017) exist in 1D Marzari et al. (2012); Alexandradinata et al. (2014); Read (2017); Alexandradinata et al. (2019), in contrast to higher dimensional systems where topological phases mean the phases with either fragile or stable band topology. For convenience, however, here we use the term “topological phases” to denote the obstructed atomic insulators in $1$D systems.

II.1 Second-order topology in 1D

Let us begin with the $1$D SSH model Hamiltonian with inversion symmetry,

[TABLE]

which is invariant under inversion symmetry $I=\sigma_{z}$ as well as chiral symmetry $S=\sigma_{y}$, that is,

[TABLE]

The energy spectrum is given by $E(k)=\pm\sqrt{\sin^{2}k+(t+\cos k)^{2}}$, which is gapped when $|t|\neq 1$. When $t=1$ ($t=-1$), the band gap closes at $k=\pi$ ($k=0$).

The insulating phases are distinguished by a topological invariant, which is nothing but the winding number of the two-component vector $(t+\cos k,\sin k)$ around the origin $(0,0)$, so that the insulating phase is topological when $|t|<1$ while it is trivial when $|t|>1$.

The topological and trivial insulators can also be distinguished by the parity configuration of occupied states at inversion-invariant momenta. In the trivial phase, the parity of the occupied band at $k=0$ ($p_{0}$) and that at $k=\pi$ ($p_{\pi}$) is $(p_{0},p_{\pi})=(-,-)$ when $t>1$ while $(p_{0},p_{\pi})=(+,+)$ when $t<-1$. In the topological phase with $|t|<1$, $(p_{0},p_{\pi})=(-,+)$. Table 1 shows the parity configurations of atomic insulators with inversion symmetry, which shows that the parity configuration of the topological insulator is equivalent to that of the atomic insulator with a $p$ orbital at the Wyckoff position $B$. On the other hand, the parity configuration of the trivial insulator is equivalent to that of the atomic insulator with a $s$ ($p$) orbital at the Wyckoff position $A$ when $t<-1$ ($t>1$) [see Fig. 1(b)].

To examine the surface states, let us construct a heterostructure composed of a topological insulator in the region with $|x|<L/2$ and trivial insulators on the sides with $|x|>L/2$. Such a heterostructure can also be described by the Hamiltonian in Eq. (1) when $t$ depends on the spatial coordinate $x$. For simplicity, we assume that $-1<t(x)<1$ for $|x|<L/2$ and $t(x)<-1$ for $|x|>L/2$. Since the band gap closes at $k=0$ at the critical point with $t=-1$, one can consider the following low-energy Dirac Hamiltonian:

[TABLE]

where $\varphi_{1}(x)\equiv 1+t(x)$ that plays the role of a bulk mass term. Then, $\varphi_{1}(x)>0$ for $|x|<L/2$ and $\varphi_{1}(x)<0$ for $|x|>L/2$. $\varphi_{1}(x)=\varphi_{1}(-x)$ due to inversion symmetry. Let us note that the symbol $\mathcal{H}_{(d,k)}$ indicates the $d$D Dirac Hamiltonian with $k$ anticommuting mass terms.

When the length $L$ is sufficiently large, there are two zero-energy states at two edges (or two domain walls), where the bulk mass term $\varphi_{1}(x)$ changes its sign [see Fig. 2(c)]. Because of inversion symmetry, the two edge states should be either occupied or unoccupied simultaneously. Thus, an inversion-symmetric ground state cannot satisfy the exact half-filling condition, which is described in Figs. 1(c) and 1(d). This phenomenon is known as the filling anomaly Benalcazar et al. (2019).

Let us now break chiral symmetry by adding an inversion-symmetric second mass term $\varphi_{2}(x)$ to $\mathcal{H}_{(1,1)}$, which gives

[TABLE]

where $\varphi_{2}(x)=-\varphi_{2}(-x)$. We assume that $\varphi_{2}(x)$ is finite only around the domain wall and call this type of mass terms as “boundary” mass terms. The Dirac Hamiltonian with $(k-1)$ boundary mass terms describes a $k$th-order topological insulator. In this respect, $\mathcal{H}_{(1,2)}$ means the Dirac Hamiltonian of $1$D second-order TI. We generalize the notion of $d$D $k$th-order TI for arbitrary $k$ in Sec. III.1.

The boundary mass term $\varphi_{2}$ makes the energy of the edges states shifted from zero. However, since inversion symmetry is maintained, the two edge states are still degenerate, and thus the filling anomaly remains. In particular, the profile of $\varphi_{2}$ depicted in Fig. 2(a) lowers the energy of the edge states. When the edge states are occupied simultaneously, the system has one extra charge relative to the half-filling, and the resulting phase is said to have the filling anomaly with $\delta=1$ [see Fig. 2(d)].

The extra charge density accumulated at a domain wall shown in Fig. 2(b) is given by the Goldstone-Wilczek formula Goldstone and Wilczek (1981)

[TABLE]

It is worth noting that Eq. (II.1) shows that the extra charge is accumulated mostly at the domain wall, even when the edge spectrum merges into the bulk spectrum. This is because the bulk mass term $\varphi_{1}$ varies only near the domain wall where the accumulated charge density is located. This property is shared by any $d$D $(d+1)$th-order inversion-symmetric TIs as shown in Sec. III.3.

The total amount of the extra charge $Q=\int dx\rho(x)$ is given by the winding number of $(\varphi_{1},\varphi_{2})$. Note that we can compactify $R^{1}$ to $S^{1}$ since $(\varphi_{1},\varphi_{2})\rightarrow(-1,0)$ at $x\rightarrow\pm\infty$, and thus the integer winding number can be defined. The boundary conditions on $\varphi_{1}$ and $\varphi_{2}$ imposed by inversion symmetry constrain $Q$ to be 1. This result is consistent with the picture based on the shifting of the energy level of the edge states by adding $\varphi_{2}$. When the bulk mass term $\varphi_{1}(x)$ changes its sign near each domain wall, \tan^{-1}\Big{(}\frac{\varphi_{2}}{\varphi_{1}}\Big{)} is abruptly changed by $\pi$. Therefore, a half-integral charge should be accumulated at each domain wall.

II.2 Third- and fourth-order topology

The idea of higher-order topology can be generalized further even to the cases with $k>(d+1)$. To illustrate the idea, let us first consider two copies of the SSH model. The relevant Hamiltonian is obtained by doubling $H_{\rm SSH}(k)$ in an inversion-symmetric way as

[TABLE]

The corresponding Dirac Hamiltonian linearized near $k=0$ is

[TABLE]

where $\varphi_{1}(x)>0$ for $|x|<L/2$ and $\varphi_{1}(x)<0$ for $|x|>L/2$, which again describes a heterostructure with two domain walls at $x=\pm L/2$. The Hamiltonian $\mathcal{H}$ is invariant under inversion symmetry $I=\tau_{x}\sigma_{z}$ and chiral symmetry $S=\tau_{z}\sigma_{0}$. When chiral symmetry $S=\tau_{z}\sigma_{0}$ is preserved, two boundary mass terms $\varphi_{2}(x)\tau_{x}\sigma_{y}+\varphi_{3}(x)\tau_{y}\sigma_{0}$ are allowed. Adding them to $\mathcal{H}$ gives

[TABLE]

where $\varphi_{2,3}(x)=-\varphi_{2,3}(-x)$. The resulting Dirac Hamiltonian describes a $1$D third-order TI belonging to class AIII.

There is one remaining boundary mass term $\varphi_{4}(x)\tau_{z}\sigma_{0}$ that breaks $S$ but keeps $I$. By adding this term to $\mathcal{H}_{(1,3)}$, we obtain the following Dirac Hamiltonian:

[TABLE]

which describes a $1$D fourth-order TI belonging to class A.

For both $\mathcal{H}_{(1,3)}$ and $\mathcal{H}_{(1,4)}$, there are neither the filling anomaly nor zero-energy states. However, the relative topology between the trivial and topological phases can still be revealed by real-space wave function.

To demonstrate the idea, let us consider a finite-size doubled SSH model composed of $N_{\rm cell}$ unit cells in which the total number of energy eigenstates is $4N_{\rm cell}$. In the case of a trivial insulator, two $s$ orbitals are located at the center of each unit cell, whereas, in the case of a topological insulator, two $p$ orbitals are located at a boundary of each unit cell. We first consider a trivial insulator with an even integer $N_{\rm cell}$ shown in Fig. 3(a). Since a unit cell is related with another unit cell via inversion symmetry, their bonding and antibonding combinations generate equal number of even- and odd-parity states. That is, $N_{+}=N_{-}=N_{\rm cell}$ where $N_{+}$ ($N_{-}$) denotes the number of occupied bands with even (odd) parity. Similarly, when $N_{\rm cell}$ is odd as shown in Fig. 3(c), all the unit cells except the one at the center form inversion-symmetric pairs, and contribute equally to $N_{+}$ and $N_{-}$. However, the central unit cell containing two $s$ orbitals at the inversion center contributes 2 only to $N_{+}$, which leads to $N_{+}=N_{\rm cell}+1$ and $N_{-}=N_{\rm cell}-1$.

Similar analysis can be done for topological insulators. Since there is no filling anomaly, a finite-size doubled SSH model can have inversion-symmetric atomic configurations as shown in Figs. 3(b) and 3(d). In contrast to the trivial insulator in which orbitals are located at the center of each unit cell, the orbitals of a topological insulator are located at the edges of each unit cell. When $N_{\rm cell}$ is even, $p$ orbitals at an edge are related with another $p$ orbitals via inversion symmetry, and their bonding and antibonding combinations contribute equally to $N_{+}$ and $N_{-}$. Since the remaining two $p$ orbitals at the inversion center contribute 2 only to $N_{-}$, we have $N_{+}=N_{\rm cell}-1$ and $N_{-}=N_{\rm cell}+1$. On the other hand, when $N_{\rm cell}$ is odd, all $p$ orbitals are related by inversion symmetry and thus $N_{+}=N_{-}=N_{\rm cell}$. The dependence of $N_{\pm}$ on the parity of $N_{\rm cell}$ is summarized in Table 2.

The four distinct cases considered above can be distinguished in terms of the net parity $\mathcal{I}$ Thorngren and Else (2018); Song et al. (2017b); Huang et al. (2017) defined as (with an additional $-1$ in our convention)

[TABLE]

Since a pair of orbitals related by inversion symmetry contribute equally to $N_{+}$ and $N_{-}$, they do not contribute to $\mathcal{I}$. Thus, it is determined by the difference in the number of $p$ and $s$ orbitals localized at the inversion center of the finite-size system. When $N_{\rm cell}$ is even, one can easily find $\mathcal{I}=0$ ($\mathcal{I}=2$) for a trivial (topological) insulator. On the other hand, when $N_{\rm cell}$ is odd, we obtain $\mathcal{I}=-2$ ($\mathcal{I}=0$) for a trivial (topological) insulator [see Fig. 3].

A similar analysis can also be applied to the $1$D second-order TI in which the filling anomaly exists. The result is summarized in Fig. 4 and Table 3. In fact, $\mathcal{I}$ can serve as a bulk topological invariant for general inversion-symmetric insulators without gapless boundary states, which is discussed in Sec. IV.5.

III Generalized higher-order topology in $d>1$ dimensions

III.1 Dirac Hamiltonian approach

Here we extend the notion of $d$D $k$th-order TI with $k>d$ to general dimensions $d>1$. We focus on the symmetry class A of the AZ symmetry classification. For this purpose, we use the Dirac Hamiltonian describing the transition between two topologically distinct phases. The Dirac Hamiltonians of inversion-symmetric TI/TSCs have been exhaustively studied in Ref. Khalaf, 2018 that we refer to construct the Dirac Hamiltonian with generalized higher-order topology.

In general, a Dirac Hamiltonian is expressed in terms of $2^{n}\times 2^{n}$ $\gamma$ matrices which satisfy the Clifford algebra $\{\gamma_{i},\gamma_{j}\}=2\delta_{ij}$ for $i,j=1,\dots,2n+1$ as

[TABLE]

where $\{\mathcal{M},\gamma_{i}\}=0$ and $\mathcal{M}^{2}=\mathds{1}_{2^{n}}$. Here, $\mathds{1}_{2^{n}}$ denotes the $2^{n}\times 2^{n}$ identical matrix. In Eq. (11), $\mathcal{M}$ is the bulk mass term because the energy eigenvalues are given by $E(\bm{k})=\pm\sqrt{\bm{k}^{2}+\lambda^{2}}$. Changing the sign of the bulk mass term $\mathcal{M}$ flips all the inversion parities at $\bm{k}=0$. We choose $\lambda>0$ ($\lambda<0$) to mean the topological (trivial) insulator. Here, we consider the cases when the inversion symmetry operator is identical to the bulk mass term, i.e., $I=\mathcal{M}$. When we choose $\mathcal{M}=\gamma_{2n+1}$, there remain $(2n-d)$ additional mass terms, $m_{1}\gamma_{d+1}+\cdots+m_{2n-d}\gamma_{2n}$, which anticommute with $\mathcal{M}$, and thus break inversion symmetry when $m_{i=1,\dots,2n-d}$ are constant.

Now we consider a topological insulator with boundaries whose global shape is inversion symmetric, and allow the position $\bm{r}$ dependence of both the bulk mass term $\lambda(\bm{r})\mathcal{M}$ and the boundary mass terms $m_{a,\bm{r}}\gamma_{d+a}$ with $a=1,\dots,2n-d$. The corresponding inversion-symmetric Dirac Hamiltonian including the boundary mass terms is given by

[TABLE]

where $M_{a}=\gamma_{d+a}$ and $\lambda(\bm{r})>0$ [$\lambda(\bm{r})<0$] inside (outside) the topological insulator. $\mathcal{H}(\bm{k},\bm{r})$ satisfies

[TABLE]

in which $m_{a,\bm{r}}=-m_{a,-\bm{r}}$.

The order $k$ of $\mathcal{H}(\bm{k},\bm{r})$ is defined as follows. At the boundary $\bm{r}$ with $\lambda(\bm{r})=0$, the energy gap is given by $2m_{s}(\bm{r})$ where $m_{s}(\bm{r})=\sqrt{\sum_{a=1}^{2n-d}m_{a,\bm{r}}^{2}}$. When $k\leq d$ where $k=2n-d+1$, all boundary mass terms $m_{a,\bm{r}}$ vanish simultaneously on $(d-k)$D domain due to the condition $m_{a,\bm{r}}=-m_{a,-\bm{r}}$ imposed by inversion symmetry. In this case, the order is given by $k$, which is the conventional way of defining the order of a topological insulator. On the other hand, when $k>d$, the boundary states are fully gapped, and thus the order $k$ cannot be well defined. In spite of this, however, we generalize the definition of the order $k$ so that $(k-1)$ is identical to the number of boundary mass terms. Namely, when a $d$D Dirac Hamiltonian has $(k-1)$ boundary mass terms, the relevant TI is a $d$D $k$th-order TI as defined in Sec. I. Let us emphasize once again that a $d$D $k$th-order TI in periodic systems can be well defined by the Dirac Hamiltonian when the relevant trivial insulator is the trivial-parity insulator, so that the sign reversal of the bulk mass term mediates the transition between the trivial-parity insulator and the $d$D $k$th-order TI.

Now we consider the Dirac Hamiltonian for $d$D $k$th-order TI with $k>d$. Since a $d$D $k$th-order TI exists only when $d+k=2l+1$ where $l=1,2,\dots$, the minimal order $k$ of a $d$D TI without gapless boundary states described by a Dirac Hamiltonian is $k=d+1$. The Dirac Hamiltonian of $d$D $(d+1)$th-order TI is

[TABLE]

Then, $d$D $k$th-order TI with $k>(d+1)$ can be obtained by superposing $2^{\frac{k-d-1}{2}}$ copies of $d$D $(d+1)$th-order TI. For instance, in 1D, we have obtained the $1$D fourth-order TI by doubling the $1$D second-order TI. Interestingly, all $d$D $k$th-order TIs with $k>d$ turn out to be a fragile phase as shown in Sec. V. The properties of the Dirac Hamiltonian of $d$D $(d+1)$th-order TI in Eq. (14) are discussed in detail in Secs. III.2 and III.3.

III.2 Homotopy of boundary mass terms

The Dirac Hamiltonian for the boundary theory $h(\bm{k},\bm{r})$ is obtained via the boundary projection of the bulk Hamiltonian Khalaf (2018); Khalaf et al. (2018); Geier et al. (2018); Trifunovic and Brouwer (2019)

[TABLE]

where $P_{+}(\bm{r})=\frac{1}{2}[\mathds{1}-i(\bm{n_{r}}\cdot\bm{\gamma})\mathcal{M}]$ denotes the projection operator to the boundary where $\bm{n}_{\bm{r}}$ is the surface normal vector at the position $\bm{r}$ and $\bm{\gamma}=(\gamma_{1},\dots,\gamma_{d})$. Thus, the boundary theory of $d$D $(d+1)$th-order TI is given by

[TABLE]

where $\bm{k}_{S}=\bm{k}-(\bm{k}\cdot\bm{n_{r}})\bm{n_{r}}$ denotes the boundary momentum and $\tilde{\gamma}_{i}$ and $\tilde{M}_{a}$ are boundary-projected gamma matrices and mass matrices. The details about the boundary projection are summarized in Appendix F. The boundary states of Eq. (16) are gapped everywhere because of $d$ distinct boundary mass terms. Since there is no gapless boundary state, the system seems to be topologically trivial. However, we show that the nontrivial higher-order band topology is manifested in the nontrivial winding structure of the boundary mass terms in real space. With a unit vector $\bm{\varphi}(\bm{r})$ defined in terms of the boundary mass terms $\vec{m}_{\bm{r}}=(m_{1,\bm{r}},\dots,m_{d,\bm{r}})$,

[TABLE]

one can consider a mapping from the boundary position $\bm{r}$ to the unit vector $\bm{\varphi}(\bm{r})$. Note that $\bm{\varphi}(\bm{r})$ is well defined since the boundary gap $m_{s}(\bm{r})=|\vec{m}_{\bm{r}}|$ is non-zero everywhere on the boundary. Also, $\bm{\varphi}(\bm{r})$ satisfies $\bm{\varphi}(\bm{r})\cdot\bm{\varphi}(\bm{r})=1$. Since both the boundary and $\bm{\varphi}(\bm{r})$ are homotopic to the $(d-1)$D sphere $S^{d-1}$, one can define the homotopy group $\pi_{d-1}(S^{d-1})=\mathbb{Z}$. The winding numbers $w_{d}$ relevant to $\pi_{d}(S^{d})=\mathbb{Z}$ are given by

[TABLE]

and similar expressions in higher dimensions. When the system is inversion symmetric, the unit vector $\bm{\varphi}(\bm{r})$ should satisfy

[TABLE]

With this configuration, the winding number is constrained to be an odd integer Outerelo et al. (2009). Let us note that even when a phase transition at the boundary occurs, the winding number modulo $2$ remains invariant due to Eq. (20). This shows that the winding number $w_{d}$ of boundary mass terms characterizes the nontrivial bulk topology. In the following section, Sec. III.3, the connection between the winding number $w_{d}$ and the bulk-boundary correspondence of $d$D $(d+1)$th-order TI is discussed in detail.

III.3 Induced U(1) current and bulk-boundary correspondence

Consider the Dirac Hamiltonian for a $d$D $(d+1)$th-order TI given in Eq. (14),

[TABLE]

where $m_{d+1,\bm{r}}=\lambda(\bm{r})$ and $M_{d+1}=\mathcal{M}$. This can be interpreted as the Hamiltonian describing a massless Dirac fermion coupled to $(d+1)$ scalar fields $M_{a}$. Then, the induced electric current $J_{d}^{\mu}$ is given by

[TABLE]

where $\mu,\mu_{i}=0,1,\dots,d$, $\partial_{0}=\partial_{t}$, $\Phi_{a}=\frac{m_{a,\bm{r}}}{m_{\bm{r}}}$, $m_{\bm{r}}^{2}=\sum_{a=1}^{d+1}m_{a,\bm{r}}^{2}$, ${\rm Area}(S^{d})=\frac{2\pi^{(d+1)/2}}{\gamma(\frac{d+1}{2})}$, and $\gamma(d)$ is the gamma function Abanov and Wiegmann (2000). Here, the sign of $s_{d}=\pm 1$ depends on the representation of gamma matrices where $\text{Tr}[\gamma_{1},\dots,\gamma_{2d+1}]=2^{d}i^{d}$. In our convention for gamma matrices, $s_{d}=-(-1)^{\frac{d(d-1)}{2}}$. The derivation of Eq. (III.3) is given in Appendix G. Let us note that $\partial_{\mu}J_{d}^{\mu}=0$. Moreover, since $\Phi_{a}\rightarrow\{0,\dots,0,-1\}$ as $\bm{r}\rightarrow\infty$, $R^{d}$ can be compactified to $S^{d}$. Then, one can regard $\int_{R^{d}}J_{d}$ as a winding number on $S^{d}$.

To compute the induced current explicitly, let us consider a finite-size $d$D $(d+1)$th-order TI with the boundary between the topological phase and the vacuum (or a trivial insulator). Then, the $d$ boundary mass terms wind along the boundary while the bulk mass term $\lambda(\bm{r})\mathcal{M}$ varies rapidly across the boundary. According to the expression of $J^{\mu}_{d}$ in Eq. (III.3), charges are accumulated in the region where the mass terms vary rapidly, that is, only near the boundary. Note that $J^{\mu}_{1}$ is also used in Sec. II.1 to explain the $1$D second-order TI [see Eq. (II.1)]. Since $J^{\mu}_{d}$ is non-zero only near the boundary, the current density at the boundary $j^{\mu}_{d-1}$ can be determined by integrating $J^{\mu}_{d}$ along $r$ taking into account $\partial_{r}m_{a,\bm{r}}\ll\partial_{r}m_{d+1,\bm{r}}$ for $a=1,\dots,d$ as

[TABLE]

where $\varphi_{a}=\frac{m_{a,\bm{r}}}{|\vec{m}_{\bm{r}}|}$ and $\vec{m}_{\bm{r}}=(m_{1,\bm{r}},\dots,m_{d,\bm{r}})$. One can also derive Eq. (III.3) rigorously in the following way. When the domain of the TI is sufficiently large, the boundary theory of $d$D $(d+1)$th-order TI given in Eq. (16), can be regarded as the Dirac Hamiltonian defined in $R^{d-1}$ with $d$ mass terms. Then, the induced current $J^{\mu}_{d-1}$ is obtained by the boundary projection, and its sign is fixed by comparison with Eq. (III.3). We get $j_{d-1}^{\mu}$ at the $(d-1)$D boundary of $d$D $(d+1)$th-order TI after compactifying $R^{d-1}$ to $S^{d-1}$ (see Appendix G).

Equation (III.3) shows that the total charge of boundary states $\int_{S^{d-1}}j^{0}_{d-1}$, which is nothing but the anomalous charge $\delta$ due to the filling anomaly, is determined by the winding of boundary mass terms as

[TABLE]

We note that the condition $\varphi_{a}(\bm{r})=-\varphi_{a}(-\bm{r})$ imposed by inversion symmetry restricts the winding to be an odd integer. The anomalous charge $\delta$ modulo two characterizes the bulk-boundary correspondence of $d$D $(d+1)$th-order TIs since any inversion-preserving process can change the boundary charge only by even integers unless it closes the bulk gap. The charge densities predicted by $j^{0}_{d-1}$ are as follows:

[TABLE]

[TABLE]

[TABLE]

It is worth noting that the induced currents shown above were also discussed in Refs. Goldstone and Wilczek, 1981; Chamon et al., 2008a, b; Ryu et al., 2009; Yao and Lee, 2010; Fukui, 2019 in various contexts. In all the previous studies, however, all the mass terms are considered as either bulk order parameters or bulk pumping parameters in time direction. In contrast to this, in our case, all $\varphi(\bm{r})$ are defined along the boundary. Thus, the physical implication of the induced currents in our work is different from that in the previous studies. Before closing this section, let us remark on the cases when the boundary projection is not well defined or the boundary is not smooth. In general, the boundary projection is not well-defined when the boundary mass terms are comparable to or larger than the bulk mass term. Nevertheless, as long as the total amount of boundary charge $\delta$ (modulo two) is preserved under any bulk gap-preserving and inversion-symmetric perturbation, our discussion on the bulk-boundary correspondence is still valid. When the boundary is not smooth, it can be treated as an extreme curvature limit of the system with smooth boundary. In Sec. V.1, we numerically show that our predictions of boundary charge are still valid even when the boundary mass terms are comparable to the bulk mass term and/or the boundary is not smooth. In addition to the total amount of the boundary charge, the distribution of the boundary charge density is also consistent with that of the winding number density.

IV Symmetry data and classification of inversion-symmetric insulators

In this section, by using symmetry data, we systematically classify inversion-symmetric insulators into three classes: atomic insulators, fragile topological insulators, and stable topological insulators. In particular, we provide a convenient scheme to diagnose fragile topology, and estimate the number of parity configurations belonging to each class. Also, we discuss the net parity of a finite-size system with open boundary and the condition for the presence of the filling anomaly.

In a $d$D inversion-symmetric system, there are $2^{d}$ distinct Wyckoff positions $W_{i}$ and inversion-invariant momenta $\bm{k}^{{\rm inv}}_{i}$ which are left invariant under inversion symmetry:

[TABLE]

where $\bm{a}$ and $\bm{G}$ are the lattice and reciprocal vectors, respectively, which satisfy $\bm{a}_{i}\cdot\bm{G}_{j}=2\pi\delta_{ij}$. We label $W_{i}$ and $\bm{k}^{{\rm inv}}_{i}$ as follows.

[TABLE]

where $\sigma=0$, $1$. More explicitly,

[TABLE]

and similarly,

[TABLE]

IV.1 Mapping between Wannier states and parity configurations

Any given parity configuration $\mathcal{B}$ can be mapped to linear combinations of inversion-symmetric Wannier states $w_{\xi}(W)$. Here, $w_{\xi}(W)$ denotes a Wannier state with parity $\xi$ localized at the Wyckoff position $W$. At the momentum $\bm{k}^{{\rm inv}}$ in the Brillouin zone, a state with parity $\xi\exp(2i\bm{k}^{{\rm inv}}\cdot W)$ can be induced by $w_{\xi}(W)$. For example, $s$ and $p$ orbitals localized at $W_{1}$ correspond to $w_{+}(W_{1})$ and $w_{-}(W_{1})$, respectively. Let us note that it is sufficient to consider symmetric Wyckoff positions $W_{i}$ only. When a orbital is located at generic point $W\neq-W$, its partner must also be located at $-W$ because of inversion symmetry. Then, the symmetric and antisymmetric combinations of the two states give $s$ and $p$ orbitals at a certain symmetric Wyckoff position, say $W_{i}$. It is worth noting that this $s$ and $p$ orbital pair located at $W_{i}$ can be deformed adiabatically to another Wyckoff position, say $W_{j}$, without breaking inversion symmetry. Such an adiabatic deformation of Wannier states can be referred as “Wannier gauge redundancy.”

For a given Wannier state representation $\bigoplus_{W,\xi}\,\mu_{W,\xi}w_{\xi}(W)$, there is a corresponding parity configuration $\mathcal{B}$ in the Brillouin zone. The equivalence relation between them Song et al. (2017a); Fang and Fu (2017); Cano et al. (2018); Bradlyn et al. (2019); Po et al. (2019); Benalcazar et al. (2019); Wieder and Bernevig (2018); Ahn et al. (2019); Song et al. (2019) is denoted as

[TABLE]

Here, $\mu_{W,\xi}$ is a rational number Po et al. (2017), i.e., $\mu_{W,\xi}\in\mathbb{Q}$, and dictates how many $w_{\xi}(W)$ are superposed in a given atomic configuration. For later convenience, we denote a parity configuration of $N_{\mathcal{B}}$ bands in $d$D that satisfies $n_{-}(\gamma=\bm{k}^{{\rm inv}}_{1})=n_{+}(\bm{k}^{{\rm inv}}\neq\gamma)=N_{\mathcal{B}}$ as $\mathcal{B}_{\gamma}(d,N_{\mathcal{B}})$. Namely, all occupied bands have odd parity at the $\gamma$ point whereas they have even parity at the other $\bm{k}^{{\rm inv}}$ in $\mathcal{B}_{\gamma}(d,N_{\mathcal{B}})$.

The equivalence relation Eq. (34) is not one to one since there is an ambiguity in choosing Wannier states, that is, the Wannier gauge redundancy:

[TABLE]

This follows from the fact that $w_{+1}(W)\oplus w_{-1}(W)$ gives $n_{+}(\bm{k}^{{\rm inv}})=n_{-}(\bm{k}^{{\rm inv}})=1$ for any $W$ and $\bm{k}^{{\rm inv}}$.

According to the recent studies Po et al. (2017); Bradlyn et al. (2017); Po et al. (2018), a parity configuration $\mathcal{B}$ can be classified into three distinct classes depending on the type of $\mu_{W,\xi}\in\mathbb{Q}$. (i) $\mathcal{B}$ is diagnosed to be a stable topological phase when the set of $\mu_{W,\xi}$ includes a fraction. Adding trivial bands corresponds to $\Delta\mu_{W,\xi}\in\mathbb{Z}$, and thus it does not change the fractional part of $\mu_{W}$. (ii) $\mathcal{B}$ is diagnosed to be a fragile topological phase when the set of $\mu_{W,\xi}$ includes at least one negative integer in any fixed Wannier gauge. (iii) $\mathcal{B}$ is diagnosed to be an (obstructed) atomic phase when all $\mu_{W,\xi}$ are integers greater than or equal to [math] in at least one fixed Wannier gauge.

For example, (i) $\mathcal{B}_{\gamma}(2,1)$ is equivalent to $\ominus\frac{1}{2}p(W_{1})\oplus\frac{1}{2}p(W_{2})\oplus\frac{1}{2}p(W_{3})\oplus\frac{1}{2}p(W_{4})$. Because of the fractional coefficients of $p$ orbitals, the phase is diagnosed to be a stable topological phase. In fact, according to the recent symmetry indicator analysis, this phase is characterized by an odd Chern number Fang et al. (2012); Po et al. (2017), which indicates its stable band topology. (ii) $\mathcal{B}_{\gamma}(2,2)$ is equivalent to $\ominus p(W_{1})\oplus p(W_{2})\oplus p(W_{3})\oplus p(W_{4})$. The negative coefficient of $p(W_{1})$ cannot be removed in any Wannier gauge, and thus the phase is diagnosed to be a fragile topological phase. (iii) $\mathcal{B}\Leftrightarrow s(W_{1})\ominus s(W_{2})\oplus 2p(W_{1})$ has a negative coefficient of $s(W_{2})$, and seems to be a fragile topological phase. However, $\mathcal{B}\Leftrightarrow p(W_{1})\oplus p(W_{2})$ up to the Wannier gauge redundancy $s(W_{1})\oplus p(W_{1})=s(W_{2})\oplus p(W_{2})$, thus $\mathcal{B}$ is diagnosed to be an atomic phase.

Now, let us compare the real-space data $\{w_{\xi}(W)\}$ and the momentum-space data $\{n_{\xi}(\bm{k}^{{\rm inv}})\}$. There are $2^{d+1}$ possible types of $w_{\xi}(W)$, counting the number of Wyckoff positions with two different parities. But, only $(2^{d}+1)$ of them are independent due to the Wannier gauge redundancy shown in Eq. (35). Namely, the Wannier gauge redundancy $s(W_{1})\oplus p(W_{1})=s(W_{2})\oplus p(W_{2})=\dots=s(W_{2^{d}})\oplus p(W_{2^{d}})$ give $(2^{d}-1)$ constraint equations. Taking into account this, one can write $(2^{d}+1)$ numbers invariant under the Wannier gauge redundancy as

[TABLE]

where $\nu_{W}=\mu_{W,-}-\mu_{W,+}$. Note that $\mu_{W,-}$ and $\mu_{W,+}$ change by the same amount under the Wannier gauge redundancy, thus, $\nu_{W}$ can serve as a topological invariant. This approach was introduced in Ref. van Miert and Ortix, 2018 for rotation-symmetric ($C_{n=2,3,4,6}$) insulators and $\nu_{W}$ indices were found explicitly. However, the possibility to interpret $\nu_{W}$ indices as the indicators for fragile topological phases was not appreciated. We show that inversion-symmetric insulators can be completely classified in terms of the indices in Eq. (36), which can also be used for the diagnosis of fragile topological phases.

On the other hand, the parity configuration $\mathcal{B}$ is characterized by $2^{d+1}$ integers,

[TABLE]

In an insulating phase, the number of occupied states $n_{\rm occ}$ is fixed to be $n_{\rm occ}=n_{+}(\bm{k}^{{\rm inv}})+n_{-}(\bm{k}^{{\rm inv}})$ for all $\bm{k}^{{\rm inv}}$. This is called the compatibility relation Kruthoff et al. (2017); Po et al. (2017); Bradlyn et al. (2017). Thus, $\mathcal{B}$ can be characterized by $(2^{d}+1)$ independent integers. One canonical choice of these $(2^{d}+1)$ integers is

[TABLE]

including $n_{-}(\bm{k}^{{\rm inv}}_{i})-n_{+}(\bm{k}^{{\rm inv}}_{i})$ with $i=1,\dots,2^{d}$. Observing that $w_{\xi}(W)$ has the parity $\xi\exp(2i\bm{k}^{{\rm inv}}\cdot W)$ at $\bm{k}^{{\rm inv}}$, one can find a mapping between $\{\nu_{W}\}$ and $\{n_{-}(\bm{k}^{{\rm inv}})-n_{+}(\bm{k}^{{\rm inv}})\}$ as

[TABLE]

where the $2^{d}\times 2^{d}$ matrix $A_{d}$ maps the real data $\{\nu_{W}\}$ and the momentum-space data $\{n_{-}(\bm{k}^{{\rm inv}})-n_{+}(\bm{k}^{{\rm inv}})\}$, and its components are given by

[TABLE]

Also, $A_{d}$ satisfies the following relations:

[TABLE]

For example, $A_{d=2}$ is given by

[TABLE]

and the expressions for $\nu_{W_{i=1,2,3,4}}$ are the same as those in Ref. van Miert and Ortix, 2018. From the fact that $n_{\pm}(\bm{k}^{{\rm inv}})$ must be a non-negative integer, the possible sets of $\{\nu_{W}\}$ are restricted. The most fundamental restriction is that 2^{d}\big{|}\sum_{j=1}^{2^{d}}[A_{d}]_{ij}\nu_{W_{j}}\big{|} must be an integer less than or equal to $n_{\rm occ}$. Below, we list some useful restrictions and inequalities derived from this:

[TABLE]

For example,

[TABLE]

cannot be realized when $n_{\rm occ}=1$ in 2D. The first set $\{\frac{1}{2},0,-\frac{1}{2},0\}$ is forbidden due to $2^{d-2}(\nu_{W_{i}}\pm\nu_{W_{j\neq i}})\in\mathbb{Z}$ while the second set $\{-1,0,1,0\}$ is forbidden because of $|\nu_{W_{i}}\pm\nu_{W_{j}}|\leq n_{\rm occ}$. The properties of $A_{d}$ and $\nu_{W}$, and the derivation of the equations from Eq. (39) to Eq. (43) are shown in Appendix B.

IV.2 Diagnosis of fragile topology

Now, we provide an efficient scheme for diagnosing the band topology in terms of $\{n_{\rm occ},\nu_{W_{1}},\dots,\nu_{W_{2^{d}}}\}$. Any given inversion-symmetric phase is an element of $[\Omega]$ where $\Omega=A,~{}F,~{}St$. Here, $[A]$, $[F]$, and $[St]$ denote a set of atomic, fragile topological, stable topological phases, respectively. Each set can be partitioned into mutually exclusive subsets: $[St]=St(St)\cup St(F)\cup St(A)$, $[F]=F(F)\cup F(A)$, and $[A]=A(A)$. The notation $\Omega_{i}(\Omega_{j})$ means that the given inversion-symmetric phase belongs to $\Omega_{i}$, but its parity configuration $\mathcal{B}$ predicts it to be in the set $\Omega_{j}$. We denote $A(A)=A$ for simplicity.

The diagnosis of stable topological phases has been well-established in $2$D and $3$D in terms of symmetry indicators Po et al. (2017); Watanabe et al. (2018). However, an efficient method to diagnose fragile phases has not been established. Here, we propose one efficient scheme. The procedure for diagnosing inversion-symmetric insulators proceeds as follows. (i) $St(St)$ is diagnosed by a criterion $\{\nu_{W_{1}},\dots,\nu_{W_{2^{d}}}\}\notin\mathbb{Z}^{2^{d}}$ where $\mathbb{Z}^{2^{d}}$ indicates the set of $2^{d}$ integers. (ii) A fragile phase can be diagnosed by a criterion $\sum_{W}|\nu_{W}|>n_{\rm occ}$, which can be obtained by the fact that an atomic phase must satisfy $\sum_{W}|\nu_{W}|\leq n_{\rm occ}$. Thus, a phase which satisfies $\sum_{W}|\nu_{W}|>n_{\rm occ}$ belongs to $St(F)\cup F(F)$. (iii) A phase which satisfies $\sum_{W}|\nu_{W}|\leq n_{\rm occ}$ cannot be distinguished from atomic phases, so it belongs to $St(A)\cup F(A)\cup A$. This two-step diagnosis procedure is summarized in the following diagram, and also in Fig. 5.

$\{\nu_{W_{1}},\dots,\nu_{W_{2^{d}}}\}\notin\mathbb{Z}^{2^{d}}$ $\Rightarrow$ $St(St)$

$\Downarrow$

$\sum_{W}|\nu_{W}|>n_{\rm occ}$ $\Rightarrow$ $St(F)\cup F(F)$

$\Downarrow$

$\sum_{W}|\nu_{W}|\leq n_{\rm occ}$ $\Rightarrow$ $St(A)\cup F(A)\cup A$

.

Let us emphasize that our diagnosis procedure goes beyond the diagnosis method proposed previously for stable topological phases Po et al. (2017); Song et al. (2018a) in that our approach can diagnose even fragile topological phases. The criterion for fragile topological and atomic phases is derived in Appendix A. It is worth noting that our classification of inversion-protected fragile TIs is complete. The band topology of class A insulators is fully characterized by inversion parities at high-symmetry points and all the Chern numbers defined in the Brillouin zone Kruthoff et al. (2017); Po et al. (2017); Shiozaki et al. (2018a). Since the Chern numbers are stable topological invariants, fragile TIs and atomic insulators can be fully characterized by using only inversion-parities. This implies that $F(A)$ is an empty set and $[F]=F(F)$.

IV.3 Results of diagnosis

We tabulate how many parity configurations are diagnosed to be atomic, fragile topological, and stable topological phases. We first fix the number of occupied bands $n_{\rm occ}$. Then, we apply the diagnosis procedure to all possible parity configurations, and count $N(A)$, $N(F)$, and $N(St)$, which are the number of parity configurations diagnosed to be atomic, fragile topological, and stable topological phases, respectively. The sum of $N(A)$, $N(F)$, and $N(St)$ must be $N_{\rm tot}$. Here, $N_{\rm tot}$ is the number of all possible parity configurations and is given by $(n_{\rm occ}+1)^{2^{d}}$. Let us note that when the ordering of energy bands is neglected, there can be $0,1,\dots,n_{\rm occ}$ states with odd parity at each inversion-invariant momentum, so there are $(n_{\rm occ}+1)^{2^{d}}$ possible parity configurations.

In 2D, as shown in Table 4, $N(A)$, $N(F)$, and $N(St)$ are given by $\frac{1}{3}N_{\rm tot}$, $\frac{1}{6}N_{\rm tot}$, and $\frac{1}{2}N_{\rm tot}$, respectively, when $n_{\rm occ}$ is sufficiently large. Namely, at least, one-sixth (one-half) of parity configurations are predicted to be fragile (stable) topological phases. On the other hand, in 3D, most of parity configurations are diagnosed to be stable topological phases as shown in Table 5. When $n_{\rm occ}$ is sufficiently large, we find that $N(A)$ is given by $\frac{1}{315}N_{\rm tot}$ but it is difficult to determine the exact values of $N(F)$ and $N(St)$ in the $n_{\rm occ}\rightarrow\infty$ limit. However, from the values of $N(A)$, $N(F)$, and $N(St)$ for various $n_{\rm occ}$, we conjecture that $N(F)=\frac{283}{10080}N_{\rm tot}$ and $N(St)=\frac{31}{32}N_{\rm tot}$, which implies that only about 3$\%$ of parity configurations are diagnosed to be fragile topological phases.

The implications of Tables 4 and 5 are as follows. First, there is no fragile TI with $n_{\rm occ}=1$ in any dimensions Alexandradinata and Höller (2018). This follows from the second inequality in Eq. (43), $|\nu_{W_{i}}\pm\nu_{W_{j}}|\leq n_{\rm occ}$. Second, a layer-stacking of $2$D fragile TIs with negligible interlayer coupling gives rise to fragile TIs in higher dimensions. For instance, in $2$D, there are eight fragile TIs with $n_{\rm occ}=2$. The $\nu_{W}$ indices of these fragile TIs are given by $\{\mp 1,\pm 1,\pm 1,\pm 1\}$ and its permutations. With these representative $\nu_{W}$ indices, the fragile TIs with $n_{\rm occ}=2$ in $3$D can be constructed in the following way. Let us consider $\{\nu_{W}\}=\{-1,1,1,1,0,0,0,0\}$ as an example of indices for a $3$D insulator. These indices correspond to a $3$D insulator obtained by a stacking of 2D layers with $\{\nu_{W}\}=\{-1,1,1,1\}$, which satisfy $\bm{d}\cdot(x,y,z)=0$ (mod 1) where $\bm{d}=(0,0,1)$ and $(x,y,z)$ denotes the position on the $2$D layers. Since there are $14$ independent $2$D layers that satisfy $\bm{d}\cdot(x,y,z)=0$ or $\frac{1}{2}$ (mod $1$) where $\bm{d}=(1,0,0),(0,1,0),(0,0,1),(1,1,0),(1,0,1),(0,1,1),(1,1,1)$, we have $8\times 14=112$ fragile TIs with $n_{\rm occ}=2$ in $3$D. Also, we note that the $1344$ fragile TIs with $n_{\rm occ}=3$ in $3$D can be obtained by adding a trivial band to $112$ fragile TIs with $n_{\rm occ}=2$. There are 16 ways to add a trivial band, which is one of $s(W_{i})$ and $p(W_{i})$ for $i=1,\dots,8$. However, four of them trivialize a 3D fragile TI with $n_{\rm occ}=2$. Accordingly, we have $112\times 12=1344$ fragile TIs with $n_{\rm occ}=3$ in 3D. From the layer-stacking picture, we conclude that the minimal number of occupied bands to have the fragile topology is two in any dimensions. This process for obtaining a higher-dimensional fragile TI by layer stacking $2$D fragile TIs can be described by using a pumping process in momentum space as discussed in Sec. VII. Third, fragile topology can be trivialized under the doubling of the unit cell. For instance, let us consider the doubling of the unit cell in all spatial directions. In this case, all the symmetric Wyckoff positions except $W_{2^{d}}$ form inversion-symmetric pairs at generic positions, thus, the resulting $\nu_{W}$ is composed of $2^{d}$ $\nu_{W_{2^{d}}}$ of the original system, that is, $\nu_{W_{i}}=\nu_{W_{2^{d}}}$ for $i=1,\dots,2^{d}$. If we repeat the same doubling process, the number of occupied bands increases while $\sum_{W}|\nu_{W}|=2^{d}|\nu_{W_{2^{d}}}|$ is fixed. Then, the result of the diagnosis for fragile TI, $\sum_{W}|\nu_{W}|>n_{\rm occ}$, can change to $\sum_{W}|\nu_{W}|\leq n_{\rm occ}$. Finally, when $n_{-}(\bm{k}^{\rm inv}_{i})$ for $i=1,\dots 2^{d}$ are less than $2^{d-1}$, at least one $\nu_{W}$ index among $\nu_{W_{2,3,\dots,2^{d}}}$ is zero, i.e., $\prod_{i=2}^{2^{d}}\nu_{W_{i}}=0$. We have verified that this is true for arbitrary number of occupied bands in $2$D and $3$D, and we expect that this is also true in any dimension $d>3$.

IV.4 Filling anomaly in general dimensions

In this section, we propose a condition for the presence or absence of the filling anomaly, independent of the geometry of finite-size systems. First, we observe that, for a given finite-size system, the number of the Wannier centers at the boundary modulo two, $N_{\partial}$, is invariant under an attaching or a detaching of unit cells to the boundary in an inversion-symmetric way. This is depicted in Fig. 6(a). Second, there are $2^{d}$ building blocks for the finite-size systems in $d$D. Any finite-size inversion-symmetric insulators can be constructed by attaching unit cells at the boundary of these building blocks. For example, in 2D, there are four building blocks composed of $1\times 1$, $2\times 1$, $1\times 2$, and $2\times 2$ unit cells [see Fig. 6(b)]. For the building block made of the $1\times 1$ unit cell, $N_{\partial}$ is given by $\nu_{W_{2}}+\nu_{W_{3}}+\nu_{W_{4}}$ modulo $2$. However, in the case of the building block composed of $1\times 2$ unit cells, $N_{\partial}$ is given by $2\nu_{W_{2}}+\nu_{W_{3}}+2\nu_{W_{4}}$ modulo $2$. Since the filling anomaly exists when $N_{\partial}$ is odd and $N_{\partial}$ depends on the building blocks, the presence or absence of the filling anomaly depends on the type of these building blocks for given $\nu_{W_{2,3,4}}$. Investigating all the building blocks in $d$D, we conclude that $\nu_{W_{i}}$ for $i=2,3,\dots,2^{d}$ must be odd integers (even integers) at the same time to have (not to have) the filling anomaly independent of the type of the building blocks. Physically, this condition is equivalent to the vanishing polarization, $\bm{P}=\sum_{W}W\nu_{W}=(0,0)$ in $2$D. In the rest of this paper, we focus on the cases where the filling anomaly is well defined independent of the type of the building blocks.

The condition for the presence of the filling anomaly puts a constraint to the number of occupied bands. As at least one $\nu_{W}$ index among $\nu_{W_{2,3,\dots,2^{d}}}$ is zero when $n_{-}(\bm{k}^{\rm inv}_{i})$ for $i=1,\dots,2^{d}$ are less than $2^{d-1}$, the filling anomaly occurs only when $n_{\rm occ}\geq 2^{d-1}$. For atomic insulators, the number of occupied bands must be greater than or equal to $(2^{d}-1)$ in order to have all $\nu_{W_{2,3,\dots,2^{d}}}$ as odd integers. Thus, there are $2^{2^{d}-1}$ minimal atomic insulators with the filling anomaly, and they have $\{\nu_{W}\}=\{0,\pm 1,\pm 1,\pm 1\}$ in $2$D and $\{0,\pm 1,\pm 1,\pm 1,\pm 1,\pm 1,\pm 1,\pm 1\}$ in $3$D. Here, all $\pm$ are independent. In contrast to atomic insulators, fragile TIs can have the filling anomaly when $n_{\rm occ}\geq 2^{d-1}$. In other words, all insulators with the filling anomaly are fragile TIs when $2^{d-1}\leq n_{\rm occ}\leq 2^{d}-2$.

IV.5 Net parity in open boundary

In this section, we discuss the net parity Trifunovic and Brouwer (2019); Song et al. (2017b); Huang et al. (2017) of finite-size inversion-symmetric insulators without gapless boundary states, which generalizes the discussion in Sec. II.2. For any inversion-symmetric insulator in $d$D which does not host gapless boundary states, the energy spectra of occupied and unoccupied states are well defined. In this case, the net parity of a finite-size system $\mathcal{I}=N_{-}-N_{+}$ in Eq. (10), characterizes the bulk band topology Trifunovic and Brouwer (2019); Song et al. (2017b); Huang et al. (2017). Here, $N_{\pm}$ denotes the number of occupied states with parity $\pm$. Now, we argue that $\mathcal{I}$ does not change if we add localized states at the boundary. For instance, when an electron is added to the boundary, its inversion partner also has to be added to preserve inversion symmetry as shown in Fig. 7(a). Then, the two electrons form a symmetric and an antisymmetric states. In this process, both $N_{+}$ and $N_{-}$ increase by one, so $\mathcal{I}$ is invariant. The only way to change $\mathcal{I}$ is to add electrons at the inversion center. However, this cannot be done in an inversion-symmetric way. Therefore, $\mathcal{I}$ is the bulk topological invariant of inversion-symmetric systems that does not change as long as the bulk gap is preserved.

As discussed in Sec. II.2, the net parity $\mathcal{I}$ is determined by the difference in the number of $p$ and $s$ orbitals at the center of the system with an open boundary, $\Delta N_{\rm center}$. When the system with an open boundary is inversion symmetric, its center must be one of the Wyckoff position $W^{*}$. Hence, $\Delta N_{\rm center}$ is given by $\nu_{W^{*}}$, and the net parity $\mathcal{I}$ is manifested as one of $\nu_{W}$ indices as shown in Fig. 7(b).

V $d$D $(d+1)$th-order TI as minimal fragile topological insulator with filling anomaly

In this section, we propose a minimal Dirac Hamiltonian for a fragile TI with filling anomaly. We assume that the number of occupied bands is $N_{\mathcal{B}}$, and all the occupied bands at the $\gamma$ point have negative parity while they have positive parity at all the other inversion-invariant momenta. Then, the corresponding parity configuration is $\mathcal{B}_{\gamma}(d,N_{\mathcal{B}})$. Since such a parity configuration can be generated by a band inversion at $\gamma$, starting from a trivial-parity insulator having only positive parity states at inversion-invariant momenta, the state corresponding to $\mathcal{B}_{\gamma}(d,N_{\mathcal{B}})$ can be described by a Dirac Hamiltonian.

The minimal dimension of the Dirac Hamiltonian can be determined by our diagnosis method of fragile topology as follows. First, we compute $2^{d}$ inversion indices $\{\nu_{W_{1}},\dots,\nu_{W_{2^{d}}}\}$ by using Eq. (39). Namely,

[TABLE]

which gives

[TABLE]

To be a fragile topological phase, $\{\nu_{W}\}$ must be a set of integers and satisfy $\sum_{W}|\nu_{W}|>N_{\mathcal{B}}$. These conditions constrain $N_{\mathcal{B}}$ to be an integer multiple of $2^{d-1}$ in any dimensions greater than one. Thus, the minimal number of occupied bands for the fragile phase is $N_{\mathcal{B}}=2^{d-1}$. Importantly, all $\nu_{W_{i}}$ for $i=2,\dots,2^{d}$ are odd integers when $N_{\mathcal{B}}=2^{d-1}$, and thus the filling anomaly exists. This is consistent with the result discussed in Sec. IV.4 where it is shown that the filling anomaly can occur when $n_{\rm occ}\geq 2^{d-1}$. Also, we can find the equivalent representation for Wannier states corresponding to $\mathcal{B}_{\gamma}(d,2^{d-1})$ by using $\nu_{W}=\mu_{W,-}-\mu_{W,+}$ such as

[TABLE]

Let us note three implications of Eqs. (46) and (47). First, the Wannier obstruction can be lifted after adding $(2^{d-1}-1)$ $p(W_{1})$ orbitals. Second, $\mathcal{B}_{\gamma}(d,N_{\mathcal{B}})$ with $N_{\mathcal{B}}<2^{d-1}$ must have stable band topology since $\nu_{W}$ must be a fractional number. Third, an integer multiple of $\mathcal{B}_{\gamma}(d,2^{d-1})$ retains the fragile topology since multiplying Eq. (46) by an integer does not change the relevant inequality, i.e., $\sum_{W}|\nu_{W}|>n_{\rm occ}$.

In Sec. III.1, it is shown that the Dirac Hamiltonian for a class A $d$D $k$th-order TI is expressed by $2^{\frac{d+k-1}{2}}\times 2^{\frac{d+k-1}{2}}$ gamma matrices. Hence, the Dirac Hamiltonian corresponding to $\mathcal{B}_{\gamma}(d,2^{d-1})$ can describe a $d$D $(d+1)$th-order TI. For example, the Dirac Hamiltonian relevant to $\mathcal{B}_{\gamma}(2,2)$ describes a $2$D third-order TI, and similarly, the Dirac Hamiltonian relevant to $\mathcal{B}_{\gamma}(3,4)$ describes a $3$D fourth-order TI. The corresponding parity configurations are shown in Fig. 8. When the number of occupied states is greater than $2^{d-1}$, the corresponding Dirac Hamiltonian describes a $d$D $k$th-order TI with $k>(d+1)$. As discussed above, all $d$D $k$th-order TIs with $k>d$ have fragile band topology.

In Secs. V.1 and V.2, we examine the charge accumulation in the $d$D $(d+1)$th-order TI and establish the bulk-boundary correspondence, which can be viewed as $2$D and $3$D generalizations of the idea discussed in Sec. II.1.

Before we proceed, let us comment on the case where the occupied states with negative parity are distributed over the multiple inversion-invariant momenta. In this case, the minimal number of occupied bands to have fragile topology is two in any dimensions which is smaller than $2^{d-1}$ for $d>2$ (see Sec. IV.3). However, these insulators cannot have the filling anomaly unless the number of occupied bands is greater than or equal to $2^{d-1}$ as discussed in Sec. IV.4. In Appendix D, we discuss how an odd number of band inversions corresponding to the $d$D $(d+1)$th-order Dirac Hamiltonian generates the filling anomaly.

V.1 Tight-binding model

First, we define $4\times 4$ gamma matrices $\gamma_{n=1,...,5}$ and $8\times 8$ gamma matrices $\check{\gamma}_{m=1,...,7}$ in order to construct the tight-binding models for class A $2$D third-order and $3$D fourth-order TIs protected by inversion symmetry. We choose that $\gamma_{n}$ and $\check{\gamma}_{n}$ with an even (odd) $n$ are real (imaginary). For example, $\{\gamma_{1},\dots,\gamma_{5}\}=\{\tau_{x}\sigma_{x},\tau_{x}\sigma_{y},\tau_{x}\sigma_{z},\tau_{y}\sigma_{0},\tau_{z}\sigma_{0}\}$ and $\{\check{\gamma}_{1},\dots,\check{\gamma}_{7}\}=\{\mu_{x}\tau_{x}\sigma_{x},\mu_{x}\tau_{x}\sigma_{y},\mu_{x}\tau_{x}\sigma_{z},\mu_{x}\tau_{y}\sigma_{0},\mu_{x}\tau_{z}\sigma_{0},\mu_{y}\tau_{0}\sigma_{0},\mu_{z}\tau_{0}\sigma_{0}\}$.

Now, we discuss how the boundary mass terms can be introduced in lattice models. Following Geier et al. (2018); Trifunovic and Brouwer (2019), we briefly review the mechanism for generating the boundary mass terms from symmetric perturbations to the bulk Hamiltonian. An inversion-symmetric Dirac Hamiltonian without any perturbation is

[TABLE]

where $I$ denotes the operator for inversion symmetry. Inversion-symmetric perturbation $H_{1}$ is given by

[TABLE]

where $M_{a}=-IM_{a}I^{-1}$ with $a=1,\dots,d$.

After projection to the boundary, $H_{1}$ becomes

[TABLE]

where $\tilde{M}_{a}=M_{a}P_{+}(\bm{r})$. Then, one can identify the boundary mass terms as $m_{a,\bm{r}}=\sum_{i=1}^{d}\Delta_{i,a}\bm{n}_{\bm{r},i}$. Since $m_{a,\bm{r}}$ is proportional to the surface normal vector $\bm{n}_{\bm{r}}$, $m_{a,\bm{r}}=-m_{a,-\bm{r}}$ is naturally implemented. Finally, the lattice regularization with the substitution $k_{i}\rightarrow\sin k_{i}$ and $\lambda\rightarrow-(d-\lambda-\sum_{i=1}^{d}\cos k_{i})$ leads to the bulk Hamiltonian

[TABLE]

which naturally contains the boundary mass terms whose explicit form appears after projection to the boundary. All the tight-binding models used in the paper are constructed by using this method. For example, the tight-binding models for class A $2$D third-order and $3$D fourth-order TIs are given by

[TABLE]

and

[TABLE]

For the tight-binding models defined above, since the boundary mass terms are given by $m_{a,\bm{r}}=\sum_{i=1}^{d}\Delta_{i,a}\bm{n}_{\bm{r},i}$, the corresponding winding numbers are determined by the sign of the determinant of $\Delta_{i,a}$. For example,

[TABLE]

where we use $m_{a=1,2}=\Delta_{1,a+2}\cos\theta+\Delta_{2,a+2}\sin\theta$ and

[TABLE]

Similarly, we obtain $w_{d}=\text{sgn}(\text{Det}(\bm{\Delta}))$. Thus, the $d$D $(d+1)$th-order TI described by the relevant Dirac Hamiltonian is characterized by $w_{d}=1$ (mod $2$).

V.2 Numerical calculations

In order to verify our theoretical prediction for charge accumulation at the boundary based on the continuum Dirac Hamiltonian, here we calculate the energy spectrum and the eigenstates of $2$D third-order and $3$D fourth-order TIs described by the lattice Hamiltonian $H_{{\rm A}(2,3)}$ and $H_{{\rm A}(3,4)}$, by considering various shapes of open boundaries. We choose square and disk geometries for $2$D third-order TI, and cube and sphere geometries for $3$D fourth-order TI. For numerical calculations, we choose $\Delta_{1,3}=0.240$, $\Delta_{1,4}=-0.196$, $\Delta_{2,3}=-0.259$, $\Delta_{2,4}=-0.204$ for the $2$D third-order TI, and $\check{\Delta}_{1,4}=-0.25$, $\check{\Delta}_{1,5}=0.05$, $\check{\Delta}_{1,6}=0.10$, $\check{\Delta}_{2,4}=-0.05$, $\check{\Delta}_{2,5}=0.30$, $\check{\Delta}_{2,6}=-0.10$, $\check{\Delta}_{3,4}=0.05$, $\check{\Delta}_{3,5}=-0.05$, $\check{\Delta}_{3,6}=0.35$ for the $3$D fourth-order TI. Then $\text{sgn}({\text{Det}(\bm{\Delta}}))=-1$ and $\text{sgn}({\text{Det}(\bm{\check{\Delta}}}))=-1$ in both cases. Also, we assume that all sublattice sites are located at the unit-cell center for simplicity.

The energy spectra in Figs. 9(a) and 9(d) show that the number of occupied bands and that of unoccupied bands are different in the topological phase ($\lambda=1$) whereas they are the same in the trivial phase ($\lambda=-1$). Figures 9(a) and 9(d) correspond to the $2$D third-order TI and the $3$D fourth-order TI, respectively. Moreover, when the states below (above) the gap are fully occupied (unoccupied) by adding an extra electron or hole to the half-filled system, we find that the extra charge accumulation is localized only at certain subregions of the boundary where the surface and the bulk mass terms change rapidly. More explicitly, it is found that the extra charge is localized at the corners of the system with a square or cube geometry, while the extra charge is spread over the whole boundary region when the system geometry has a disk or sphere shape, as shown in Fig. 9.

The difference in the number of occupied and unoccupied states in $d$D $(d+1)$th-order TIs can be understood by examining the change in the energy spectrum of a $d$D $d$th-order TI with chiral symmetry when chiral-symmetry-breaking perturbations are added. Specifically, a $d$D $d$th-order TI with chiral symmetry has two zero-energy modes related by inversion. Adding chiral-symmetry-breaking perturbation, the zero-energy modes are pushed either upward or downward until they are merged to the unoccupied or occupied states, leading to a $d$D $(d+1)$th-order TI without in-gap states. As long as inversion symmetry is preserved, the two in-gap states should be degenerate during this process. The filling anomaly requires that an extra electron or hole is added to the half-filled system, and the states below (above) the gap are fully occupied (unoccupied), to preserve inversion symmetry. Then, one can examine the extra charge accumulation relative to the half-filling (see Sec. II.1 for the related discussion in $1$D second-order TIs).

However, to understand the distribution of the extra charge accumulation at the boundary, we need the additional information of the induced current $j^{\mu}_{d}$ as explained in Sec. III.3. To confirm the charge accumulation predicted by the induced current in Sec. III.3, we first numerically solve the lattice model for a $2$D third-order TI and compute the extra charge density $\Delta\rho(x,y)$ at the position $(x,y)$ relative to the half-filling ($\rho_{0}=2$). In both the square and disk geometries, the extra charge density $\Delta\rho(x,y)$ is defined as

[TABLE]

where $\psi_{i}(x,y)$ is the $i$th energy eigenstate and $N_{\rm occ}$ is the number of occupied states. Since total extra charge is $-1$ according to $\delta=\int j^{0}_{1}=\text{sgn}(\text{Det}(\bm{\Delta}))=-1$, the total sum of $\Delta\rho(x,y)$ is equal to $-1$. Namely,

[TABLE]

where $N_{x}$ and $N_{y}$ denote the number of lattice sites in the $x$ and $y$ directions, respectively.

In the case of the square geometry, the extra charge is mostly localized at the corners as shown in Fig. 9(b). This is simply because the boundary mass terms vary abruptly near the corners. To demonstrate the spatial distribution of the extra charge accumulation, we divide the system into four quadrants. Then we compare the total extra charge in each quadrant, $Q^{\rm num.}_{i=1,2,3,4}$ with the total induced charge density $j^{0}_{1}$ along the boundary of each quadrant, $Q_{i}$. Because of inversion symmetry, $Q^{\rm num.}_{1}=Q^{\rm num.}_{3}$ and $Q^{\rm num.}_{2}=Q^{\rm num.}_{4}$. Explicitly, $Q^{\rm num.}_{1,2}$ are defined as

[TABLE]

According to the numerical data displayed in Table 6, the deviation of $Q_{i}$ from $Q^{\rm num.}_{i}$ is only about one-percent. This is because when the bulk gap is sufficiently large, our formula for the induced current becomes more accurate due to the negligible contributions of higher-order terms in the derivative expansion. More detailed information about how the induced currents are computed is in Appendix G. On the other hand, in contrast to the case of the square geometry, the extra charge is spread over the whole boundary in the disk geometry as shown in Fig. 9(c). This delocalized charge accumulation is also consistent with the profile of $j^{0}_{1}(\theta)$ in Eq. (26).

Now, let us consider the $3$D fourth-order TI. Both in the cube and sphere geometries, the extra charge density $\Delta\rho(x,y,z)$ at the position $(x,y,z)$ compared to half-filling ($\rho_{0}=4$) is defined as

[TABLE]

where $\psi_{i}(x,y,z)$ is the $i$th energy eigenstate and $N_{\rm occ}$ is the number of occupied states. Also, the total accumulated charge is $1$ according to $\delta=\int j^{0}_{2}=-\text{sgn}(\text{Det}(\bm{\check{\Delta}}))=1$, and thus the total sum of $\Delta\rho(x,y,z)$ is equal to $1$. Namely,

[TABLE]

where $N_{i=x,y,z}$ denote the number of lattice sites in the $i$ direction.

In the case of the cube geometry, the extra charge is mostly accumulated at the corners where the boundary mass terms vary rapidly, as shown in Fig. 9(e). To demonstrate the spatial distribution of the extra charge accumulation, we divide the system into eight octants. Then, we compare the total extra charge in each octant, $Q^{\rm num.}_{i=1,\dots,8}$ with the total induced charge density $j^{0}_{2}$ along the boundary of each octant, $Q_{i}$. Inversion symmetry imposes $Q^{\rm num.}_{1}=Q^{\rm num.}_{7}$, $Q^{\rm num.}_{2}=Q^{\rm num.}_{8}$, $Q^{\rm num.}_{3}=Q^{\rm num.}_{5}$, and $Q^{\rm num.}_{4}=Q^{\rm num.}_{6}$. Explicitly, $Q^{\rm num.}_{1,2,3,4}$ are defined as

[TABLE]

The results of numerical calculations are displayed in Table 7, which shows that the deviation of $Q_{i}$ from $Q^{\rm num.}_{i}$ is about 10$\%$ percent. This is because our formula for the induced current in Eq. (III.3) get corrections from the higher-order terms when the bulk gap is small. When the bulk gap becomes larger, the deviation of $Q_{i}$ from $Q^{\rm num.}_{i}$ becomes smaller. In contrast to the case of the cube geometry, the extra charge is spread over the whole boundary in the sphere geometry as shown in Fig. 9(f). This delocalized charge accumulation is also consistent with the profile of $j^{0}_{2}(\theta,\phi)$ in Eq. (III.3). These numerical studies for $2$D third-order and $3$D fourth-order TIs in various open boundaries demonstrate the validity of the theory based on the induced current.

VI Wilson loop of fragile TI

As we have shown in Secs. IV and V, the inversion parities completely characterize inversion-protected fragile topological insulators. On the other hand, there is an alternative powerful method, called the Wilson loop method Yu et al. (2011); Alexandradinata et al. (2014, 2016); Fidkowski et al. (2011); Benalcazar et al. (2017a, b), that can be used to diagnosis the Wannier obstruction. Here, we demonstrate the consistency between the inversion-parity analysis and the Wilson loop method to complete the analysis of the fragile band topology.

The Wilson loop is an operator defined along a periodic direction in the Brillouin zone, whose spectral winding along the other periodic direction indicates the Wannier obstruction on the 2D subspace in the Brillouin zone Soluyanov and Vanderbilt (2011); Yu et al. (2011); Fidkowski et al. (2011). When the Wilson loop spectrum is gapped, the nested Wilson loop can be defined by treating the Wilson loop spectrum like the energy spectrum Benalcazar et al. (2017a, b). As shown below, the spectral winding of the Wilson loop and the nested Wilson loop spectra contain the essential characteristics of fragile TIs in three dimensions.

For completeness, in this section, we investigate the fragile topology of $d$D $(d+1)$th-order TIs by using the Wilson loop method and show that the Wilson loop analysis is consistent with the diagnosis based on the parity eigenvalues. First, we observe that the Wilson loop spectrum of fragile TIs has a non-zero relative winding Alexandradinata et al. (2014). Then, we show that the relative winding of the Wilson loop can be unwound by adding appropriately chosen trivial bands. In dimensions greater than two, however, the number of additional orbitals required to trivialize the fragile topology does not match with the number of additional orbitals required to unwind the relative winding of the Wilson loop spectrum. We find that their number difference counts the number of orbitals required to unwind the relative winding of the nested Wilson loop spectrum. Comparing the Wilson loop and nested Wilson loop spectra of $2$D third-order and $3$D fourth-order TIs, we show that the fragile topology of $3$D insulators is characterized by a non-zero relative winding of the nested Wilson loop, and its unwinding requires extra trivial bands to be added. We show that this idea is generally valid for $d$D $k$th-order TIs with $d\geq 3$ and $k>d$. In all dimensions $d$ and for all order $k$, the number of trivial bands needed to unwind all nested Wilson loops is consistent with the number of bands needed to trivialize the fragile topology determined by parity eigenvalues, which is indicated in Eq. (47).

VI.1 Fragile band topology and the relative winding number

The fragile band topology of a $d$D $(d+1)$th-order TI is manifested by the winding of the (nested) Wilson loop spectrum. The Wilson loop $W_{\bm{k}_{0}+\bm{G}\leftarrow\bm{k}_{0}}$ indicates the parallel transport of Bloch wave functions along a non-contractible loop in the Brillouin zone, from $\bm{k}_{0}$ to $\bm{k}_{0}+\bm{G}$ where $\bm{G}$ is a reciprocal lattice vector Yu et al. (2011); Alexandradinata et al. (2014, 2016). In a tight-binding model, the Bloch wave function is replaced by the eigenstate $|n,\bm{k}\rangle$ of the tight-binding Hamiltonian $H(\bm{k})$, which satisfies $H(\bm{k})|n,\bm{k}\rangle=E_{n}(\bm{k})|n,\bm{k}\rangle$. The Wilson loop is defined in terms of $|n,\bm{k}\rangle$ as

[TABLE]

where $\bm{k}_{i}=\bm{k}_{0}+\frac{i}{N}\bm{G}$ with an integer $i$, and $n_{\rm occ}$ denotes the number of occupied states. And the projector into the occupied states $P(\bm{k})$ is defined by $P(\bm{k})=\sum_{n=1}^{n_{\rm occ}}\lvert{n,\bm{k}}\rangle\langle{n,\bm{k}}\rvert$. Note that the Wilson loop spectrum $\{{\theta}(\bm{k_{0}})\}$, $W_{\bm{k}_{0}+\bm{G}\leftarrow\bm{k}_{0}}\lvert{\theta(\bm{k}_{0})}\rangle=e^{i\theta(\bm{k}_{0})}\lvert{\theta(\bm{k}_{0})}\rangle$, is independent of the start point $\bm{k}_{0}$. Instead of the Wilson loop spectrum, we sometimes use the equivalent term, the Wilson bands in the following discussions. Also, since we mostly consider the Wilson loop defined along the $k_{1}$ direction, we will denote $\{\theta_{1}(k_{2},k_{3},\dots)\}$ as the Wilson loop spectrum. The definition and properties of the Wilson loop is summarized in Appendix H.

Here, let us review the Wilson loop of inversion-symmetric systems studied comprehensively in Ref. Alexandradinata et al. (2014). For simplicity, let us consider a $2$D plane $(k_{1},k_{2})$ which is left invariant under the inversion. The Wilson loop spectrum satisfies $\{\theta_{1}(k_{2})\}=\{-\theta_{1}(-k_{2})\}$ and can exhibit chiral and/or relative winding. We always assume that there is no chiral winding since it represents the non-zero Chern number and we are interested in the insulators which are not stable TIs. Figure 10 shows the Wilson loop spectra for a $2$D third-order TI and a $3$D fourth-order TI. In Fig. 10(a) relevant to the $2$D third-order TI, one Wilson band has the winding number $1$ while the other one has the winding number $-1$. Their winding number difference is called the relative winding number $N_{W}$, and $N_{W}=1$ in this case. Similarly, Fig. 10(b) has $N_{W}=2$ while Figs. 10(c) and 10(d) have $N_{W}=0$.

Due to the constraint $\{\theta_{1}(k_{2})\}=\{-\theta_{1}(-k_{2})\}$, the Wilson loop spectrum at $\bar{k}_{2}=0,\pi$ consists of $\theta_{1}(\bar{k}_{2})=0$ or $\pi$ which are inversion symmetric by itself, or $(-\theta_{0},\theta_{0})$ which are inversion symmetric as a pair. Thus, we define $n^{W}_{0}(\bar{k}_{2})$, $n^{W}_{\pi}(\bar{k}_{2})$, and $n^{W}_{c}(\bar{k}_{2})$, which indicate the number of Wilson loop spectrum with eigenvalues $\theta(\bar{k}_{2})=0$, $\pi$, and $(-\theta_{0},\theta_{0})$, respectively. For example, the parity configuration $\mathcal{B}_{\gamma}(2,2)$ shown in Fig. 8(a) is mapped to $\{n^{W}_{\theta}(0)\}=(n^{W}_{0}(0),n^{W}_{\pi}(0),n^{W}_{c}(0))=(0,2,0)$ and $\{n^{W}_{\theta}(\pi)\}=(n^{W}_{0}(\pi),n^{W}_{\pi}(\pi),n^{W}_{c}(\pi))=(2,0,0)$. Then the relative winding number $N_{W}$ is given by Alexandradinata et al. (2014)

[TABLE]

which is unambiguously determined by $n_{\xi}(\bm{k}^{{\rm inv}})$ where $n_{\xi}(\bm{k}^{{\rm inv}})$ denotes the number of the occupied bands at $\bm{k}^{\rm inv}=(0,0),(\pi,0),(0,\pi)$, or $(\pi,\pi)$ with inversion parity $\xi=\pm 1$. The mapping between the inversion parities of occupied states $\{n_{\xi}(\bm{k}^{{\rm inv}})\}$ and the Wilson loop spectrum $\{n^{W}_{\theta}(\bar{k}_{2})\}$ is studied for arbitrary number of occupied states $n_{\rm occ}$ Alexandradinata et al. (2014), as summarized in Appendix. H.6.

Equation (VI.1) implies that the relative winding $N_{W}$ is fragile, i.e., $N_{W}$ decreases as appropriately chosen trivial bands are added. To see this, let us consider the lattice model for a $d$D $(d+1)$th-order TI in Eq. (52) with $\lambda=1$, which has $2^{d-1}$ occupied eigenstates whose parity configuration corresponds to $\mathcal{B}_{\gamma}(d,N_{\mathcal{B}}=2^{d-1})$. We focus on the $(k_{1},k_{2})$ plane with $k_{i>2}=0$. It is straightforward to show that the Wilson loop spectrum of $\mathcal{B}_{\gamma}(d,2^{d-1})$ is given by $\{n^{W}_{\theta}(0)\}=\{n^{W}_{0}(0),n^{W}_{\pi}(0),n^{W}_{c}(0)\}=\{0,2^{d-1},0\}$ and $\{n^{W}_{\theta}(\pi)\}=\{n^{W}_{0}(\pi),n^{W}_{\pi}(\pi),n^{W}_{c}(\pi)\}=\{2^{d-1},0,0\}$. Then, Eq. (VI.1) gives the relative winding number $N_{W}=2^{d-2}$. Now, let us examine how $N_{W}$ changes as trivial bands are added. In the $2$D unit cell in real space, which is reciprocal to the $(k_{1},k_{2})$ plane, there are four Wyckoff positions: $A=\bm{0}$, $B=\frac{1}{2}\bm{a}_{1}$, $C=\frac{1}{2}\bm{a}_{1}+\frac{1}{2}\bm{a}_{2}$, $D=\frac{1}{2}\bm{a}_{2}$ where $\bm{a}_{1,2}$ are relevant lattice vectors. Then, adding a $p(B)$, $p(C)$, or $p(D)$ orbital does not change $N_{W}$. Here, we use the expression $p(W)$ to denote the $p$ orbital at the Wyckoff position $W$. On the other hand, adding a $p(A)$ decreases $N_{W}$ by $1$. In general, the relative winding of the Wilson loop can be unwound by adding $2^{d-2}$ $p(A)$ orbitals for a $d$D $(d+1)$th-order TI. The above results are summarized in Table 8.

We verify numerically the fragile nature of the relative winding number discussed above. Figure 11 shows that the Wilson loop spectrum of $2$D third-order TIs changes as expected when $s(A)$ or $p(A)$ orbitals are added. Namely, adding a $p(A)$ orbital opens a gap in the Wilson loop spectrum while adding a $s(A)$ orbital cannot. The details of numerical calculations are summarized in Appendix. J. Since the Wilson bands are divided into two sectors ($P_{1,2}$) after we add a $p(A)$ orbital as shown in Fig. 11(b), the Berry phase $\tilde{\Phi}$ of each sector can be calculated. $P_{1,2}$ regroup the Wilson bands in an inversion-symmetric manner, and thus have quantized Berry phases $\tilde{\Phi}[P_{1,2}]$ in the $k_{2}$ direction. We find $\tilde{\Phi}[P_{1}]=\pi$ and $\tilde{\Phi}[P_{2}]=(0,\pi)$. On the other hand, let us note that the average $\Phi$ of the Wilson loop eigenvalues in $P_{1}$ and $P_{2}$ sectors along the $k_{1}$ direction are [math] in both cases. As the two Wilson bands in $P_{2}$ have a fixed symmetric value $\theta_{1}(k_{2})=\pi$ at $k_{2}=0$, we interpret $\Phi[P_{2}]=0$ in a way that the Wannier centers of the two Wilson bands in the $\bm{a}_{1}$ direction are $\frac{1}{2}$ in both cases. Thus, we conclude that the Wannier centers of three bands in $P_{1}$ and $P_{2}$ are $(0,\pi)$, $(\pi,0)$, and $(\pi,\pi)$, respectively. These values support that the Wannier state representation for the $2$D third-order TI is $\ominus p(A)\oplus p(B)\oplus p(C)\oplus p(D)$ Wieder and Bernevig (2018). This is also consistent with Eq. (47) and the corner charge accumulation Benalcazar et al. (2019) in the square boundary shown in Fig. 9(b).

Let us note that the fragility of $N_{W}=1$ Wilson bands has already been studied in Ref. Wang et al., 2018; Wieder and Bernevig, 2018. However, the fragility of $N_{W}=2$ Wilson bands in a $3$D fourth-order TI and in its higher-dimensional generalization has more involved structures and has not been studied yet. In the following, we show that, when $d>2$, adding $2^{d-2}$ $p$(A) orbitals are not sufficient to trivialize the fragile topology of $d$D $(d+1)$th-order TI even though it makes all the relative windings of the Wilson loops to be zero. This implies the nontrivial winding of the nested Wilson loops, which is addressed in the next section.

VI.2 Nested Wilson loop and trivialization of fragile topology

In this section, we discuss the fragility of the Wilson loop of the $3$D fourth-order TI. Symmetry data analysis shows that the trivialization of the $3$D fourth-order TI requires at least three $p(A)$ orbitals although the relative winding of the Wilson bands can be gapped with two $p(A)$ orbitals. To understand the origin of the mismatch in the number of required orbitals, we investigate the nested Wilson loop of the $3$D fourth-order TI with two or three additional $p(A)$ orbitals.

VI.2.1 Adding one $p(A)$ orbital

Adding one $p(A)$ orbital reduces the relative winding from $N_{W}=2$ to $N_{W}=1$ as shown in Figs. 12(a)-12(c). Thus, the winding of the Wilson loop persists, and thus the nontrivial band topology still exists.

VI.2.2 Adding two $p(A)$ orbitals

Adding two $p(A)$ orbitals reduces the relative winding from $N_{W}=2$ to $N_{W}=0$ as shown in Figs. 12(d)-12(f). Six Wilson bands are divided into two groups, $P_{\rm in}$ and $P_{\rm out}$, in an inversion-symmetric manner. Then, we can compute the nested Wilson loops of $P_{\rm in}$ and $P_{\rm out}$, which show non-zero nested Chern numbers $C_{\rm nest}=-2$ and $C_{\rm nest}=2$, respectively, as shown in Figs. 13(a) and 13(b). This can be explained by the fact that a $3$D fourth-order TI can be considered as two copies of the class A $3$D second order TI, i.e., the axion insulator, which is characterized by an odd nested Chern number Wieder and Bernevig (2018). However, let us note that an even value of the nested Chern number itself does not guarantee the nontrivial bulk topology. This is because only the parity of $C_{\rm nest}$ remains robust when a band crossing between Wilson bands in $P_{\rm in}$ and those in $P_{\rm out}$, or a Wilson band transition, happens while keeping the bulk gap Wieder and Bernevig (2018). Since two $p(A)$ orbitals cannot trivialize the band topology according to the symmetry data analysis, the nontrivial bulk topology must be manifested in the nested Wilson loop. We show that the nontrivial band topology appears in the form of the relative winding in the nested Wilson loop. After the band crossing between Wilson bands, the nested Wilson loop of $P_{\rm in}$ shows the nested relative winding number $N_{W,{\rm nest}}=1$ while the nested Wilson loop of $P_{\rm out}$ is completely unwound as shown in Figs. 13(c) and 13(d). In the nested Wilson loop spectra for $P_{\rm in}$ sector shown in Figs. 13(a) and 13(c), we have $|C_{\rm nest}|+2N_{W,{\rm nest}}=2$. This means that the fragile band topology manifests as the nontrivial winding of the nested Wilson loop spectrum encoded in either $C_{\rm nest}$ or $N_{W,{\rm nest}}$.

By generalizing the above results, we propose an indicator $N_{d,{\rm nest}}$ which characterizes the fragile band topology. First, let us define $N_{d,{\rm nest}}$ in a way similar to its counterpart $N_{d}$ in Eq. (VI.1) is defined for the Wilson loop:

[TABLE]

where $n_{\pi}^{W_{\rm nest}}(\bar{k}_{3})$ and $n_{c}^{W_{\rm nest}}(\bar{k}_{3})$ are the number of eigenvalues equal to $\pi$ and that equal to a generic value, respectively, in the nested Wilson loop spectrum computed at $\bar{k}_{3}=0,~{}\pi$. The indicator $N_{d,{\rm nest}}$ is invariant under a band crossing between the Wilson bands, and thus it determines the band topology of a given nested Wilson loop spectrum. We note that there may be one or more gaps in the Wilson loop spectrum, and all possible nested Wilson loop spectra have to be considered. Let us distinguish the indicator $N_{d,{\rm nest}}$ into the following ways: i) When $N_{d,{\rm nest}}$ is an odd integer, the nested Wilson loop spectrum shows an odd $C_{\rm nest}$. This implies the stable band topology. ii) When $N_{d,{\rm nest}}$ is a positive even integer, $|C_{\rm nest}|+2N_{W,{\rm nest}}=N_{d,{\rm nest}}$ holds for $|C_{\rm nest}|<N_{d,{\rm nest}}$. Thus, the nested Wilson loop shows the nested relative winding number equal to $N_{W,{\rm nest}}=\frac{1}{2}N_{d,{\rm nest}}$ even when the band crossing between the Wilson bands change $C_{\rm nest}$ to zero. iii) When $N_{d,{\rm nest}}$ is a negative even integer or zero, the nested Wilson loop can be unwound by a band crossing between the Wilson bands. The rules ii) and iii) can be applied to the nested Wilson loop spectrum of $P_{\rm in}$ and $P_{\rm out}$ which show $N_{d,{\rm nest}}=2$ and $-2$, respectively, as shown in Figs. 13(a)-13(d). More details of our claim and the relevant proof are provided in Appendix. I.

VI.2.3 Adding three $p(A)$ orbitals

Adding the third $p(A)$ orbital to the case considered above makes little changes in the Wilson loop spectra. For instance, the third $p(A)$ orbital does not change the nested Chern number $C_{\rm nest}$ as shown in Figs. 13(e) and 13(f). However, after the band crossing between Wilson bands in $P_{\rm in}$ and those in $P_{\rm out}$, the windings in the nested Wilson loops for both $P_{\rm in}$ and $P_{\rm out}$ have been unwound as shown in Figs. 13(g) and 13(h). This can be inferred from the fact that the nested Wilson loop spectrum of $P_{\rm in}$ and $P_{\rm out}$ in Figs. 13(e)-13(h) show $N_{d,{\rm nest}}=0$ and $-2$, respectively. Therefore, we conclude that the fragile TI is characterized by the nontrivial relative winding of the nested Wilson loop.

Since the Wannier obstruction is fully relaxed, we can find the Wannier state representations corresponding to the cases shown in Figs. 13(g) and 13(h). One interesting feature of the nested Wilson loops in Fig. 13(g) is that even the nest Wilson bands from $P_{\rm in}$ can be divided further into two different groups of nested Wilson bands, $P_{\rm in}(P_{\rm in})$ and $P_{\rm out}(P_{\rm in})$, which is possible since there is no winding or relative winding in the spectrum. Similarly, Fig. 13(h) shows that the nested Wilson bands from $P_{\rm out}$ can be divided into $P_{\rm in}(P_{\rm out})$ and $P_{\rm out}(P_{\rm out})$ (see Appendix H for more details). Then, the corresponding averaged eigenvalues $\Phi$ and the Berry phases $\tilde{\Phi}$ are given by $\Phi[P_{\text{in,out}}]=\Phi[P_{\text{in,out}}(P_{\text{in,out}})]=0$, $\tilde{\Phi}[P_{\rm in}(P_{\rm in})]=\pi$, $\tilde{\Phi}[P_{\rm out}(P_{\rm in})]=(0,\pi)$, $\tilde{\Phi}[P_{\rm in}(P_{\rm out})]=(0,\pi)$, and $\tilde{\Phi}[P_{\rm out}(P_{\rm out})]=(0,\pi)$. Here, $\Phi[P_{i}]$ ($\Phi[P_{i}(P_{j})]$) indicates the average eigenvalues of the Wilson band $P_{i}$ [the nested Wilson band $P_{i}(P_{j})$] while $\tilde{\Phi}[P_{i}]$ ($\tilde{\Phi}[P_{i}(P_{j})]$) denotes the Berry phase of the Wilson band $P_{i}$ [the nested Wilson band $P_{i}(P_{j})$]. The corresponding Wannier centers of the seven occupied states are $(0,\pi,0)$, $(\pi,0,0)$, $(\pi,\pi,0)$, $(0,0,\pi)$, $(0,\pi,\pi)$, $(\pi,0,\pi)$, and $(\pi,\pi,\pi)$, respectively. This supports that the Wannier state representation for the $3$D fourth-order TI is $\bigoplus_{i=2}^{8}\,p(W_{i})\ominus 3p(W_{1})$, which is consistent with Eq. (47). Also, the corner charge accumulation in the cubic geometry shown in Fig. 9(e) can be related to the Wannier state representation in the same way as in $2$D.

VI.3 Generalization to higher dimensions and higher-order

The idea presented in the previous section can be generalized further to $d$D $k$th-order TI with $d\geq 3$ and $k>d$. The key observation is that only the inner Wilson bands in $P_{in}$ of the trivialized $l$th nested Wilson loop have the relative winding number in their $(l+1)$th nested Wilson loop spectrum. We assume that the nested Chern numbers are always zero, which can always be achieved via the band crossing between Wilson bands.

Explicitly, according to the discussion in Sec. V, a $d$D $k$th-order TI has an equivalent Wannier state representation given by

[TABLE]

where $d+k=2n+1$ and $n=d,d+1,\dots$ for class A systems. Thus, the $d$D $k$th-order TI can be trivialized by adding $2^{n-d}(2^{d-1}-1)$ $p(W_{1})$ orbitals according to Eq. (67). Now, we give a heuristic argument suggesting that $2^{n-d}\times 2^{d-l-2}$ $p(W_{1})$ orbitals trivialize the relative winding of $l$th nested Wilson loop for $l=0,1,\dots,d-2$. Here, the 0th nested Wilson loop means the conventional non-nested Wilson loop. First, we consider the Wilson loop of $d$D $k$th-order TI defined along the $k_{d}$ direction. As shown in Sec. VI.1, it has the relative winding number $2^{n-d}\times 2^{d-2}$. Thus, adding $2^{n-d}\times 2^{d-2}$ $p(W_{1})$ orbitals unwinds the Wilson loop completely. Then the resulting Wannier state representation is

[TABLE]

Since there is no relative winding in the Wilson loop, the Wilson bands can be divided into two groups $P_{\rm in}$ and $P_{\rm out}$. Here, $P_{in}$ consist of the nearly flat Wilson bands near $\theta_{1}=0$ which come from the additional trivial orbitals. This indicates that the Wannier centers of the Wilson bands in $P_{in}$ in the $\bm{a}_{d}$ direction are all zero. Then, the Wannier state representation of the Wilson bands in $P_{in}$ is equivalent to

[TABLE]

where we projected out the Wyckoff positions $W_{i}$ ($i=2^{d-1}+1,\dots,2^{d}$) whose coordinate in the $\bm{a}_{d}$ direction is $1/2$. Since the Wannier state representation in Eq. (69), is equivalent to that for the $(d-1)$D $(k-1)$th-order TI, we can repeat a similar procedure to trivialize the first nested Wilson loop. Repeating this procedure recursively, we conclude that the inner Wannier bands $P_{in}$ of the $l$th nested Wilson loop are equivalent to the $(d-l-1)$D $(k-l-1)$th-order TI which is characterized by the relative winding number $2^{n-d}\times 2^{d-l-3}$. This completes our heuristic proof. Since adding $2^{n-d}\times 2^{d-l-2}$ $p(W_{1})$ orbitals can unwind the $l$th nested Wilson loop for $l=0,1,\dots,d-2$, the total number of $p(W_{1})$ orbitals required for complete trivialization is

[TABLE]

In this way, we can understand the total number of additional orbitals required to trivialize the fragile topology in terms of the relative winding numbers in all possible nested Wilson loops.

VII General pumping process

In this section, we describe the pumping process of inversion-symmetric TIs. The main claim is that the band topology of the ($d+1$)D $(k-1)$th-order TI can be understood in terms of the adiabatic pumping process between the $d$D $k$th-order TI and the trivial-parity insulator along the $(d+1)$th momentum direction. In particular, we discuss how the higher-order TIs with hinge states, i.e., $(d+1)$D $d$th-order TIs, can be constructed via the pumping process of $d$D $(d+1)$th-order TIs. This can be considered as the generalization of the idea proposed in Ref. Wieder and Bernevig, 2018.

VII.1 General pumping process

Let us consider an adiabatic pumping process of the parity configurations of $d$D $k$th-order TIs along the $k_{d+1}$ direction between $k_{d+1}=0$ and $k_{d+1}=\pi$ planes. We label the symmetry data in the $k_{d+1}=0$ plane as $\{n_{-}(\bm{k}^{{\rm inv}}_{1})-n_{+}(\bm{k}^{{\rm inv}}_{1}),\dots,n_{-}(\bm{k}^{{\rm inv}}_{2^{d}})-n_{+}(\bm{k}^{{\rm inv}}_{2^{d}})\}$ and those in the $k_{d+1}=\pi$ plane as $\{n_{-}(\bm{k}^{{\rm inv}}_{2^{d}+1})-n_{+}(\bm{k}^{{\rm inv}}_{2^{d}+1}),\dots,n_{-}(\bm{k}^{{\rm inv}}_{2^{d+1}})-n_{+}(\bm{k}^{{\rm inv}}_{2^{d+1}})\}$, respectively. At each $d$D subspace, the momentum-space symmetry data can be mapped to $2^{d}$ $\nu_{W}(k_{d+1})$ indices for $d$D inversion-symmetric insulator by using Eq. (39) as

[TABLE]

and

[TABLE]

Let us note that, thanks to the simple structure of inversion symmetry, the matrix $A_{d}$ in Eq. (39), which connects the symmetry data in real and momentum spaces, satisfies the following simple recursion relation:

[TABLE]

which is proved in Appendix. B. Then, $2^{d+1}$ $\nu_{W}$ indices for $(d+1)$D inversion-symmetric insulator are given by

[TABLE]

where $i=1,\dots,2^{d}$. Thus, the resulting phase can be diagnosed by using $\{\nu_{W}(0)\}$ and $\{\nu_{W}(\pi)\}$.

Equation (VII.1) shows that any fragile TI or atomic insulators in $(d+1)$D can be understood as a pumping process between fragile TI or atomic insulators. Since all $\nu_{W}$ indices for $(d+1)$D fragile TI or atomic insulator are integers, all $\nu_{W_{i}}(0)$ and $\nu_{W_{i}}(\pi)$ are also integers according to Eq. (VII.1). This means that the parity configurations at $k_{d+1}=0$ and $\pi$ correspond to either fragile TIs or atomic insulators.

In Sec. IV.3, it is mentioned that, when $n_{\rm occ}=2$, the fragile TI in higher dimensions can be obtained by a layer stacking of 2D fragile TI, and this process can also be described by a momentum-space pumping of fragile TIs. To illustrate this idea, let us consider a trivial pumping process in which the parity configurations at $k_{d+1}=0$ and $\pi$ are identical to each other, i.e., $\nu_{W_{i}}(0)=\nu_{W_{i}}(\pi)$ for $i=1,\dots,2^{d}$. In this situation, $\nu_{W}$ indices for $(d+1)$D inversion-symmetric insulator are given by $\nu_{W_{i}}=\nu_{W_{i}}(0)$ and $\nu_{W_{i+2^{d}}}=0$ for $i=1,\dots,2^{d}$. This shows that a fragile TI with $n_{\rm occ}=2$ in $3$D, $\{\nu_{W}\}=\{-1,1,1,1,0,0,0,0\}$, for example, can be understood as a trivial pumping process of $2$D fragile TI with $\{\nu_{W}(0)\}=\{\nu_{W}(\pi)\}=\{-1,1,1,1\}$ in the $k_{3}$ direction, consistent with the previous description.

VII.2 Hamiltonian mapping

Here we show that an adiabatic pumping between the $d$D $k$th-order TI and the trivial-parity insulator leads to the ($d+1$)D $(k-1)$th-order TI:

($d+1$)D $(k-1)$th-order TI

=

Pumping between a $d$D $k$th-order TI and the trivial-parity insulator

.

The Hamiltonians for a class A $d$D $k$th-order TI and a $(d+1)$D $(k-1)$th-order TI are given by

[TABLE]

and

[TABLE]

where the number of $\gamma$ matrices is $q=d+k=2n+1$ and their dimension is $2^{n}\times 2^{n}$. The bulk mass term $\mathcal{M}$ and the inversion symmetry operator $I$ are given by $\mathcal{M}=I=\gamma_{q}$. The parity configurations corresponding to $H_{{\rm A}(d,k)}$ and $H_{{\rm A}(d+1,k-1)}$ are $\mathcal{B}_{\gamma}(d,2^{n-1})$ and $\mathcal{B}_{\gamma}(d+1,2^{n-1})$, respectively. The $\nu_{W}$ indices for $\mathcal{B}_{\gamma}(d,2^{n-1})$ are given by

[TABLE]

while those for the trivial parity configuration of $2^{n-1}$ bands with only positive parities are

[TABLE]

Thus, the pumping between $\mathcal{B}_{\gamma}(d,2^{n-1})$ and the trivial parity configuration gives $2^{d+1}$ $\nu_{W}$ indices $2^{n-d-1}\{1-2^{d},1,1,\dots,1\}$, which is exactly the $\nu_{W}$ indices of $\mathcal{B}_{\gamma}(d+1,2^{n-1})$.

Also, one can explicitly construct a mapping between the Hamiltonians given by

[TABLE]

where $H_{{\rm A}(d,k)}(\bm{k},\lambda)$ is the Hamiltonian for a $d$D $k$th-order TI in Eq. (75). The map $\varphi_{\uparrow}[H]$ increases the dimension of the Hamiltonian $H$ by one, and describes the pumping process of the $d$D $k$th-order TI. Here $\lambda$ is a parameter governing the topological phase transition between the trivial and topological phases. $(d+k)$ should be an odd integer in order that the $d$D $k$th-order and $(d+1)$D $(k-1)$th-order TIs exist. $H_{{\rm A}(d,k)}(\bm{k},\lambda)$ is in the trivial (topological) phase when $\lambda<0$ ($0<\lambda<2$). For any $k_{d+1}$, the map $\varphi_{\uparrow}[H_{{\rm A}(d,k)}](\bm{k},k_{d+1},\lambda)=H_{{\rm A}(d+1,k-1)}$ preserves the bulk gap, since

[TABLE]

The nontrivial topology of $\varphi_{\uparrow}[H_{{\rm A}(d,k)}]$ originates from the band topology of $H_{{\rm A}(d,k)}$. $\varphi_{\uparrow}[H_{{\rm A}(d,k)}]$ is in the topological (trivial) phase of $H_{A(d,k)}$ at $k_{d+1}=0$ ($k_{d+1}=\pi$) when $0<\lambda<2$. Hence, the map $\varphi_{\uparrow}[H_{{\rm A}(d,k)}]$ describes the desired pumping process.

One can consider the following examples as an application of the pumping picture. For instance, a pumping of a $1$D second-order TI described by the SSH model without chiral symmetry leads to the $2$D first-order TI, the Chern insulator. A pumping of a $2$D third-order TI gives a $3$D second-order TI, the axion insulator Wieder and Bernevig (2018). Finally, a pumping of a $3$D second-order TI gives a $4$D first-order TI, the $4$D Chern insulator. In the case of other AZ classes, the Hamiltonian mapping is more restricted. For instance, in class AI, a pumping of a $2$D third-order TI produces $\mathbb{Z}_{2}$ monopole nodal line semimetal Morimoto and Furusaki (2014); Fang et al. (2015); Zhao et al. (2016); Zhao and Lu (2017); Bzdušek and Sigrist (2017); Song et al. (2018c); Ahn et al. (2018). In the case of class AII, a pumping of a $2$D third-order TI generates the doubled strong TI characterized by the $\mathbb{Z}_{4}$ symmetry indicator Po et al. (2017); Fang and Fu (2017); Khalaf et al. (2018). Finally, let us mention several related works discussing the pumping of a $2$D third-order TI. References van Miert and Ortix (2018); Trifunovic et al. (2019) discussed a charge pumping process of roto-inversion-symmetric systems. Reference Lee et al. (2019) discussed a fractional charge bound to a vortex, which can be understood as a charge-pumping process, in the systems with various crystalline symmetries including inversion symmetry. Although the relevant symmetries are different from ours, all these pumping processes are essentially the pumping of a 2D third-order TI, or 2D second-order TI when the chiral symmetry is present, to a 3D second-order TI.

VIII Conclusion

In this work, we systematically study the bulk and boundary properties of fragile TIs protected by inversion symmetry. In particular, we focus on the class A inversion-symmetric insulators which can be characterized by $(2^{d}+1)$ invariants, $n_{\rm occ}$ and $\{\nu_{W}\}$, defined in terms of symmetry data in momentum space. Using these invariants, stable TIs, fragile TIs, and atomic insulators are distinguished. We propose an efficient method for diagnosing fragile topological phases protected by inversion symmetry, and estimate the number of parity configurations corresponding to atomic insulators, fragile TIs, and stable TIs, respectively. Considering that the symmetry indicator cannot distinguish atomic insulators and fragile TIs, our diagnosis method goes beyond the recently proposed diagnosis schemes.

We have extended the notion of higher-order topology by generalizing the concept of $d$D $k$th-order TI to the cases with $k>d$, which is normally defined for $k\leq d$. We show that $d$D $k$th-order TIs with $k>d$ have fragile band topology. In contrast to the case of stable TIs where the relative topology between a trivial and a topological insulators can be described by a Dirac Hamiltonian unambiguously, the relative topology between a fragile TI and an atomic insulator is more subtle. In particular, the Dirac Hamiltonian approach to the fragile TI is meaningful only when we specify trivial insulator taking into account the parity of occupied states at all inversion-invariant momenta. In our case, we define the trivial-parity insulator having positive-parity occupied states at all inversion-invariant momenta as a reference trivial insulator. Then, the band inversion described by the Dirac Hamiltonian gives a fragile TI unambiguously. In fact, once we have the full symmetry data at all inversion-invariant momenta, we can completely analyze the fragile band topology protected by inversion symmetry which we have demonstrated rigorously.

Moreover, we derive the bulk-boundary correspondence of $d$D $(d+1)$th-order TI, and identify a minimal fragile phase with the filling anomaly as the $d$D $(d+1)$th-order TI. In other words, a fragile TI corresponding to the $d$D $(d+1)$th-order phase can have nontrivial boundary charge. We also explain the topological origin of the charge accumulation on the boundary. We show how the mass winding, the induced current, the filling anomaly, and corner charges are related to each other in a unified way. In the case of $d$D $k$th-order TIs with $k>(d+1)$ although there is no nontrivial charge accumulation, the net parity $\mathcal{I}$, serves as a bulk topological invariant for for any $k>d$. In fact, the net parity characterizes the bulk topology of all inversion-symmetric insulators without gapless boundary states. Also, from the careful study of Wilson loops and nested Wilson loops, we show that fragile topological phases are characterized by the relative winding in all possible Wilson loops. Finally, we discuss the implication of the pumping process of inversion-symmetric insulators.

We expect that all the concepts developed in this paper, such as the mass winding, topological induced currents, the net parity, and the diagnosis of fragile topology based on symmetry data, can be applied to other topological crystalline insulators in general. For example, the mass winding approach can be extended to fragile topological or atomic phases with rotation symmetries, considering the recent progress reported in Refs. Benalcazar et al., 2017b, a; Liu et al., 2019a; Benalcazar et al., 2019. We leave the extension of our theory to more general cases of fragile TIs protected by various space groups and magnetic space groups for future studies.

Note added—. Recently, we noticed that a paper by Song et al. Song et al. (2019) appeared on arXiv. They classified all fragile topological phases indicated by symmetry eigenvalues for the class AII and $230$ space groups. For inversion-symmetric cases, i.e., space group 2, linear inequalities for diagnosing fragile TIs are proposed in Ref. Song et al., 2019 which is equivalent to our nonlinear inequality in Fig. 5.

Acknowledgements.

We thank Sungjoon Park for his participation in an early stage of this project. Y.H. and J.A. were supported by IBS-R009-D1. 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-20190008), and the POSCO Science Fellowship of POSCO TJ Park Foundation (Grant No. 0426-20180002). This work was supported in part by the U.S. Army Research Office under Grant No. W911NF-18-1-0137.

Appendix A Derivation of diagnosis criterion

In this appendix, we derive our criterion $\sum_{W}|\nu_{W}|>n_{\rm occ}$ that distinguishes the fragile topological and atomic phases. In the following, we adopt two assumptions: (i) the band topology is not stable topological, i.e., it is either a fragile topological phase or an atomic phase, and (ii) we consider only the band topology protected by inversion symmetry, such that inversion parity determines the band topology exactly.

Let us first prove that the phases with $\sum_{W}|\nu_{W}|>n_{\rm occ}$ cannot be atomic insulator. Any atomic phase has $\mu_{W}$ integer greater than or equal to zero. Then, the following inequality holds:

[TABLE]

Thus, a phase which satisfies $\sum_{W}|\nu_{W}|>n_{\rm occ}$ cannot be atomic phase.

Now, let us prove that insulators satisfying $\sum_{W}|\nu_{W}|\leq n_{\rm occ}$ are atomic insulators. In other words, we can construct an equivalent Wannier state representation with $\mu_{W,\pm}\geq 0$. From the $n_{\rm occ}$ occupied states, we first construct $\sum_{W}|\nu_{W}|$ Wannier states,

[TABLE]

where $s_{W}=\text{sgn}(\nu_{W})$. By construction, we have $\mu_{W,\xi}=(\frac{s_{W}-\xi}{2})\nu_{W}=(|\nu_{W}|-\xi\nu_{W})/2\geq 0$. Then, we consider the remaining $n_{\rm occ}-\sum_{W}|\nu_{W}|$ states. Since the Wannier state representation in Eq. (82) satisfies $\mu_{W,-}-\mu_{W,+}=\nu_{W}$ already, adding the remaining states should not change $\nu_{W}$ at each $W$. It can be done by locating pairs of even- and odd-parity states at the Wyckoff positions. It is possible to do so because $n_{\rm occ}-\sum_{W}|\nu_{W}|=\sum_{W,\xi}(1+\xi s_{W})\mu_{W,\xi}$ and $(1\pm s_{W})$ are even integers. Thus, we can locate $\frac{1}{2}(n_{\rm occ}-\sum_{W}|\nu_{W}|)$ pairs of $s$ and $p$ orbitals at any $W$ without changing $\nu_{W}$.

In conclusion, an insulator satisfying $\sum_{W}|\nu_{W}|>n_{\rm occ}$ is fragile TI, and vice versa.

Appendix B Properties of $A_{d}$ and $\nu_{W}$

B.1 Representation of symmetric Wannier state

Let us consider a Wannier state $w_{\xi}(W)$ centered at the Wyckoff position $W$ in the unit cell labeled by the lattice vector $\bm{R}$. We denote such states as $\lvert{w_{\xi,\bm{R}}(W)}\rangle$. Inversion symmetry acts on $\lvert{w_{\xi,\bm{R}}(W)}\rangle$ as

[TABLE]

where $\xi=\pm 1$. Then, the corresponding Bloch wave function $\lvert{\psi_{\xi}^{W}(\bm{k})}\rangle$ is

[TABLE]

where $N_{\rm cell}$ is the number of unit cells. Its transformation property under inversion is inherited from $\lvert{w_{\xi,\bm{R}}(W)}\rangle$,

[TABLE]

This means that the Bloch wave function $\lvert{\psi^{W}_{\xi}(\bm{k})}\rangle$ constructed by $w_{\xi}(W_{i})$ has parity $\xi\eta_{ij}$ at inversion-invariant momentum $\bm{k}^{{\rm inv}}_{j}$, where we define $\eta_{ij}\equiv\xi e^{-2i\bm{k}^{{\rm inv}}_{j}\cdot W_{i}}$. In other words, $|w_{\xi,\bm{R}}(W_{i})\rangle$ generates $n_{\pm}(\bm{k}^{{\rm inv}}_{j})=\frac{1}{2}(1\pm\xi\eta_{ij})$ Bloch states with parity $\pm 1$ at $\bm{k}^{\rm inv}_{j}$.

B.2 Properties of $\eta_{ij}$

Since $\bm{k}^{{\rm inv}}_{i}=2\pi W_{i}$ and $\frac{2}{\pi}\bm{k}^{{\rm inv}}_{i}\cdot W_{j}\in\mathbb{Z}$, $\eta_{ij}=e^{-2i\bm{k}^{\rm inv}_{i}\cdot W_{j}}$ satisfies

[TABLE]

and

[TABLE]

Let first us prove $\sum_{i=1}^{2^{d}}\eta_{ij}=2^{d}\delta_{1j}$. When $j=1$, all $\eta_{i1}$ are equal to $1$ because $W_{1}=0$ by definition. Thus, $\sum_{i=1}^{2^{d}}\eta_{i1}=\sum_{i=1}^{2^{d}}1=2^{d}$. When $j\neq 1$, $\eta_{ij}$ always has a partner $\eta_{i^{\prime}j}=-\eta_{ij}$ for some $i^{\prime}$. Therefore, $\sum_{i=1}^{2^{d}}\eta_{ij}$ is zero. Combining the above relations, we obtain $\sum_{i=1}^{2^{d}}\eta_{ij}=2^{d}\delta_{j1}$.

In a similar manner, one can prove that $\sum_{l=1}^{2^{d}}\eta_{il}\eta_{lj}=\sum_{l=1}^{2^{d}}e^{-2i(\bm{k}^{{\rm inv}}_{i}-\bm{k}^{{\rm inv}}_{j})\cdot W_{l}}=2^{d}\delta_{ij}$. The proof is trivial when $i=j$. When $i\neq j$, either $\bm{k}^{\rm inv}_{i}$ or $\bm{k}^{\rm inv}_{j}$ is non-zero. Let us assume that $\bm{k}^{\rm inv}_{i}\neq 0$ without loss of generality. Then, $\eta_{il}\eta_{lj}=e^{-2i(\bm{k}^{\rm inv}_{i}-\bm{k}^{\rm inv}_{j})\cdot W_{l}}$ always has a partner $\eta_{il^{\prime}}\eta_{l^{\prime}j}=-\eta_{il}\eta_{lj}$ for some $l^{\prime}$. The summation of $\eta_{il}\eta_{lj}$ over $l$ thus vanishes.

B.3 Derivation of the mapping $A_{d}$

We derive the matrix $A_{d}$ that relates the number of the Wannier sates in real space to the number of occupied states of even and odd parity in momentum space. As we show in the first subsection, $|w_{\xi,\bm{R}}(W_{i})\rangle$ generates $n_{\pm}(\bm{k}^{{\rm inv}}_{j})=\frac{1}{2}(1\pm\xi\eta_{ij})$ Bloch states with parity $\pm 1$ at $\bm{k}^{\rm inv}_{j}$, a set of Wannier states $\sum_{W,\xi}\,\mu_{W,\xi}w_{\xi}(W)$ produces

[TABLE]

Thus, we obtain

[TABLE]

Equivalently, $\nu_{W_{i}}=\sum_{j=1}^{2^{d}}2^{-d}\eta_{ij}\left[n_{-}(\bm{k}^{{\rm inv}}_{i})-n_{+}(\bm{k}^{{\rm inv}}_{i})\right]$ because $\eta^{-1}=2^{-d}\eta$ in the matrix notation as shown in Eqs. (86) and (87). We introduce a new notation $A_{d}$ by

[TABLE]

where

[TABLE]

From Eqs. (86) and (87), one can find that $A_{d}$ satisfies

[TABLE]

For later use, we write Eq. (90) in another form,

[TABLE]

by using $n_{+}(\bm{k}^{{\rm inv}}_{i})+n_{-}(\bm{k}^{{\rm inv}}_{i})=n_{\rm occ}$.

Also, $A_{d+1}$ and $A_{d}$ have a recursion relation,

[TABLE]

This can be easily proven by using $\frac{1}{\pi}(\bm{k}^{{\rm inv}}_{i+2^{d}}-\bm{k}^{{\rm inv}}_{i})=2(W_{i+2^{d}}-W_{i})=(0,\dots,0,1)$ for $i=1,\dots,2^{d}$.

B.4 Properties of $\nu_{W}$

As $n_{\xi}(\bm{k}^{\rm inv})$ must be an integer which satisfies $0\leq n_{\xi}(\bm{k}^{{\rm inv}})\leq n_{\rm occ}$, $\{\nu_{W}\}$ cannot take any values. Most importantly, $\{\nu_{W}\}$ satisfy

[TABLE]

where $n$ is an integer less than or equal to $n_{\rm occ}$, according to Eq. (B.3). Also, $\nu_{W}$ indices have the following properties:

[TABLE]

Proof of $2^{d-1}\nu_{W_{i}}\in\mathbb{Z}$.

As $n_{-}(\bm{k}^{{\rm inv}})$ is an integer greater than or equal to zero and $A_{d}$ satisfies $2^{d}[A_{d}]_{ij}=\pm 1$, $\nu_{W_{i}}=\sum_{j=1}^{2^{d}}2[A_{d}]_{ij}n_{-}(\bm{k}^{{\rm inv}}_{j})$ (mod $1$) in Eq. (93) implies

[TABLE]

Proof of $2^{d-2}(\nu_{W_{i}}\pm\nu_{W_{j}})\in\mathbb{Z}$.

From $2^{d-1}\left([A_{d}]_{ij}\pm[A_{d}]_{i^{\prime}j^{\prime}}\right)=0$ or $\pm 1$, it follows that $\nu_{W_{i}}\pm\nu_{W_{i^{\prime}}}=\sum_{j=1}^{2^{d}}2\left([A_{d}]_{ij}\pm[A_{d}]_{i^{\prime}j}\right)n_{-}(\bm{k}^{{\rm inv}}_{j})$ for $i,i^{\prime}\neq 1$, which implies

[TABLE]

when $i,i^{\prime}\neq 1$. Similarly, $\nu_{W_{1}}\pm\nu_{W_{i}}=\sum_{j=1}^{2^{d}}2\left([A_{d}]_{1j}\pm[A_{d}]_{ij}\right)n_{-}(\bm{k}^{{\rm inv}}_{j})-n_{\rm occ}$ for $i\neq 1$ implies

[TABLE]

Thus, we obtain $2^{d-2}(\nu_{W_{i}}\pm\nu_{W_{j}})\in\mathbb{Z}$ for $i\neq j$. This completes the proof because when $i=i^{\prime}$ the property $b$ reduces to the property $a$, which we prove above.

Proof of $|\nu_{W_{i}}|\leq n_{\rm occ}$.

Let us first prove $|\nu_{W_{1}}|\leq n_{\rm occ}$. From Eqs. (92) and (93), we obtain

[TABLE]

Thus, $|\nu_{W_{1}}|\leq n_{\rm occ}$. We can prove $|\nu_{W_{i}}|\leq n_{\rm occ}$ for $i\neq 1$ in a similar way. For $i\neq 1$, the number of $[A_{d}]_{ij}=\pm\frac{1}{2^{d}}$ is equal. Let us relabel the $j$ indices such that $[A_{d}]_{ij}=\frac{1}{2^{d}}$ for $j=1,\dots,2^{d-1}$ and $-\frac{1}{2^{d}}$ for $j=2^{d-1}+1,\dots,2^{d}$. Then, we obtain

[TABLE]

Proof of $|\nu_{W_{i}}\pm\nu_{W_{j\neq i}}|\leq n_{\rm occ}$.

Let us first prove $|\nu_{W_{1}}\pm\nu_{W_{i\neq 1}}|\leq n_{\rm occ}$. We can relabel the $j$ indices for given $i$ such that $[A_{d}]_{ij}=\pm\frac{1}{2^{d}}$ for $j=1,\dots,2^{d-1}$ and $\mp\frac{1}{2^{d}}$ for $j=2^{d-1}+1,\dots,2^{d}$. Then, $2^{d-1}\left([A_{d}]_{1j}\pm[A_{d}]_{ij}\right)=1$ for $j=1,\dots,2^{d-1}$ and [math] for $j=2^{d-1}+1,\dots,2^{d}$. Thus,

[TABLE]

When $i,i^{\prime}\neq 1$ and $i\neq i^{\prime}$, we can relabel the $j$ indices such that 2^{d-1}\Big{(}[A_{d}]_{ij}\pm[A_{d}]_{i^{\prime}j}\Big{)}=1 for $j=1,\dots,2^{d-2}$, $-1$ for $j=2^{d-2}+1,\dots,2^{d-1}$, and [math] for $j=2^{d-1},\dots,2^{d}$. The relabeling of $j$ indices is achieved as follows. Since both $\bm{k}^{{\rm inv}}_{i}$ and $\bm{k}^{{\rm inv}}_{i^{\prime}}$ are inversion-invariant momenta, so is $\bm{k}^{{\rm inv}}_{i^{\prime}}-\bm{k}^{{\rm inv}}_{i}$. Thus, we can write it as $\bm{k}^{{\rm inv}}_{i^{\prime}}-\bm{k}^{{\rm inv}}_{i}=\frac{\Delta\bm{G}}{2}$ where $\Delta\bm{G}$ is a reciprocal lattice vector. Then, $\eta_{ij}\pm\eta_{i^{\prime}j}=\eta_{ij}(1\pm\exp(i\Delta\bm{G}\cdot W_{j}))$, thus $\eta_{ij}\pm\eta_{i^{\prime}j}=0$ for $j=2^{d-1}+1,\dots,2^{d}$ after the relabeling of $j$ indices. From $\sum_{j}\eta_{ij}\pm\eta_{i^{\prime}j}=0$, $\eta_{ij}\pm\eta_{i^{\prime}j}=\pm 1$ are equally distributed, say $\eta_{ij}\pm\eta_{i^{\prime}j}=1$ for $j=1,\dots,2^{d-2}$ and $-1$ for $j=2^{d-2}+1,\dots,2^{d-1}$. Finally,

[TABLE]

Proof of $|\sum_{W}\nu_{W}|\leq n_{\rm occ}$.

Using Eq. (92), we obtain

[TABLE]

Appendix C Relations between $\nu_{W}$ indices and symmetry indicator

Symmetry indicators are constructed from symmetry eigenvalues in momentum space, and these distinguish stable topological phases from fragile topological and atomic phases Po et al. (2017). Here, we relate $\nu_{W}$ indices with inversion symmetry indicators.

In $d$ dimensions, one can define the $\mathbb{Z}_{2^{d-1}}$ symmetry indicator $z_{2^{d-1}}$ for inversion symmetry,

[TABLE]

To see why it is defined modulo $2^{d-1}$, we connect it to $\nu_{W}$ indices. Let us recall that an atomic orbital localized at $W_{i}$ with inversion parity $\xi$ has parity $\xi\exp(2i\bm{k}^{{\rm inv}}_{j}\cdot W_{i})$ at inversion-invariant momentum $\bm{k}^{{\rm inv}}_{j}$. Therefore, we have

[TABLE]

where Eqs. (89) and (87) are used in the first and second lines. Since all fragile topological and atomic phases have $\nu_{W}$ indices as integers so that they have trivial $z_{2^{d-1}}$ indicator, $z_{2^{d-1}}$=0, this indicator is nontrivial only for stable topological phases.

While only $\nu_{W_{1}}$ determines the strong index $z_{2^{d-1}}$, other $\nu_{W}$ indices are related to the weak indices $z^{(w)}$. For example, in two dimensions, $\nu_{W_{i}}\pm\nu_{W_{j}}$ are constrained to be integers as shown in Appendix B. Thus, all stable topological phases detected by $\nu_{A}=\frac{1}{2}$ modulo $1$ have $\nu_{B,C,D}$ as half integers. Thus, $\nu_{B,C,D}$ do not form additional symmetry indicators. This follows from the fact that there is no weak topology in two dimensions because there is no Wannier obstruction coming from the $1$D band topology: all $1$D insulators are atomic insulators. On the other hand, we have three $\mathbb{Z}_{2}$ weak indices in three dimensions, $z^{(w)}_{i=1,2,3}$ Fu and Kane (2007). These generate seven distinct weak topological phases characterized by $(z^{(w)}_{1},z^{(w)}_{2},z^{(w)}_{3})=(0,0,1),(0,1,0),(1,0,0),(1,1,0),(1,0,1),(0,1,1),(1,1,1)$. The weak symmetry indicators,

[TABLE]

are expressed by $\nu_{W_{1}}$, $\nu_{W_{2}}$, $\nu_{W_{3}}$, and $\nu_{W_{5}}$ as follows:

[TABLE]

Note that symmetry-indicated stable topological phases can be either insulators or semimetals Po et al. (2017). For example, in three dimensions, the phase is gapless when the strong indicator $z_{4}$ is odd Turner et al. (2012); Hughes et al. (2011).

Appendix D Filling anomaly from inversion of minimal number of bands corresponding to $d$D $(d+1)$th-order Dirac Hamiltonian

In this appendix, we discuss how an odd number of band inversions corresponding to $d$D $(d+1)$th-order Dirac Hamiltonian generates the filling anomaly. Here, we mean a band inversion at inversion-invariant momentum $\bm{k}^{\rm inv}_{i}$ by a change of parity configuration by $\Delta n_{-}(\bm{k}^{\rm inv}_{i})=\pm 2^{d-1}$. A change of $\nu_{W}$ indices due to a band inversion at $\bm{k}^{\rm inv}_{i}$, $\bm{\Pi}(\bm{k}^{\rm inv}_{i})$, is readily obtained from Eq. (39):

[TABLE]

Here, $\pm$ corresponds to $\Delta n_{-}(\bm{k}^{\rm inv}_{i})=\pm 2^{d-1}$. For example, band inversions at $\bm{k}^{\rm inv}_{1}$ and $\bm{k}^{\rm inv}_{2}$ in $3$D correspond to $\bm{\Pi}(\bm{k}^{\rm inv}_{1})=\pm\{1,1,1,1,1,1,1,1\}$ and $\bm{\Pi}(\bm{k}^{\rm inv}_{2})=\pm\{1,-1,1,-1,1,-1,1,-1\}$, respectively. Since $\Pi_{i}(\bm{k}^{\rm inv}_{j})=\pm 1$, changes in $\nu_{W_{i}}$ for $i=1,\dots,2^{d}$ are odd (even) integers simultaneously when the band inversions occur odd (even) number of times.

First, let us suppose that we have an insulator without the filling anomaly, i.e., $\nu_{W_{i}}$ are even integers for $i=2,\dots,2^{d}$. After an odd (even) number of band inversions, $\nu_{W_{i}}$ become odd (even) integers for $i=2,\dots,2^{d}$ and this implies the presence (absence) of the filling anomaly. Likewise, an odd (even) number of band inversions implies the absence (presence) of the filling anomaly when we start from an insulator with the filling anomaly, i.e., $\nu_{W_{i}}$ are odd integers for $i=2,\dots,2^{d}$.

Conversely, any inversion-symmetric insulators with (without) the filling anomaly become atomic insulator without (with) the filling anomaly after an odd (even) number of band inversions. For example, when $n_{\rm occ}=6$, $\{\nu_{W}\}=\{-1,1,3,1,1,-1,1,1\}$ has the following parity configuration.

[TABLE]

We always re-order the energy levels so that odd-parity states are atop even-parity states. After the band inversions at $\bm{k}^{\rm inv}_{1,2,5}$, $\nu_{W}$ indices changes by $-\sum_{i=1,2,5}\bm{\Pi}(\bm{k}^{\rm inv}_{i})$, thus $\{\nu_{W}\}=\{-4,0,0,0,0,0,0,2\}$:

[TABLE]

Now the parity configuration is equivalent to $4s(W_{1})\oplus 2p(W_{8})$, which is an atomic insulator without the filling anomaly. Hence, the number of band inversions modulo $2$, which is equal to $N_{\rm BI}=\sum_{i=1}^{2^{d}}\left\lfloor\frac{n_{-}(\bm{k}^{\rm inv}_{i})}{2^{d-1}}\right\rfloor$ (mod $2$), determines the filling anomaly, i.e., $N_{\rm BI}=0$ ($N_{\rm BI}=1$) for the absence (presence) of the filling anomaly. Here, $\left\lfloor x\right\rfloor$ denotes the greatest integer less than or equal to $x$.

Appendix E Layer construction of inversion-symmetric TIs and TSCs

Inversion-symmetric $k$th-order TIs/TSCs for $k\leq d$ in $d$ dimensions can be constructed by staking first-order TIs/TSCs in $(d-k+1)$ dimensions Khalaf (2018). Based on the layer construction, higher-order topological phases protected by inversion symmetry are neatly summarized in Table 9, which is reproduced from Ref. Khalaf, 2018. On the other hand, it was pointed out in Refs. Trifunovic and Brouwer, 2019; Song et al., 2017b; Huang et al., 2017 that topological phases without gapless boundary states are given by the layer construction starting from [math]D states. In this way, any $d$D systems are constructed by stacking [math]D states (see Fig. 14). Here, we extended the layer construction in Ref. Khalaf, 2018 to the case for $k>d$, where no gapless boundary states appear. We focus on the case where inversion-symmetry operator $I$ satisfies the commutation relations with $T$ and $P$ in the same way as the Hamiltonian, i.e., $[I,T]=0$ and $\{I,P\}=0$. When $I$ satisfies commutation relations different from that of $H$, only order of $k\leq 2$ or 3 can be realized depending on the AZ classes, and we refer readers to Khalaf (2018); Trifunovic and Brouwer (2019) for discussions about such cases.

The layer construction can be employed to construct the Hamiltonian for generalized $k$th-order TIs/TSCs with $k>d$ as was done for $k\leq d$ in Ref. Khalaf, 2018. Before introducing a dimensional raising mapping of Hamiltonian based on the layer construction, we adopt the following convention:

$\mathcal{H}_{(d,k)}$ and $I_{(d,k)}$ denote the Dirac Hamiltonian and inversion symmetry operator of $d$D $k$th-order TI/TSC. 2. 2.

$T_{(d,k)}$, $P_{(d,k)}$, and $S_{(d,k)}$ denote time-reversal, particle-hole, and chiral symmetry operators for $d$D $k$th-order TI/TSC. These non-spatial symmetries exist or not according to the AZ class. 3. 3.

$\mathcal{M}_{(d,k)}$ denotes the bulk mass term and satisfies

[TABLE]

Given a Hamiltonian $\mathcal{H}_{(d,k)}$, $\mathcal{H}_{(d+1,k+1)}$ is obtained by the dimensional raising mapping

[TABLE]

where $m_{k,\bm{r}}=-m_{k,-\bm{r}}$ denotes the boundary mass term. Note that the boundary mass terms, $m_{a,\bm{r}}$, for $a=1,\dots,k-1$, in $\mathcal{H}_{(d,k)}$ are deformed to be odd under inversion symmetry in $(d+1)$D once the dimensional raising mapping Eq. (126) is done. The symmetry operators of $(d+1)$D $(k+1)$th-order TI/TSC are given by

[TABLE]

thus, the dimensional raising mapping (126) does not change the AZ class. Note that Eq. (126) is understood as a continuum limit of

[TABLE]

where $\mathcal{M}_{(d+1,k+1)}=\mathcal{M}_{(d,k)}\otimes\tau_{z}$. Accordingly, the Hamiltonian for $d$D $(d+n)$th-order TI/TSC, $\mathcal{H}_{(d,d+n)}$ for $n>0$ and $d\geq 1$, is constructed by a “[math]D $n$th-order” Hamiltonian, $\mathcal{H}_{(0,n)}$ by using the dimensional raising mapping repeatedly:

[TABLE]

We define a [math]D $n$th-order Hamiltonian $\mathcal{H}_{(0,n)}$,

[TABLE]

where $\bm{R}$ is regarded as a set of auxiliary parameters that are promoted to $\bm{r}$ after the construction of $\mathcal{H}_{(d,d+n)}$ is done. Thus, $m_{a,\bm{R}}$ transforms under the symmetry operators,

[TABLE]

Now, let us regard $m_{a,\bm{R}}$ as momentum-like variables. Then, $T$ and $P$ can be regarded as space-time inversion symmetry and particle-hole symmetry combined with inversion symmetry. Hence, the [math]D Hamiltonian $\mathcal{H}_{(0,n)}$ is identical to the inversion-symmetric $(n-1)$D Dirac Hamiltonian in the presence/absence of $IT$, $IP$, and $S$ Bzdušek and Sigrist (2017). Therefore, $n$ is allowed when $n+s=1$ (mod $2$) in complex AZ classes, and $n+s=1,2,3,5$ (mod $8$) in real AZ classes where $s=0,1$ (mod $2$) and $s=0,\dots,7$ (mod $8$) label the complex and real AZ classes, respectively. This is summarized in Table 9, and it is identical to the classification table for $k\geq d$ in Ref. Khalaf, 2018. Therefore, our analysis generalizes Table 9 to arbitrary non-negative integers $k$ and $d$. Here, the inversion-protected $d$D $k$th-order topology has a $\mathbb{Z}_{2}$ classification, because a two-copy of $d$D $k$th-order TIs is not a $d$D $k$th-order TI anymore. Instead, it has a higher order. In our Dirac Hamiltonian construction in class A, a two-copy of $d$D $k$th-order TIs gives $d$D $(k+2)$th-order TI. For example in class AII, a two-copy of $d$D $(d+1)$th-order TIs gives $d$D $(d+5)$th-order TI.

Appendix F Boundary projection

In this appendix, we review the boundary projection introduced in Refs. Khalaf, 2018; Khalaf et al., 2018; Geier et al., 2018; Trifunovic and Brouwer, 2019. Consider the Dirac Hamiltonian in $d$ dimensions, $\mathcal{H}(\bm{r})$:

[TABLE]

Here, $\gamma_{i=1,\dots,d}$ are gamma matrices and $\mathcal{M}$ is bulk mass term. gamma matrices $\gamma_{i}$ and the bulk mass term $\mathcal{M}$ satisfy $\{\gamma_{i},\gamma_{j}\}=2\delta_{ij}$ and $\{\gamma_{i},\mathcal{M}\}=0$. We consider a system that is in the topological phase. Let us take $\lambda>0$ in the topological phase and $\lambda<0$ in the trivial phase. Then, $\lambda$ is positive inside the system while $\lambda$ is negative outside the system. We set $\lambda(\bm{r})$ such that it approaches to $+1$ as $\bm{r}$ approaches the center of the system, $\bm{r}=0$, and it approaches to $-1$ as $\bm{r}$ goes to infinity. On the boundary, $\lambda=0$ because $\lambda$ changes sign there.

For our analysis of the boundary theory of the Dirac Hamiltonian, we define $(d-1)$ orthonormal tangent basis vector $\bm{e}_{\bm{r},L=1,\dots,d-1}$ and the normal vectors $\bm{n}_{\bm{r}}$. On the tangent space at $\bm{r}$, the Dirac Hamiltonian in Eq. (132) is expressed as

[TABLE]

Here, $\bm{k}_{S}$ denotes the boundary momentum normal to $\bm{n}_{\bm{r}}$, $\bm{k}_{S}=\bm{k}-(\bm{k}\cdot\bm{n}_{\bm{r}})\bm{n}_{\bm{r}}$. Also, we use the vector notation for gamma matrices, $\bm{\gamma}=(\gamma_{1},\dots,\gamma_{d})$, and $x_{\perp}$ denotes the normal coordinate at $\bm{r}$ on the boundary, $\delta x_{\perp}=\bm{n}_{\bm{r}}\cdot\delta\bm{r}$. Note that the kinetic term, the first term in Eq. (133), can be expressed in terms of $\bm{e}_{\bm{r},L}$:

[TABLE]

where $(\bm{k}_{S})_{L}=\bm{k}_{S}\cdot\bm{e}_{\bm{r},L}$.

When the system is sufficiently large, $\mathcal{H}_{\perp}(\bm{k}_{S},\bm{r})$ has eigenstate $\psi(\bm{k}_{S},\bm{r})$ that is exponentially localized on the boundary. One can show that the Schrödinger equation,

[TABLE]

is solved by the following ansatz:

[TABLE]

where we define projection operators $P_{\pm}(\bm{r})=\frac{1}{2}(\mathds{1}\mp i(\bm{n_{r}}\cdot\bm{\gamma})\mathcal{M})$ is projection operator. The choice of projection operator, $P_{+}(\bm{r})$ or $P_{-}(\bm{r})$, depends on the profile of $\lambda(\bm{r})$. In our case, since $\lambda=-1$ outside the system, $P_{+}(\bm{r})$ is the proper choice. $P_{\pm}(\bm{r})$ satisfies the following relation:

[TABLE]

With the help of the ansatz in Eqs. (136) and (137), the eigenvalue problem in Eq. (135) is reduced to the form,

[TABLE]

where we define the boundary Hamiltonian $h(\bm{k_{S}},\bm{r})\equiv P_{+}(\bm{r})\mathcal{H}_{\perp}(\bm{k}_{S},\bm{r})P_{+}(\bm{r})$.

Now, we consider the $k$th-order topology, which is defined by the presence of $(k-1)$ odd-parity mass terms $m_{a,\bm{r}}M_{a}$ with $a=1,\dots,k-1$, where $M_{a}$ satisfy $\{\gamma_{i},M_{a}\}=0$, $\{\mathcal{M},M_{a}\}=0$, and $\{M_{a},M_{b}\}=2\delta_{ab}$. Inversion symmetry under $I=\mathcal{M}$ requires that $m_{a,\bm{r}}$ is odd under inversion, i.e., $m_{a,\bm{r}}=-m_{a,-\bm{r}}$. Those mass terms break translation symmetry because of their $\bm{r}$ dependence, and we assume that they are finite only near the boundary.

After the boundary mass terms are included, the boundary Hamiltonian becomes

[TABLE]

The eigenvalue problem for $h(\bm{k_{S}},\bm{r})$ can be solved in the following way. The boundary-projected gamma matrices and mass matrices $\tilde{\gamma}_{L}$ and $\tilde{M}_{a}$, are defined by the projection,

[TABLE]

and satisfy $\{\tilde{\gamma}_{L},\tilde{\gamma}_{M}\}=2\delta_{LM}P_{+}(\bm{r})$, $\{\tilde{M}_{a},\tilde{M}_{b}\}=2\delta_{ab}P_{+}(\bm{r})$ and $\{\tilde{\gamma}_{L},\tilde{M}_{a}\}=0$. Hence, energy eigenvalues are given by $E(\bm{k}_{S},\bm{r})=\pm\sqrt{\bm{k}_{S}^{2}+\sum_{a=1}^{k-1}m_{a,\bm{r}}^{2}}$.

So far, the boundary-projected mass terms $m_{a,\bm{r}}\tilde{M}_{a}$ are generated by the position-dependent mass terms $m_{a,\bm{r}}M_{a}$. Alternatively, $m_{a,\bm{r}}\tilde{M}_{a}$ can be generated by a constant perturbation $\mathcal{H}_{\rm pert}$ Geier et al. (2018); Trifunovic and Brouwer (2019),

[TABLE]

After the boundary projection, $\mathcal{H}_{\rm pert}$ is reduced to

[TABLE]

We use this method in numerical calculations.

Appendix G Induced current

In this appendix, we derive the induced current $J^{\mu}_{d}$ in Eq. (III.3). Let us consider the Dirac Hamiltonian coupled to $(d+1)$ slowly varying background scalar fields $m_{a,\bm{r}}$,

[TABLE]

where $\gamma$ denotes gamma matrix and $M_{a}=\gamma_{a+d}$. gamma matrices satisfy $\{\gamma_{i},\gamma_{j}\}=2\delta_{ij}$ and are chosen to be $\gamma_{1}\gamma_{2}\dots\gamma_{2d+1}=\chi_{d}\mathds{1}_{2^{d}}$ where $\chi_{d}=i^{d}$. Thus, $\text{Tr}[\gamma_{i_{1}}\gamma_{i_{2}}\dots\gamma_{i_{2d+1}}]=\epsilon_{i_{1}i_{2}\dots i_{2d+1}}\text{Tr}[\chi_{d}\mathds{1}_{2^{d}}]=\epsilon_{i_{1}i_{2}\dots i_{2d+1}}\chi_{d}2^{d}$ where $\epsilon_{12\dots 2d+1}=1$. The Lagrangian is given by $\mathcal{L}=i\Psi^{\dagger}\partial_{t}\Psi-\bar{\Psi}\mathcal{H}\Psi$ where $\Psi$ denotes the Dirac spinor and $\bar{\Psi}=\Psi^{\dagger}\gamma^{0}$:

[TABLE]

where we use the Feynman slash notation, $\not{A}=\sum_{\mu=0}^{d}A_{\mu}\gamma^{\mu}$, and define new gamma matrices $\gamma^{0}=\gamma_{2d+1}$ and $\gamma^{i\neq 0}=\gamma^{0}\gamma_{i}$. $\gamma^{\mu}$ satisfies $\{\gamma^{\mu},\gamma^{\nu}\}=-2\eta^{\mu\nu}$ where $\eta_{\mu\nu}=\text{diag}(-1,1,\dots,1)$ for $\mu,\nu=0,1,\dots,2d$, and $\text{Tr}[\gamma^{\mu_{1}}\gamma^{\mu_{2}}\dots\gamma^{\mu_{2d+1}}]=\epsilon^{\mu_{1}\mu_{2}\dots\mu_{2d+1}}\chi_{d}(-1)^{d}2^{d}$, where $\epsilon^{012\dots 2d}=1$.

We calculate the current density $\langle{j^{\mu}_{d}(x)}\rangle$ by using

[TABLE]

where $G(x,y)=i\langle{T\Psi(x)\bar{\Psi}(y)}\rangle$ is the propagator. Here, we denote $x^{\mu}=(x^{0},\bm{x})$ as $x$ for simplicity. The propagator $G(x,y)$ satisfies

[TABLE]

Regarding $G(x,y)$ as a function of $x$ and $(x-y)$, we Fourier transform the propagator by Midorikawa (1985)

[TABLE]

Note that translation symmetry is broken by $m_{a,\bm{r}}$, so $G(x,y)$ can not be a function of $(x-y)$ only. Then, we expand $\tilde{G}(x,p)$ by its spatial derivative terms:

[TABLE]

Here, we define $\tilde{G}_{0}(x,p)$:

[TABLE]

where $m_{\bm{x}}^{2}=\sum_{a=1}^{d+1}m_{a,\bm{x}}^{2}$ and $M_{\bm{x}}=m_{d+1,\bm{x}}-\sum_{a=1}^{d}m_{a,\bm{x}}\gamma^{a+d}$. The lowest-order contribution comes from $k=d$ term in Eq. (G). Up to the lowest-order in the derivative, we have

[TABLE]

From

[TABLE]

we obtain

[TABLE]

where $\hat{m}_{a,\bm{x}}=\frac{m_{a,\bm{x}}}{m_{\bm{x}}}$ and ${\rm Area}(S^{d})=\frac{2\pi^{(d+1)/2}}{\gamma(\frac{d+1}{2})}$. Hence, the induced current $\langle{j^{\mu}_{d}(x)}\rangle$ in $d$ dimensions (denoted as $J^{\mu}_{d}$ in the main text) is given by the winding number density of mass terms $m_{a,\bm{r}}$ where $a=1,\dots,d$.

The induced current density can be calculated either using the bulk Hamiltonian or the boundary Hamiltonian. There is a caveat because the convention for the gamma matrices may change by the projection, which can lead a different choice of $\chi_{d}$ on the boundary while we take $\chi_{d}=i^{d}$ in the bulk. Therefore, we need to carefully calculate $\tilde{\chi}_{d-1}\equiv 2^{-(d-1)}\text{Tr}[\tilde{\gamma}_{1}\dots\tilde{\gamma}_{d-1}\tilde{M}_{1}\dots\tilde{M}_{d}]$ on the boundary, where $\tilde{\gamma}$ and $\tilde{M}$ are gamma matrices and mass matrices projected to boundary defined in Appendix. F. It turns out that $\tilde{\chi}_{d-1}=-\chi_{d-1}$ as we show now.

Now, let us turn to the boundary theory when a domain between trivial and topological phases is formed. The boundary Hamiltonian $\bm{h}(\bm{k}_{S},\bm{r})$ is

[TABLE]

which is obtained by the projection to the boundary, introduced in Appendix F. Since $h(\bm{k}_{S},\bm{r})$ has $(d-1)$ kinetic terms and $d$ mass terms, the induced current $\langle{j^{\mu}_{d-1}(x)}\rangle$ on the $(d-1)$D boundary is given by repeating the above calculation. The minus sign can be traced in the following way. From the definition of $\tilde{\gamma}_{L=1,\dots,d-1}=(\bm{e}_{\bm{r},L}\cdot\bm{\gamma})P_{+}(\bm{r})$ and $\tilde{M}_{a=1,\dots,d}=M_{a}P_{+}(\bm{r})$, where $\bm{e}_{\bm{r},L}$ are orthonormal bases at $\bm{r}$, we show $\text{Tr}[\tilde{\gamma}_{1}\dots\tilde{\gamma}_{d-1}\tilde{M}_{1}\dots\tilde{M}_{d}]=-2^{d-1}\chi_{d-1}$ explicitly:

[TABLE]

where we define $\bm{e}_{\bm{r},0}\equiv\bm{n}_{\bm{r}}$ in the third line, $\text{Det}(\bm{e}_{\bm{r}})$ is the determinant of $\bm{e}_{\bm{r}}$ in the matrix notation.

Finally, let us mention that $d$D $(d+1)$th-order phase is allowed for classes A, AI, AII, BDI, and D. However, when the chiral or particle-hole symmetry is exists, the induced current is zero. Thus, charge accumulation on boundary occurs only for classes A, AI, and AII.

Appendix H Review of the Wilson loop

In this appendix, we review the Wilson loop and its symmetry transformation. We closely follow Ref. Alexandradinata et al., 2016. In addition, we introduce a more general method which exploits parallel transport and symmetry transformation of Berry connection.

H.1 Tight-binding notation

In the tight-binding model, Hamiltonian is spanned by a basis set $\lvert{\bm{R}\,\alpha}\rangle$. We denote $\bm{R}$ and $\bm{x}_{\alpha}$ as the unit cell position and the sublattice position of the orbital $\alpha$:

[TABLE]

Here, $\lvert{\bm{k}\,\alpha}\rangle$ is the Fourier transformation of $\lvert{\bm{R}\,\alpha}\rangle$,

[TABLE]

where $N_{\rm cell}$ is the number of unit cells in the system. These basis states are orthonormal: $\langle{\bm{R}\,\alpha|\bm{R}^{\prime}\,\beta}\rangle=\delta_{\bm{R},\bm{R}^{\prime}}\delta_{\alpha\beta}$ and $\langle{\bm{k}\,\alpha|\bm{k}^{\prime}\,\beta}\rangle=\delta_{\bm{k},\bm{k}^{\prime}}\delta_{\alpha\beta}$.

When the orbitals are located at $\bm{x}_{\alpha}\neq 0$, $\lvert{\bm{k}\,\alpha}\rangle$ is not periodic in the Brillouin zone.

[TABLE]

where $V(\bm{k})_{\alpha\beta}=e^{-i\bm{k}\cdot\bm{x}_{\alpha}}\delta_{\alpha\beta}$ (note that we take a convention different from the one in Ref. Alexandradinata et al., 2016). Diagonalizing $H(\bm{k})_{\alpha\beta}$, we obtain energy eigenstates $\lvert{n,\bm{k}}\rangle_{\alpha}$, where

[TABLE]

The energy eigenstates $\lvert{n,\bm{k}}\rangle$ satisfy

[TABLE]

where $n_{\rm tot}$ is the rank of Hamiltonian $H(\bm{k})_{\alpha\beta}$, which is also the number of orbitals: $\sum_{\alpha}1=n_{\rm tot}$. Since $H(\bm{k}+\bm{G})=V(\bm{G})\,H(\bm{k})\,V(\bm{G})^{-1}$, the “periodic gauge” is defined by

[TABLE]

H.2 Symmetry transformation

Consider a symmetry operation $\hat{\sigma}$: $\bm{r}\rightarrow D_{\sigma}\bm{r}+\bm{d}$. $D_{\sigma}$ is representation of the point-group part of $\hat{\sigma}$, which is an element of orthogonal group, $O(d)$. It acts on the basis $\lvert{\bm{R}\,\alpha}\rangle$ as follows:

[TABLE]

where $\tilde{\bm{R}}=D_{\sigma}(\bm{R}+\bm{x}_{\alpha})+\bm{d}-\bm{x}_{\tilde{\alpha}}$. Since $U(\sigma)$ acts on sublattice degrees of freedom, $U(\sigma)$ is bijective and independent of $\bm{R}$. Using the fact that $\bm{k}\cdot\bm{x}=D_{\sigma}\bm{k}\cdot D_{\sigma}\bm{x}$ which follows from $D_{\sigma}\in O(d)$, we obtain the symmetry transformation of $|\bm{k},\alpha\rangle$ under $\hat{\sigma}$:

[TABLE]

where we define $U_{\sigma}(\bm{k})=e^{-iD_{\sigma}\bm{k}\cdot\bm{d}}\,U(\sigma)$. From the condition that $\hat{H}$ is symmetric under $\hat{\sigma}$, $\hat{H}=\hat{\sigma}\hat{H}\hat{\sigma}^{-1}$, we obtain

[TABLE]

At high-symmetry momenta $\bm{k}^{\rm inv}_{\sigma}=D_{\sigma}\bm{k}^{\rm inv}_{\sigma}+\bm{G}(\bm{k}^{\rm inv}_{\sigma})$, $H(\bm{k}^{\rm inv}_{\sigma})$ and $V(\bm{G}(\bm{k}^{\rm inv}_{\sigma}))\,U(\bm{k}^{\rm inv}_{\sigma})$ commute,

[TABLE]

For inversion symmetry, the inversion-invariant momenta are given by $\bm{k}^{\rm inv}=-\bm{k}^{\rm inv}$ (mod $\bm{G}$).

From Eq. (163), we can find an energy eigenstate at $D_{\sigma}\bm{k}$, $|n,D_{\sigma}\bm{k}\rangle$, as a symmetry partner of $|n,\bm{k}\rangle$ up to unitary transformation $B_{\sigma}$, the sewing matrix:

[TABLE]

$B_{\sigma}$ appears here because a symmetry transformation can lead to a linear combination of energy eigenstates when they are degenerate. It is straightforward to show that the sewing matrix $B_{\sigma}(\bm{k})$,

[TABLE]

is a unitary matrix, $B_{\sigma}(\bm{k})\in U(n_{\rm occ})$, and is periodic in the Brillouin zone:

[TABLE]

Now, we investigate the transformation of the non-Abelian Berry connection $[a_{i}(\bm{k})]_{nm}$ defined by the occupied states $|n,\bm{k}\rangle$:

[TABLE]

where $\nabla_{i}$ denotes $\frac{\partial}{\partial k_{i}}$. Note that $a_{i}(\bm{k})$ is Hermitian.

[TABLE]

The Berry connection transforms under the gauge transformation $\lvert{n,\bm{k}}\rangle\rightarrow\lvert{n,\bm{k}}\rangle_{G}=\lvert{m,\bm{k}}\rangle\,G(\bm{k})_{mn}$ of the occupied states as

[TABLE]

From Eqs. (165) and (169), we obtain the symmetry transformation of $a_{i}(\bm{k})$ under $\hat{\sigma}$,

[TABLE]

where $\mathds{1}$ denotes an identity matrix in the occupied subspace. The results Eqs. (171) and (172) are used to derive the symmetry transformations of the Wilson loop.

H.3 Wilson line

Let us consider a line $\mathcal{L}:\bm{k}^{i}\rightarrow\bm{k}^{f}$ which is parameterized by $\lambda$. Start and end points $\bm{k}^{i,f}$ are given by $\bm{k}^{i}=\bm{k}(\lambda=0)$ and $\bm{k}^{f}=\bm{k}(\lambda=1)$ as shown in Fig. 15. Then, the Wilson line $[W_{\mathcal{L}}]_{nm}$ is defined by the path-ordered exponential of the Berry connection $[a_{i}(\bm{k})]_{nm}$ along the line $\mathcal{L}$:

[TABLE]

For convenience, we choose a smooth gauge along $\mathcal{L}$. Along the line $\mathcal{L}$, the covariant derivative $D_{\lambda}$ is defined as

[TABLE]

The Wilson line $W_{\mathcal{L}_{\lambda}}$ is the solution of

[TABLE]

with the boundary conditions $W_{\mathcal{L}_{\lambda}=0}=\mathds{1}$ and $W_{\mathcal{L}_{\lambda=1}}=W_{\mathcal{L}}$. Here, $\mathcal{L}_{\lambda}:\bm{k}^{i}\rightarrow\bm{k}(\lambda)$ and $\mathcal{L}_{\lambda=1}=\mathcal{L}$. With the help of Eqs. (H.3) and (175), various properties of the Wilson line are derived.

First, $W_{\mathcal{L}}$ is unitary: $W_{\mathcal{L}}^{\dagger}=W_{\mathcal{L}}^{-1}$. Also note that $W_{\mathcal{L}}^{\dagger}=W_{\mathcal{L}_{R}}$ where $W_{\mathcal{L}_{R}}$ is the “reversed” Wilson loop which is defined along the path-reversed line $\mathcal{L}_{R}:\bm{k}^{f}\rightarrow\bm{k}^{i}$.

Second, the gauge transformation of the Wilson line $W^{G}_{\mathcal{L}_{\lambda}}$ is derived from Eq. (171). From Eq. (175) for $W^{G}_{\mathcal{L}_{\lambda}}$,

[TABLE]

Hence,

[TABLE]

The solution of the above equation is given by

[TABLE]

with the boundary condition $W_{\mathcal{L}_{\lambda=0}}=\mathds{1}$.

Finally, we obtain the Wilson line $W_{\sigma\mathcal{L}_{\lambda}}$, which is the symmetry transformation of $W_{\mathcal{L}_{\lambda}}$ under $\hat{\sigma}$, where $\sigma\mathcal{L}_{\lambda}:D_{\sigma}\bm{k}^{i}\rightarrow D_{\sigma}\bm{k}^{f}$ as shown in Fig. 15. From Eqs. (172) and (175), we derive

[TABLE]

in a similar manner to the one used for deriving the gauge transformation. For completeness, we derive the Wilson line $W_{\Theta\mathcal{L}}$, which is the symmetry transformation of the Wilson line $W_{\mathcal{L}_{\lambda}}$ under antiunitary symmetry $\hat{\Theta}=\hat{\sigma}\mathcal{K}$, where $\mathcal{K}$ denotes the complex-conjugation operator. $\hat{\Theta}$ acts on $\bm{r}$ and $\bm{k}$ as $\hat{\Theta}:\bm{r}\rightarrow D_{\sigma}\bm{r}+\bm{d}$ and $\bm{k}\rightarrow-D_{\sigma}\bm{k}$. Thus, $\Theta{\mathcal{L}}:-D_{\sigma}\bm{k}^{i}\rightarrow-D_{\sigma}\bm{k}^{f}$ is a line with a start point $-D_{\sigma}\bm{k}^{i}$ and an end point $-D_{\sigma}\bm{k}^{f}$. The Wilson line $W_{\Theta\mathcal{L}}$ is given by

[TABLE]

where $B_{\Theta}(\bm{k})$ is sewing matrix for $\hat{\Theta}$:

[TABLE]

and $\hat{\Theta}\lvert{\bm{k}\,\alpha}\rangle=\lvert{-D_{\sigma}\bm{k}\,\beta}\rangle\,\Theta(\bm{k})_{\beta\alpha}$.

H.4 Wilson Loop

From Eq. (178), we conclude that the spectrum of $W_{\mathcal{L}}$ is gauge invariant when $\bm{k}^{i}=\bm{k}^{f}$ (mod $\bm{G}$). In the Brillouin zone, consider a non-contractible loop $\mathcal{L}:\bm{k}_{0}\rightarrow\bm{k}_{0}+\bm{G}$. For numerical calculations, we implement the discretized form,

[TABLE]

where $\bm{k}_{i}=\bm{k}_{0}+\frac{i}{N}\bm{G}$ and $P(\bm{k})=\sum_{n=1}^{n_{\rm occ}}\lvert{n,\bm{k}}\rangle\langle{n,\bm{k}}\rvert$ is the projector for the occupied states. However, we stick to formulate the Wilson loop in the continuum form

[TABLE]

to exploit the covariant derivative $D_{\lambda}$ and symmetry transformation of the Berry connection $a(\bm{k})$. The eigenstates $\lvert{\theta(\bm{k}_{0})}\rangle$ of the Wilson loop $W_{\bm{k}_{0}+\bm{G}\leftarrow\bm{k}_{0}}$ are given by $W_{\bm{k}_{0}+\bm{G}\leftarrow\bm{k}_{0}}\lvert{\theta(\bm{k}_{0})}\rangle=e^{i\theta(\bm{k}_{0})}\lvert{\theta(\bm{k}_{0})}\rangle$. Because of the periodic gauge, the Wilson loop is periodic,

[TABLE]

for any reciprocal vector $\bm{G}^{\prime}$. Also, the Wilson loop spectrum $\{\theta(\bm{k}_{0})\}$ is independent of the start point $\bm{k}_{0}$,

[TABLE]

where $\bm{k}_{0},\bm{k}^{\prime}_{0}\in\mathcal{L}:\bm{k}_{0}\rightarrow\bm{k}_{0}+\bm{G}$.

Now, let us discuss a consequence of symmetry $\hat{\sigma}$ on the Wilson loop. Although Eqs. (179) and (180) applies to arbitrary symmetries and lines, we restrict our attention to inversion $I$ and $k_{1}$-directed loop $\mathcal{L}_{1}$ in three dimensions, $\mathcal{L}_{1}:(k_{1,0},k_{2},k_{3})\rightarrow(k_{1,0}+2\pi,k_{2},k_{3})$ for the discussion relevant to this paper. Since inversion symmetry is symmorphic, one can find a coordinate system with $\bm{d}=0$. Applying Eq. (179),

[TABLE]

where $\bm{k}^{(1)}_{0}=(k_{1,0},k_{2},k_{3})$. With help of Eqs. (184) and (H.4), $W_{I\mathcal{L}_{1}}$ satisfies $W_{I\mathcal{L}_{1}}=W_{\mathcal{L}_{1}}^{-1}$ and its spectrum is $\{-\theta^{a}_{1}(-\bm{k}^{(1)}_{0})\}$. Also, Eq. (186) means that the spectrum of $W_{I\mathcal{L}_{1}}$, $\{-\theta^{a}_{1}(-\bm{k}^{(1)}_{0})\}$, is identical to that of $W_{\mathcal{L}_{1}}$, $\{\theta^{a}_{1}(\bm{k}^{(1)}_{0})\}$. Thus, $\lvert{\theta^{a}_{1}(\bm{k}^{(1)}_{0})}\rangle$ and $\lvert{-\theta^{a^{\prime}}_{1}(-\bm{k}^{(1)}_{0})}\rangle$ are paired by inversion symmetry, where $\lvert{\theta^{a}_{1}(\bm{k}^{(1)}_{0})}\rangle$ is an eigenstate of $W_{\mathcal{L}_{1}}$ with an eigenvalue $e^{i\theta^{a}_{1}(\bm{k}^{(1)}_{0})}$:

[TABLE]

This can be viewed as a unitary analogy of particle-hole symmetry. We define the sewing matrix $B^{(1)}_{I}(\bm{k}^{(1)}_{0})$ as

[TABLE]

Now, we define eigenstates $\lvert{u^{a}_{1}(\bm{k}^{(1)}_{0})}\rangle$ of $\hat{W}_{\mathcal{L}_{1}}=\lvert{n,\bm{k}^{(1)}_{0}}\rangle W_{\mathcal{L}_{1}}\langle{n,\bm{k}^{(1)}_{0}}\rvert$ with eigenvalue $e^{i\theta^{a}_{1}(\bm{k}^{(1)}_{0})}$, where

[TABLE]

H.5 Nested Wilson loop

When the Wilson bands are gapped, a nested Wilson loop can be defined Benalcazar et al. (2017a, b). We first choose an inversion-symmetric subset of the Wilson bands by defining a projector $\tilde{P}^{(1)}$ preserving inversion symmetry as shown in Fig. 16. The Wilson bands $\lvert{\theta^{a}_{1}(\bm{k}^{(1)}_{0})}\rangle$ selected by $\tilde{P}^{(1)}$ are labeled by $a=1,\dots,n_{\tilde{P}^{(1)}}$. To define the nested Wilson loop, we choose the loop $\mathcal{L}_{2}$ along the $k_{2}$ direction: $\mathcal{L}_{2}:(k_{1,0},k_{2,0},k_{3})\rightarrow(k_{1,0},k_{2,0}+2\pi,k_{3})$. In terms of the Berry connection $[a^{(1)}_{i}(\bm{k})]_{ab}$,

[TABLE]

the nested Wilson loop is defined by

[TABLE]

Repeating the derivation of the symmetry transformation of the Wilson loop, we obtain

[TABLE]

where $\bm{k}^{(2)}_{0}=(k_{1,0},k_{2,0},k_{3})$. Thus, $\lvert{\theta^{A}_{2\leftarrow 1}(\bm{k}^{(2)}_{0})}\rangle$ and $\lvert{-\theta^{A^{\prime}}_{2\leftarrow 1}(-\bm{k}^{(2)}_{0})}\rangle$ are paired by inversion symmetry where $\lvert{\theta^{A}_{2\leftarrow 1}(\bm{k}^{(2)}_{0})}\rangle$ is eigenstate of $W^{(1)}_{\mathcal{L}_{2}}$ with eigenvalue $e^{i\theta^{A}_{2\leftarrow 1}(\bm{k}^{(2)}_{0})}$:

[TABLE]

Thus, sewing matrix is defined as

[TABLE]

As for the nested Wilson loop, we define eigenstates $\lvert{u^{A}_{2\leftarrow 1}(\bm{k}^{(2)}_{0})}\rangle$ of $\hat{W}^{(1)}_{\mathcal{L}_{2}}=\lvert{n,\bm{k}^{(2)}_{0}}\rangle W^{(1)}_{\mathcal{L}_{2}}\langle{n,\bm{k}^{(2)}_{0}}\rvert$ with eigenvalue $e^{i\theta^{A}_{2\leftarrow 1}(\bm{k}^{(2)}_{0})}$, where

[TABLE]

When the nested Wilson bands are gapped, the doubly nested Wilson loop can be defined after we project the nested Wilson bands with an inversion-symmetric projector $\tilde{P}^{(2)}$. The nested Wilson bands $\lvert{\theta^{A}_{2\leftarrow 1}(\bm{k}^{(2)}_{0})}\rangle$ selected by $\tilde{P}^{(2)}$ are labeled by $A=1,\dots,n_{\tilde{P}^{(2)}}$. For the doubly nested Wilson loop, the loop $\mathcal{L}_{3}$ is chosen to be $k_{3}$-directed line, $\mathcal{L}_{3}:(k_{1,0},k_{2,0},k_{3,0})\rightarrow(k_{1,0},k_{2,0},k_{3,0}+2\pi)$. If we define a Berry connection $[a^{(2)}_{i}(\bm{k}^{(2)}_{0})]_{AB}$,

[TABLE]

the doubly nested Wilson loop can be defined by

[TABLE]

Repeating this procedure, the $l$th nested Wilson loop $W^{(l)}_{\mathcal{L}}$ can be defined for an integer.

H.6 Mapping between inversion parities and Wilson loop spectrum

Here, we review the mapping between parity configuration of the occupied states and the Wilson loop spectrum. We recommend to refer Alexandradinata et al. (2014) for the comprehensive study of the mapping. The mapping takes two steps.

Step 1: Mapping from parity configuration to Wilson loop spectrum.

Without loss of generality, we can focus on the $(k_{1},k_{2})$ plane in $d$D Brillouin zone. In this plane, inversion symmetry in two dimensions, $(k_{1},k_{2})\rightarrow(-k_{1},-k_{2})$, is induced by inversion symmetry in $d$ dimensions. Then, we can define the number of occupied states at $(\bar{k}_{1},\bar{k}_{2})=-(\bar{k}_{1},\bar{k}_{2})$ with parity $\xi$, $n_{\xi}(\bar{k}_{1},\bar{k}_{2})$. At $\bar{k}_{2}=0$ or $\pi$, one can define the fewest number $n_{s}(\bar{k}_{2})$ among $\{n_{\xi}(0,\bar{k}_{2}),n_{\xi}(\pi,\bar{k}_{2})\}$. When $n_{s}(k_{s})=n_{\xi_{s}}(k_{s},\bar{k}_{2})$, the mapping from $\{n_{\xi}(\bar{k}_{1},\bar{k}_{2})\}$ to the Wilson loop spectrum $\{\theta_{1}(\bar{k}_{2})\}$ is given by the following rules:

[TABLE]

where $n^{W}_{0}(\bar{k}_{2})$, $n^{W}_{\pi}(\bar{k}_{2})$ and $n^{W}_{c}(\bar{k}_{2})$ are the number of Wilson loop spectrum $\{\theta_{1}(\bar{k}_{2})\}$ with eigenvalues $\theta_{1}(\bar{k}_{2})=0$, $\pi$ and $(-\theta_{0},\theta_{0})$, respectively. Here, we denote $\theta_{0}$ as a generic value not equal to [math] or $\pi$.

Step 2: Mapping from Wilson loop spectrum to relative winding number.

With $\{n^{W}_{\theta}(\bar{k}_{2})\}=(n^{W}_{\xi_{s}}(\bar{k}_{2}),n^{W}_{-\xi_{s}}(\bar{k}_{2}),n^{W}_{c}(\bar{k}_{2}))$ for a given Wilson loop spectrum, $N_{d}$ is defined as

[TABLE]

When $|C|<N_{d}$ where $C$ is the Chern number defined in the $(k_{1},k_{2})$ plane, $N_{d}$ is given by $N_{d}=2N_{W}+|C|$. If $|C|\geq N_{d}$, the relative winding number $N_{W}$ is zero. Also note that the Chern number $C$ and $N_{d}$ satisfies

[TABLE]

Since a non-zero Chern number implies the stable band topology, let us assume that $C=0$. Then, the Wilson loop is unwound if $N_{d}\leq 0$. Otherwise, $N_{d}>0$, the relative winding number $N_{W}$ is given by $N_{W}=\frac{1}{2}N_{d}$.

Appendix I Characterization of fragile band topology in the nested Wilson loop spectrum

In this Appendix, we provide the proof for the claim that the nested Wilson loop spectrum characterized by a positive even integer $N_{d,{\rm nest}}$, which is defined as

[TABLE]

implies the fragile band topology. Here, $n^{W_{\rm nest}}_{\pi}(\bar{k}_{3})$ and $n^{W_{\rm nest}}_{c}(\bar{k}_{3})$ are the number of eigenvalues $\theta_{2\leftarrow 1}(\bar{k}_{3})$, which is equal to $\pi$ and takes a generic value, respectively, in the nested Wilson loop spectrum at $\bar{k}_{3}=0$ and $\pi$.

Before proving the claim, we note that the nested Wilson loop spectrum satisfies similar relations discussed in Appendix H.6 since the Wilson loop and the nested Wilson loop spectra are constrained by the inversion symmetry in the same way. First, we can assign an inversion parity to each Wilson loop eigenstate at inversion-invariant momenta $\bm{k}^{\rm inv}$. At $\bm{k}^{\rm inv}$, $|\theta_{1}(\bm{k}^{\rm inv})\rangle$ is an inversion eigenstate with definite parity when $\theta_{1}(\bm{k}^{\rm inv})=0$ and $\pi$. On the other hand, $|\theta_{1}(\bm{k}^{\rm inv})\rangle$ and $|-\theta_{1}(\bm{k}^{\rm inv})\rangle$ are paired by inversion symmetry when $\theta_{1}(\bm{k}^{\rm inv})$ is a generic value. This pair forms inversion-eigenstates with even and odd parities. The parities of the Wilson loop eigenstates at $\bm{k}^{\rm inv}$ are encoded in the sewing matrix $B^{(1)}_{I}(\bm{k}^{\rm inv})$ in Eq. (188). By repeating Step 1 in Appendix H.6, we obtain an eigenvalue configuration in the nested Wilson loop spectrum $\{n_{\theta}^{W_{\rm nest}}(\bar{k}_{3})\}=(n_{0}^{W_{\rm nest}}(\bar{k}_{3}),n_{\pi}^{W_{\rm nest}}(\bar{k}_{3}),n_{c}^{W_{\rm nest}}(\bar{k}_{3}))$ at $\bar{k}_{3}=0$ and $\pi$. This configuration uniquely determines $N_{d,{\rm nest}}$ in Eq. (I). Second, the nested Chern number $C_{\rm nest}$, the nested relative winding number $N_{W,{\rm nest}}$ and $N_{d,{\rm nest}}$ satisfy relations similar to ones for the Chern number $C$, the relative winding number $N_{W}$ and $N_{d}$ presented in Step 2 in Appendix H.6.

We begin our proof by showing that a band crossing between the Wilson bands does not change $N_{d,{\rm nest}}$ while it changes the nested Chern number $C_{\rm nest}$ only by an even integer. When a gap closing and reopening in the Wilson loop spectrum happens at a generic momentum, it does not change the parity configuration of the Wilson loop eigenstates. As the parity configuration determines an eigenvalue configuration in the nested Wilson loop spectrum, $N_{d,{\rm nest}}$ is invariant. If the Wilson bands are inverted at an inversion-invariant momentum $\bm{k}^{\rm inv}$, only pairs of Wilson bands $|\theta_{1}(\bm{k}^{\rm inv})\rangle$ and $|-\theta_{1}(\bm{k}^{\rm inv})\rangle$ are involved in the Wilson band inversion. Thus, there is no net exchange of inversion-parities, and $N_{d,{\rm nest}}$ is invariant. Consequently, the Wilson band crossing does not change $N_{d,{\rm nest}}$. On the other hand, the Wilson band crossing can change $C_{\rm nest}$ by an even integer since it always happens pairwise due to inversion symmetry Wieder and Bernevig (2018).

For an integer-valued $N_{d,{\rm nest}}$, there are three cases to be considered: $N_{d,{\rm nest}}$ is equal to i) an odd integer, ii) a positive even integer and iii) a negative even integer or 0.

1. i)

When $N_{d,{\rm nest}}$ is an odd integer, Eq. (200) implies that $C_{\rm nest}$ is also an odd integer. Since a Wilson band crossing cannot change an odd $C_{\rm nest}$ to zero, the nested Wilson loop spectrum indicates the stable band topology. 2. ii)

When $N_{d,{\rm nest}}$ is a positive even integer, Eq. (200) implies that $C_{\rm nest}$ is also an even integer. For even $C_{\rm nest}$, a Wilson band crossing can change $C_{\rm nest}$ to zero. However, the nested Wilson loop spectrum cannot be gapped and it shows the relative winding since $N_{d,{\rm nest}}=2N_{W,{\rm nest}}+|C_{\rm nest}|$ holds when $|C_{\rm nest}|<N_{d,{\rm nest}}$. 3. iii)

When $N_{d,{\rm nest}}$ is a negative even integer or zero, the nested relative winding number $N_{W,{\rm nest}}$ is always zero. This is because the relative winding exists only when $|C_{\rm nest}|<N_{d,{\rm nest}}$. Thus, the nested Wilson loop spectrum can be gapped when a Wilson band crossing happens.

Appendix J Numerical details of the Wilson loop calculation

In this appendix, we summarized the numerical details in Sec. VI.

J.1 $2$D third-order TI

The tight-binding Hamiltonian of $2$D third-order TI is given by

[TABLE]

where $\{\gamma_{1},\dots,\gamma_{5}\}=\{\tau_{x}\sigma_{x},\tau_{x}\sigma_{y},\tau_{x}\sigma_{z},\tau_{y}\sigma_{0},\tau_{z}\sigma_{0}\}$ and $\gamma_{5}=-\gamma_{1}\dots\gamma_{4}$. For numerical calculations, $\Delta_{1,3}=0.240$, $\Delta_{1,4}=-0.196$, $\Delta_{2,3}=-0.259$, and $\Delta_{2,4}=-0.204$ are used. $H_{{\rm A}(2,3)}(\bm{k})$ is inversion symmetric, i.e., $I_{0}\,H_{{\rm A}(2,3)}(\bm{k}){I_{0}}^{-1}=H_{{\rm A}(2,3)}(-\bm{k})$, with inversion symmetry operator $I_{0}$:

[TABLE]

The energy bands of $H_{{\rm A}(2,3)}(\bm{k})$ are shown in Fig. 17(b).

J.1.1 $2$D third-order TI $+$ $s(A)$

Now, we add an $s$ orbital localized at the Wyckoff position $A=(0,0)$. On-site energy $h_{1}$ and an interaction Hamiltonian between the original four bands and $s$ orbital are chosen as

[TABLE]

Thus, the perturbation Hamiltonian $H_{c}$ is given by

[TABLE]

where $0_{4}$ denotes $4\times 4$ null matrix.

Hence, Hamiltonian $H_{+s(A)}$ and inversion symmetry operator $I$ for the five bands are given by

[TABLE]

and $I\,H_{+s(A)}(\bm{k})\,I^{-1}=H_{+s(A)}(-\bm{k})$. The energy bands of $H_{+s(A)}(\bm{k})$ are shown in Fig. 17(c).

J.1.2 $2$D third-order TI $+$ $p(A)$

We add a $p$ orbital localized at the Wyckoff position $A=(0,0)$. On-site energy $h_{1}$ and an interaction Hamiltonian between the original four bands and $p$ orbital are chosen as

[TABLE]

Thus, the perturbation Hamiltonian $H_{c}$ is given by

[TABLE]

Hence, Hamiltonian $H_{+p(A)}$ and inversion symmetry operator $I$ for the five bands are given by

[TABLE]

and $I\,H_{+p(A)}(\bm{k})\,I^{-1}=H_{+p(A)}(-\bm{k})$. The energy bands of $H_{+p(A)}(\bm{k})$ are shown in Fig. 17(d).

J.2 $3$D fourth-order TI

The tight-binding Hamiltonian of $3$D fourth-order TI is given by

[TABLE]

where $\{\check{\gamma}_{1},\dots,\check{\gamma}_{7}\}=$ $\{\mu_{x}\tau_{x}\sigma_{x},\mu_{x}\tau_{x}\sigma_{y},\mu_{x}\tau_{x}\sigma_{z},\mu_{x}\tau_{y}\sigma_{0},$ $\mu_{x}\tau_{z}\sigma_{0},\mu_{y}\tau_{0}\sigma_{0},\mu_{z}\tau_{0}\sigma_{0}\}$ and $\check{\gamma}_{7}=i\check{\gamma}_{1}\dots\check{\gamma}_{6}$. For numerical calculations, $\check{\Delta}_{1,4}=-0.25$, $\check{\Delta}_{1,5}=0.05$, $\check{\Delta}_{1,6}=0.10$, $\check{\Delta}_{2,4}=-0.05$, $\check{\Delta}_{2,5}=0.30$, $\check{\Delta}_{2,6}=-0.10$, $\check{\Delta}_{3,4}=0.05$, $\check{\Delta}_{3,5}=-0.05$, and $\check{\Delta}_{3,6}=0.35$ are used. $H_{{\rm A}(3,4)}(\bm{k})$ is inversion symmetric, i.e., $I_{0}\,H_{{\rm A}(3,4)}(\bm{k}){I_{0}}^{-1}=H_{{\rm A}(3,4)}(-\bm{k})$, with inversion symmetry operator $I_{0}$,

[TABLE]

Energy bands of $H_{{\rm A}(3,4)}(\bm{k})$ are shown in Fig. 18(b).

J.2.1 $3$D fourth-order TI $+$ $3p(A)$

Now, we add three $p$ orbitals localized at the Wyckoff position $A=(0,0,0)$. For this, we add six orbitals at generic Wyckoff positions $\bm{x}_{9,\dots,14}$.

[TABLE]

Note that the original eight orbitals are localized at $\bm{x}_{1,\dots,8}=(0,0,0)$. The Hamiltonian for the additional six bands is chosen to be

[TABLE]

and it is invariant under the inversion symmetry operator $I_{c}$:

[TABLE]

Diagonalizing $h_{1}$ we obtain three $s$ and three $p$ orbitals localized at $A$. Three $s$ orbitals are degenerate with an energy eigenvalue $-24$, and three $p$ orbitals have energy eigenvalues $-14$, $-16$, and $-18$, respectively. Now, we couple the additional six bands to the original eight bands according to $h_{c}$:

[TABLE]

Thus, the perturbation Hamiltonian $H_{c1}$ is given by

[TABLE]

To sum up, the Hamiltonian $H$ and inversion symmetry operator $I$ for the 14 bands are given by

[TABLE]

and $I\,H(\bm{k})\,I^{-1}=H(-\bm{k})$. The energy bands of $H(\bm{k})$ are shown in Fig. 18(c). In order to study the effect of adding $p$ orbitals, we define projectors $P_{5}$, $P_{6}$, and $P_{7}$ as the projections to the occupied states with an additional $p(A)$ orbital, additional two $p(A)$ orbitals, and additional three $p(A)$ orbitals, respectively [see Fig. 18(c) for the definition of $P_{5,6,7}$].

J.2.2 Wilson band transition

The nested Wilson loop spectrum depends on the topology of Wilson bands. Thus, a bulk-gap-preserving perturbation can change the topology of the nested Wilson loop through a band crossing between Wilson bands, which we call a Wilson band transition. To investigate how the nested Wilson loop spectrum changes after Wilson band transition systematically, we choose a new coupling Hamiltonian ${h_{2}}^{\prime}$ as

[TABLE]

Thus, new perturbation Hamiltonian $H_{c2}$ is given by

[TABLE]

Now, we define an interpolating Hamiltonian $H_{c}(q)$ and $H(\bm{k},q)$ as

[TABLE]

and find critical values $q_{c}$ where the Wilson bands close the gap. The Wilson bands with the projector $P_{6}$ close the gap at $q_{c}=0.49$ and the band crossing points are $(k_{2},k_{3},\theta_{1}(k_{2},k_{3}))\simeq\pm(1.054\pi,0.0099\pi,-0.0117\pi)$ as shown in Figs. 19(c) and 19(d). For Wilson bands with the projector $P_{7}$, the Wilson band transition occurs at $q_{c}=0.52$ and the gap closes at $(k_{2},k_{3},\theta_{1}(k_{2},k_{3}))\simeq\pm(1.040\pi,0.0083\pi,-0.0114\pi)$. After the Wilson band transition, the nested Chern number changes from $\pm 2$ to [math] as shown in Fig. 13.