Higher-order topological phases in crystalline and non-crystalline systems: a review

Under mirror symmetry, $M_x: (x,y) \rightarrow (-x,y)$ or $M_y : (x,y) \rightarrow (x,-y)$, we see that the operator $\hat{Q}_{xy}$ in equation (12) transforms into $-\hat{Q}_{xy}$. Thus, we have $q_{xy} \rightarrow -q_{xy}$. If the system is invariant under the mirror transformation, we have $q_{xy} = -q_{xy} \mod 1$ so that q_xy takes the quantized value of either 0 or $1/2$. Similarly, we see that the four-fold rotation $C_4 : (x,y) \rightarrow (-y,x)$ also flips the sign of q_xy . Thus, the C₄ symmetry enforces that $q_{xy} = -q_{xy} \mod 1$, indicating that q_xy is quantized in the presence of the C₄ symmetry.

Apart from crystalline symmetries, the quadrupole moment can be quantized by internal symmetries alone. In the following, we present the proof of quantization of q_xy (13) protected by chiral symmetry by following [134, 135]. Let H be a generic single-particle Hamiltonian in real space preserving chiral (sublattice) symmetry, that is, there exists a unitury operator Π such that ${\Pi} H {\Pi}^{-1} = -H$. If we consider an occupied state $|\psi_n\rangle$ with energy E_n , then we can obtain an unoccupied state by applying Π to the occupied one, that is, ${\Pi}|\psi_n\rangle$ with energy $-E_n$. Since the energy spectrum is symmetric about the zero energy, we consider the half-filling case with $n_a = 2 n_c$. Then, we have $q_{xy}^{(0)} = \frac{1}{2}\sum_{j = 1}^{n_a} x_j y_j/(L_x L_y) = \frac{1}{4\pi}\mathrm{Im}\log \det \hat{D} \mod 1$. The set $\{{\Pi}|\psi_1\rangle,{\Pi}|\psi_2\rangle,\cdots, {\Pi}|\psi_{n_c}\rangle\}$ therefore constitutes the unoccupied states. These unoccupied states can be represented by $U_{u} = \left({\Pi}|\psi_1\rangle,{\Pi}|\psi_2\rangle,\cdots, {\Pi}|\psi_{n_c}\rangle\right) = {\Pi}U_{o}$. The quadrupole moment (13) can be written as

$$\begin{equation} q_{xy} = \frac{1}{2\pi}\mathrm{Im}\log \left[\det \left(U_{o}^\dagger \hat{D} U_{o} \right) \sqrt{\det \hat{D}^\dagger}\right]. \end{equation} \tag{ 14 }$$

Because the unitary matrix $\hat{D}$ commutes with the chiral (sublattice) symmetry transformation Π with $\Pi^\dagger \Pi = 1$, i.e. $\left[\hat{D},\Pi\right] = 0$,

$$\begin{align} q_{xy}& = \frac{1}{2\pi}\mathrm{Im}\log \left[\det \left(U_{u}^\dagger \Pi \hat{D} \Pi^\dagger U_{u} \right) \sqrt{\det \hat{D}^\dagger}\right] \nonumber \\ & = \frac{1}{2\pi}\mathrm{Im}\log \left[\det \left(U_{u}^\dagger \hat{D} U_{u} \right) \sqrt{\det \hat{D}^\dagger}\right]. \end{align} \tag{ 15 }$$

Next, we will prove that $2q_{xy} = 0 \mod 1$.

Proof. We define an $n_a \times n_a$ unitary matrix $U_{t} = \left(U_{o},U_{u}\right)$. It can be easily seen that $\det \hat{D}^\dagger = \det \left(U_t^\dagger \hat{D}^\dagger U_t \right)$. Then, we have the following relations

$$\begin{align} 4\pi q_{xy}& = \mathrm{Im}\log \left[\det \left(U_{o}^\dagger \hat{D} U_{o} \right) \sqrt{\det \hat{D}^\dagger}\right] \nonumber \\ &\quad + \mathrm{Im}\log \left[\det \left(U_{u}^\dagger \hat{D} U_{u} \right) \sqrt{\det \hat{D}^\dagger}\right] \nonumber \\ & = \mathrm{Im}\log \left[ \det \left(U_{o}^\dagger \hat{D} U_{o} \right) \det \left(U_{u}^\dagger \hat{D} U_{u}\right)\right] \nonumber \\ &\quad +\mathrm{Im}\log \det \left(U_t^\dagger \hat{D}^\dagger U_t \right) \nonumber \\ & = \mathrm{Im}\log \det \begin{pmatrix} U_{o}^\dagger \hat{D} U_{o} & U_{o}^\dagger \hat{D} U_{u} \\ 0 & U_{u}^\dagger \hat{D} U_{u} \\ \end{pmatrix} \nonumber \\ &\quad +\mathrm{Im}\log \det \begin{pmatrix} U_{o}^\dagger \hat{D}^\dagger U_{o} & U_{o}^\dagger \hat{D}^\dagger U_{u} \\ U_{u}^\dagger \hat{D}^\dagger U_{o} & U_{u}^\dagger \hat{D}^\dagger U_{u} \\ \end{pmatrix} \nonumber \\ & = \mathrm{Im}\log \det \begin{pmatrix} \mathbb{I} & 0 \\ U_{u}^\dagger \hat{D} U_{u} U_{u}^\dagger \hat{D}^\dagger U_{o} & U_{u}^\dagger \hat{D} U_{u} U_{u}^\dagger \hat{D}^\dagger U_{u} \\ \end{pmatrix} \nonumber \\ & = \mathrm{Im}\log \det \left(U_{u}^\dagger \hat{D} U_{u} U_{u}^\dagger \hat{D}^\dagger U_{u} \right) \nonumber \\ & = \mathrm{Im}\log \det \left(U_{u}^\dagger \hat{D} U_{u} \right) + \mathrm{Im}\log \left[ \det \left(U_{u}^\dagger \hat{D} U_{u} \right)\right]^{*} \nonumber \\ & = 0 \mod 2\pi. \end{align} \tag{ 16 }$$

In the derivation, the following orthonormal properties are used, $U_{o}^\dagger U_{o} = U_{u}^\dagger U_{u} = \mathbb{I}$, $U_{o}^\dagger U_{u} = 0$ and $U_{o}U_{o}^\dagger + U_{u}U_{u}^\dagger = \mathbb{I}$. □

Therefore, we get the conclusion that $2q_{xy} = 0 \mod 1$, namely, q_xy is quantized to 0 or $1/2$ up to an integer. One can also apply the above procedure to prove that the octupole moment is quantized in 3Ds in the presence of chiral symmetry. Similar conclusion can be generalized to systems with particle-hole symmetry [135]. Since the chiral symmetry alone can protect the quantization of the quadrupole moment, we can have quadrupole TIs in disordered systems lacking crystalline symmetry [134, 135, 226].

In general, we consider a Bloch Hamiltonian $H(\textbf{k})$ which has $N_{\mathrm{orb}}$ degrees of freedom and $N_{\mathrm{occ}}$ occupied bands with Bloch functions $|u_\textbf{k}^{n}\rangle$ for $n = 1,\cdots,N_{\mathrm{occ}}$. The Wilson loop is closely related to the position operator projected into the occupied space [254]. In the following, we will focus on the Wilson loop along k_x , $\mathcal{W}_{x, \textbf{k}}$, and the Wilson loop $\mathcal{W}_{y, {\textbf{k}}}$ along k_y can be evaluated similarly. The Wilson loop $\mathcal{W}_{x, \textbf{k}}$ starting from the base point $\textbf{k} = (k_x,k_y)$ of the loop in Brillouin zone, is defined under full periodic boundary conditions as follows

$$\begin{equation} \mathcal{W}_{x, {\textbf{k}}} = F_{x,\textbf{k}+\left(N_x-1\right)\mathbf{\delta k_x}} \cdots F_{x,\textbf{k} + \mathbf{\delta k_x}} F_{x,{\textbf{k}}}, \end{equation} \tag{ 17 }$$

where $[F_{x,{\textbf{k}}}]^{mn} = \left\langle u^{m}_\mathbf{k+{\delta} k_x} \right| \left. u^n_\textbf{k} \right\rangle$, for $\mathbf{\delta k_x} = (2\pi/N_x,0)$. Since the Wilson loop is unitary in the thermodynamic limit, it can be defined as an exponential of a Hermitian matrix

$$\begin{equation} \mathcal{W}_{x, {\textbf{k}}}\equiv e^{iH_{\mathcal{W}_x}\left({\textbf{k}}\right)}. \end{equation} \tag{ 18 }$$

Then, we refer to $H_{\mathcal{W}_x}({\textbf{k}})$ as the Wannier Hamiltonian. It has been illustrated that the Wannier Hamiltonian's spectrum has a close connection with the spectrum and topology of the boundary perpendicular to the x direction [255].

The Wilson loop operator can be diagonalized as

$$\begin{equation} \mathcal{W}_{x, {\textbf{k}}} |\nu^{\,j}_{x,{\textbf{k}}}\rangle = e^{i 2\pi \nu^{\,j}_x\left(k_y\right)} |\nu^{\,j}_{x,{\textbf{k}}}\rangle, \end{equation} \tag{ 19 }$$

where $|\nu^{\,j}_{x,{\textbf{k}}}\rangle$ is the eigenvectors with components $[\nu^{\,j}_{x,{\textbf{k}}}]^n$ for $n = 1,\cdots,N_{\mathrm{occ}}$. The Wannier Hamiltonian $H_{\mathcal{W}_{x}}(\textbf{k})$ has eigenvalues $2\pi \nu^{\,j}_{x}(k_y)$, $j = 1,\cdots, N_{\mathrm{occ}}$, which only depend on k_y of the base point k. The eigenvalues $\nu^{\,j}_{x}(k_y)$ defined modulo 1 are referred to as the Wannier centers and compose the Wannier bands.

For the BBH model with two occupied bands, the phases $\nu^{j = \pm}_x(k_y)$ of the Wilson loop $\mathcal{W}_{x, \textbf{k}}$ over two occupied bands form a pair of opposite values $\nu_x^-(k_y) = -\nu_x^+(k_y) \mod 1$ due to the mirror symmetry M_x , leading to the vanishing of the bulk polarizations. The two Wannier bands $\nu^{j = \pm}_x(k_y)$ have a Wannier gap at both $\nu_x = 0$ and $\nu_x = 1/2$ over the Brillouin zone. Therefore, we can have a 1D gapped Wannier Hamiltonian $H_{\mathcal{W}_x}(k_y)$ for $k_y \in [-\pi,\pi]$, which may host nontrivial band topology.

When the Wannier bands $\nu^{\,j}_x(k_y)$ are gapped for $k_y \in (-\pi,\pi]$, we can define two Wannier sectors separated by the gap. For quadrupole insulators with two anticommuting mirror symmetries M_x and M_y , the Wannier bands $\nu^{\,j}_x(k_y)$ are gapped for $k_y \in (-\pi,\pi]$ at both $\nu_x = 0$ and $\nu_x = 1/2$. Then, we can define two Wannier sectors $\nu^-_x \in (0,1/2)$ and $\nu^+_x \in (1/2,1)$. In the BBH case, there is only one band in each Wannier sector.

The topological invariant for a Wannier sector ν_x can be evaluated based on the nested Wilson loops. With the Wannier bands $\nu^{l}_x(k_y)$ and Wilson loop eigenstates $|\nu^l_{x,{\textbf{k}}}\rangle$, the nested Wilson loop is a Wilson loop of one subspace of Wannier bands along k_y . Specifically, for the two Wannier sectors $\nu^{\pm}_x$ split by the Wannier gap, we define the Wannier band basis in the Wannier sector $\nu_x^-$ as

$$\begin{equation} |w^l_{x,{\textbf{k}}}\rangle = \sum_{n = 1}^{N_{\mathrm{occ}}}|u^n_\textbf{k}\rangle \left[\nu^l_{x,{\textbf{k}}}\right]^n \end{equation} \tag{ 20 }$$

for $l \in 1\ldots N_W$. N_W is the number of the Wannier bands in the sector $\nu_x^-$. This basis obeys $\left\langle w^l_{x,{\textbf{k}}} \right| \left. w^{l^{^{\prime}}}_{x, {\textbf{k}}} \right\rangle = \delta_{l l^{^{\prime}}}$. Then, the nested Wilson loop along y for the Wannier sector ν_x over the Wannier band basis is defined as

$$\begin{equation} \tilde{\mathcal{W}}^{\nu_x}_{y, k_x} = F^{\nu_x}_{y,\textbf{k}+\left(N_y-1\right)\mathbf{\delta k_y}} \cdots F^{\nu_x}_{y,\textbf{k} + \mathbf{\delta k_y}} F^{\nu_x}_{y,{\textbf{k}}}, \end{equation} \tag{ 21 }$$

where $[F^{\nu_x}_{y,{\textbf{k}}}]^{rs} = \left\langle w^r_{x,\textbf{k}+\mathbf{\delta k_y}} \right| \left. w^s_{x,{\textbf{k}}} \right\rangle$ with $r,s \in 1 \ldots N_W$ over all Wannier bands in the Wannier sector ν_x , and $\mathbf{\delta k_y} = (0, 2\pi/N_y)$. The total polarization of the Wannier sector ν_x is

$$\begin{equation} p^{\nu^\pm_x}_y = \frac{1}{N_x}\sum_{k_x} \frac{1}{2\pi}\mathrm{Im}\log \det \left[\tilde{\mathcal{W}}^\pm_{y, k_x} \right], \end{equation} \tag{ 22 }$$

which, in the thermodynamic limit, becomes an integral

$$\begin{equation} p^{\nu_x}_y = -\frac{1}{\left(2\pi\right)^2} \int_{\textrm{BZ} } {\text{Tr}}\, \left[\tilde{\mathcal{A}}^{\nu_x}_{y,{\textbf{k}}}\right] \textrm{d}^2{\textbf{k}}, \end{equation} \tag{ 23 }$$

where $\tilde{\mathcal{A}}^{\nu_x}_{y,{\textbf{k}}}$ is the Berry connection of Wannier bands ν_x having components

$$\begin{align} \left[\tilde{\mathcal{A}}^{\nu_x}_{y,{\textbf{k}}}\right]^{l l^{^{\prime}}} = -i \langle w^l_{x,{\textbf{k}}} | \partial_{k_y} |w^{l^{^{\prime}}}_{x,{\textbf{k}}}\rangle, \end{align} \tag{ 24 }$$

where $l,l^{^{\prime}} \in 1 \ldots N_W$ run over the Wannier bands in Wannier sector ν_x [34, 35].

Under the mirror reflection M_y and M_x , the Wannier-sector polarizations satisfy $p^{\nu_x}_y \equiv -p^{\nu_x}_y \mod 1$ and $p^{\nu_y}_x \equiv -p^{\nu_y}_x \mod 1$, respectively [34, 35]. Hence, $p^{\nu_x}_y$ and $p^{\nu_y}_x$ are quantized to 0 or $1/2$. For the topological nontrivial phase of the BBH model (2), we have $(p^{\nu_x}_y, p^{\nu_y}_x) = (1/2,1/2)$ in consistent with the quadrupole moment $q_{xy} = 1/2$. The transition between the phase with different Wannier-sector polarizations is driven by the gap closing of the Wannier bands.

The mirror symmetries can also simplify the computation of the topological invariants. Since the 1D Wannier Hamiltonian $H_{\mathcal{W}_{x}}(\textbf{k})$ is invariant under the mirror reflection M_y inherited from the bulk, the Wannier-sector polarization can be obtained from the reflection eigenvalues of Wannier band basis (20) at mirror invariant momenta [34, 35].

Consider a generic Bloch Hamiltonian $H(\textbf{k})$ with a symmorphic lattice symmetry g: $\textbf{r}\rightarrow D_g\textbf{r}$, described by $g H(\textbf{k}) g^\dagger = H(D_g\textbf{k})$, where g is a unitary operator for the symmetry. We now define a Wilson line following a path C in momentum space from k_i to k_f as

$$\begin{equation} \mathcal{W}_{\textbf{k}_f \leftarrow \textbf{k}_i} = F_{\textbf{k}_f-\delta \textbf{k}}\cdots F_{\textbf{k}_i+\delta \textbf{k}}F_{\textbf{k}_i}, \end{equation} \tag{ 25 }$$

where $[F_\textbf{k}]^{mn} = \langle u_\mathbf{k+\delta k}^m|u_\textbf{k}^n\rangle$ with $m,n$ labeling the occupied bands and $\delta\textbf{k} = (\textbf{k}_f - \textbf{k}_i)/N$. We take $N\rightarrow \infty$ in the thermodynamic limit. The Wilson line operator transforms under the symmetry as [34, 35]

$$\begin{equation} B_{g,\textbf{k}_f} \mathcal{W}_{\textbf{k}_f \leftarrow \textbf{k}_i} B_{g,\textbf{k}_i}^\dagger = \mathcal{W}_{D_g \textbf{k}_f \leftarrow D_g \textbf{k}_i}. \end{equation} \tag{ 26 }$$

Here a unitary sewing matrix $B_{g,\textbf{k}}^{mn} = \langle u_{D_g \textbf{k}}^m| g |u_\textbf{k}^n\rangle$ is defined. The matrix reflects the connection between eigenstates at k with eigenstates at $D_g \textbf{k}$. In particular, a Wilson loop along the contour $\mathcal{C}$ at the base point k satisfies

$$\begin{equation} B_{g,\textbf{k}} \mathcal{W}_{\mathcal{C},\textbf{k}} B_{g,\textbf{k}}^\dagger = \mathcal{W}_{D_g \mathcal{C},D_g {\textbf{k}}} . \end{equation} \tag{ 27 }$$

Now we focus on the topology of the Wannier band $\nu_x(k_y)$ from the Wilson loop $\mathcal{W}_{x,\textbf{k}}$ in mirror-symmetric systems. The results for $\nu_y(k_y)$ can be obtained analogously. Under the mirror symmetry M_x , we have the transformation for the Wilson loop as

$$\begin{equation} B_{m_x,\textbf{k}_f} \mathcal{W}_{x,\textbf{k}} B_{m_x,\textbf{k}_i}^\dagger = \mathcal{W}_{-x,M_x\textbf{k}}, \end{equation} \tag{ 28 }$$

where x and −x label the direction that a Wilson loop is obtained. We then define a Wannier Hamiltonian for the Wilson loop $\mathcal{W}_{x,(k_x = 0,k_y)}$. The Wannier Hamiltonian in equation (18) is a multivalued function due to the logarithm, and thus we need to redefine it relative to a branch cut ε [67],

$$\begin{equation} H_{\mathcal{W}_x}^\epsilon \left(k_y\right) \equiv -i\log_{\epsilon}\left(\mathcal{W}_{x,\textbf{k}}\right), \end{equation} \tag{ 29 }$$

where we take $\log_{\epsilon} e^{i\phi} = i\phi$, for $\epsilon\unicode{x2A7D} \phi\lt\epsilon+2\pi$. With the symmetries, it is proved that the Wannier Hamiltonian should satisfy the following constraints [67],

$$\begin{align} &B_{m_x,\left(0,k_y\right)}H_{\mathcal{W}_x}^\epsilon \left(k_y\right)B_{m_x,\left(0,k_y\right)}^\dagger = -H_{\mathcal{W}_x}^{-\epsilon}\left(k_y\right) +2\pi I_{N_{\textrm{occ}}} \nonumber \\ &B_{m_x,\left(0,k_y\right)}H_{\mathcal{W}_x}^\epsilon \left(k_y\right)B_{m_x,\left(0,k_y\right)}^\dagger = -H_{\mathcal{W}_x}^{-\epsilon+2\pi}\left(k_y\right) +4\pi I_{N_{\textrm{occ}}}. \end{align} \tag{ 30 }$$

The Wannier bands obey the periodicity from the relation $\nu_x \equiv \nu_x \mod 1$, reminiscent of the effective Hamiltonian of Floquet TIs [256, 257]. It is thus possible that the Wannier bands under open boundary exhibit in-gap modes at both $\nu_x = 0$ and $\nu_x = 1/2$ in a cylinder geometry. As a result, the Wannier-sector polarization are trivial so that they fail to identify the topology of the Wannier bands. The quadrupole TIs with this type of Wannier bands are called anomalous quadrupole TIs [43]. They resemble the anomalous topological phases in a periodically driven system, where the edge states persist even though the traditional topological invariant of the effective Hamiltonian vanishes [256].

To correctly identify the topology of Wannier bands, Yang and coworkers introduce a topological invariant which can fully characterize the change of the Wannier band topology due to the Wannier gap closing [67]. Specifically, a Wilson line with respect to ε is introduced,

$$\begin{equation} \mathcal{W}_{x,k_x\leftarrow 0}^{\epsilon}\left(k_y\right) \equiv \mathcal{W}_{x,k_x\leftarrow 0}\left(k_y\right) \exp\left[-iH_{\mathcal{W}_x}^{\epsilon}\left(k_y\right)\frac{k_x}{2\pi}\right], \end{equation} \tag{ 31 }$$

where $\mathcal{W}_{x,k_x\leftarrow 0}(k_y)\equiv \mathcal{W}_{(k_x,k_y)\leftarrow (0,k_y)}$. It is shown that $\mathcal{W}_{x,2\pi\leftarrow 0}^{\epsilon}(k_y) = I_{N_{\textrm{occ}}}$, and one arrives at the fact that $\mathcal{W}_{x,k_x\leftarrow 0}^{\epsilon}(k_y)$ is periodic versus k_x . The symmetry further imposes constraints on the Wilson line with respect to ε for $k_x = \pi$ and $\epsilon = 0,\pi$ such that [67]

$$\begin{align} &B_{m_x,\left(\pi,k_y\right)}\mathcal{W}_{x,\pi \leftarrow 0}^{\epsilon = 0}\left(k_y\right) B_{m_x,\left(0,k_y\right)}^{\dagger} = -\mathcal{W}_{\pi\leftarrow 0}^{\epsilon = 0}\left(k_y\right), \nonumber \\ &B_{m_x,\left(\pi,k_y\right)}\mathcal{W}_{x,\pi \leftarrow 0}^{\epsilon = \pi}\left(k_y\right) B_{m_x,\left(0,k_y\right)}^{\dagger} = \mathcal{W}_{\pi\leftarrow 0}^{\epsilon = \pi}\left(k_y\right). \end{align} \tag{ 32 }$$

For the quadrupole insulator respecting two anticommuting mirror symmetries with $N_{\textrm{occ}} = 2$, it can be shown that there exists a basis in which the sewing matrices along the mirror invariant lines take the form [67]

$$\begin{equation} B_{m_x,\left(\pi,k_y\right)} = B_{m_x,\left(0,k_y\right)} = \left( \begin{array}{cc} 1 & 0 \\ 0 & -1 \\ \end{array} \right). \end{equation} \tag{ 33 }$$

As a result, the Wilson loops are $2\times 2$ matrices taking forms as

$$\begin{align} \mathcal{W}_{x,\pi \leftarrow 0}^{\epsilon = 0}\left(k_y\right)& = \left( \begin{array}{cc} 0 & U_+^{\epsilon = 0}\left(k_y\right) \\ U_-^{\epsilon = 0}\left(k_y\right) & 0 \\ \end{array} \right), \end{align} \tag{ 34 }$$

$$\begin{align} \mathcal{W}_{x,\pi \leftarrow 0}^{\epsilon = \pi}\left(k_y\right)& = \left( \begin{array}{cc} U_+^{\epsilon = \pi}\left(k_y\right) & 0 \\ 0 & U_-^{\epsilon = \pi}\left(k_y\right) \\ \end{array} \right). \end{align} \tag{ 35 }$$

The topological invariant for the Wannier gaps at $\epsilon = 0,\pi$ can thus be defined as a winding number

$$\begin{equation} W_{\nu_x}^{\epsilon} = \frac{1}{2\pi i}\int_0^{2\pi} \textrm{d}k_y\partial_{k_y}\log U_{+}^{\epsilon}\left(k_y\right). \end{equation} \tag{ 36 }$$

Since the Wilson line depends on the eigenstates, the winding number is not gauge invariant. Thus, when one varies the system parameter to observe the topological phase transitions changing the winding number, the global phase of wave functions are required to be continuous along mirror invariant lines in Brillouin zone. The numerical procedure is described in detail in [67].

It has been shown that the change of the winding number can reflect the Wannier gap closing at either ν = 0 or $\nu = \pm 1/2$. The winding number only responds to the closure of the bulk energy gap or the Wannier gap rather than the edge energy gap closing.

The edge polarization $p_{x}^{\mathrm{edge~}\beta}$ ($p_{y}^{\mathrm{edge~}\alpha}$) describes the edge polarization per unit length along the x (y) direction at the boundaries perpendicular to y (x). The Greek letters $\beta = \pm y$ ($\alpha = \pm x$) label the y-normal (x-normal) edges with the sign identifying the top or bottom (left or right) boundaries. We now review how the edge polarization is calculated based on the Wilson loop in a cylinder geometry following [34, 35]. Specifically, for a 2D insulator with $N_x\times N_y$ unit cells, we evaluate the polarization p_x along x with respect to the position along y. We take the periodic boundary condition along x and the open boundary condition along y. The system can be treated as a pseudo-1D Hamiltonian with $N_{\mathrm{orb}} N_y$ internal orbitals and N_x unit cells along x, and k_x remains a good quantum number. We can compute the Wilson loop $\mathcal{W}_x$ for the pseudo-1D system with $N_{\mathrm{occ}} N_y$ occupied bands.

With the obtained Wilson loop of the pseudo-1D system, we can construct the hybrid Wannier functions as [35]

$$\begin{align} |\Psi^j_{R_x}\rangle = \frac{1}{\sqrt{N_x}} \sum_{n}\sum_{k_x} \left[ \nu^{\,j}_{k_x} \right]^n e^{-i k_x R_x} |n,k_x\rangle, \end{align} \tag{ 37 }$$

where $\left[ \nu^{\,j}_{k_x} \right]^n$ is the nth component of the jth eigenstate $|\nu^{\,j}_{k_x}\rangle$ for the Wilson loop with the Wannier center $\nu^{\,j}_{x}$ for $j = 1,\cdots,N_{\mathrm{occ}}N_y$, and $|n,k_x\rangle$ is the nth eigenstates of the pseudo-1D Hamiltonian at k_x . One can calculate the density of the hybrid Wannier functions (37) to extract the spatially resolved polarization [35]. The density is given by

$$\begin{equation} \rho^{\,j}\left(R_y\right) = \frac{1}{N_x} \sum_{n, k_x, \alpha} \left| \left[u^n_{k_x}\right]^{R_y, \alpha}\left[\nu^{\,j}_{k_x}\right]^n\right|^2, \end{equation} \tag{ 38 }$$

where $[u^n_{k_x}]^{R_y, \alpha}$ is the corresponding component of the nth eigenstate of the pseudo-1D Hamiltonian. Then, we can calculate the polarization along x at each site R_y as [35]

$$\begin{align} p_x\left(R_y\right) = \sum_j \rho^{\,j}\left(R_y\right) \nu_x^j. \end{align} \tag{ 39 }$$

In the calculation, we need to break the mirror symmetries infinitesimally to split the degeneracy of two Wannier centers at $\pm 1/2$ so that they gives opposite polarizations localized at opposite edges. One then can compute the edge polarization by summing $p_x(R_y)$ over a half of the system along y, i.e.

$$\begin{align} p_x^{\mathrm{edge~}-y}& = \sum_{R_y = 1}^{N_y/2}p_x\left(R_y\right) \nonumber\\ & = -p_x^{\mathrm{edge~}+y} = -\sum_{R_y = N_y/2+1}^{N_y}p_x\left(R_y\right). \end{align} \tag{ 40 }$$

It is quantized to 0 or 0.5 up to an integer. One can similarly calculate the edge polarization along y.

Because of the mirror symmetry, the nontrivial Wannier bands can have mid-gap edge modes at 0 or $\pm 1/2$. However, only the edge modes at $1/2$ are responsible for nontrivial edge polarizations. Specifically, if the Wannier bands ν_x under open boundaries perpendicular to y have edge modes at $\pm 1/2$, these modes give rise to the edge-localized polarizations along x. Thus, the winding number in equation (36) calculated based on bulk wave functions can characterize the change of the corresponding edge polarization due to bulk gap closing or Wannier band gap closing. However, the edge polarization will also change when the edge energy gap closes and reopens. The edge energy gap closure is usually associated with the change of the quadrupole moment. If the correspondence is true, then one can take into account the effects of the edge energy gap closure by defining the change of edge polarization (a $\mathbb{Z}_2$ quantity) with respect to some system parameter γ as [67]

$$\begin{align} p_x^{\mathrm{edge}}\left(\gamma_1\right)-p_x^{\mathrm{edge}}\left(\gamma_0\right) &= e \left\{\left[W_{\nu_x}^{\epsilon = \pi}\left(\gamma_1\right)-W_{\nu_x}^{\epsilon = \pi}\left(\gamma_0\right) \right. \right. \nonumber \\ &\quad \left.\left. -\Delta N_{q,x} \right]/2 \mod 1 \right\}. \end{align} \tag{ 41 }$$

Here one uses $\Delta N_{q,x}$ to count the number of times that the change of the quadrupole moment occurs from the y-normal edge energy gap closing, as one varies the system parameter γ from γ₀ to γ₁. Due to the fact that the right sides only depend on the bulk property, one may conclude that the bulk topolgy determines the edge polarization. In a system with mirror symmetries, one can also evaluate $\Delta N_{q,x}$ by counting the number of times of the change in a parity (the eigenvalue of the reflection symmetry operator along x) at the high-symmetric points $k_x = 0$ or $k_x = \pi$ for a state localized at a y-normal boundary. If the system is in a trivial phase with $p_x^{\mathrm{edge}}(\gamma_0) = 0$ at γ₀, the formula in equation (41) can be reduced to

$$\begin{equation} p_x^{\mathrm{edge}}\left(\gamma_1\right) = e \left\{\left[W_{\nu_x}^{\epsilon = \pi}\left(\gamma_1\right) -\Delta N_{q,x} \right]/2 \mod 1 \right\} \end{equation} \tag{ 42 }$$

by choosing a gauge such that $W_{\nu_x}^{\epsilon = \pi}(\gamma_0) = 0$ in the topologically trivial phase. It has been shown that the edge polarization calculated based on the formula (42) coincides with the Wannier spectrum and the edge polarization computed using the hybrid Wannier functions in a cylinder geometry [67].

This is the simplest case beyond the BBH model. Figure 4 illustrates the energy spectrum of the lattice model in a geometry with open boundaries along both directions with respect to φ. We clearly see the appearance of four zero-energy states in the gap highlighted by the blue lines when $\phi \in (-\frac{\pi}{2}, -\frac{\pi}{4})$ and $(\frac{\pi}{4}, \frac{\pi}{2})$. These states are spatially localized at the four corners of the 2D lattice [118]. As we have demonstrated in the previous section, the presence of zero-energy corner states is resulted from the topology of the 1D off-diagonal AAH models. For $H_{\textrm{AAH},x}(k_x)$ with $\alpha_x = 1/2$ and $H_{\textrm{AAH},y}(k_y)$ with $\alpha_y = 1/4$, zero-energy edge states appear when $\phi \in (-\frac{\pi}{2}, \frac{\pi}{2})$ and $\phi \in (-\frac{\pi}{2},-\frac{\pi}{4}) \cup (\frac{\pi}{4}, \frac{\pi}{2})$, respectively. The regions' intersection thus corresponds to the topological parameter region observed in figure 4(a), corresponding to the topological invariant $\nu_{2D} = 1$ with $w_x = w_y = 1$.

Zoom In Zoom Out Reset image size

Since the system respects chiral symmetry, one can also use the quadrupole moment to characterize its topology given that there are four zero-energy states (one at each corner). Note that the quadrupole moment is a $\mathbb{Z}_2$ topological invariant (similar to the Berry phase), and thus can only distinguish between the phases with even (including zero) number of corner modes at each corner and the ones with old number. In figure 4(b), we see that the quadrupole moment is equal to $1/2$ in the topological region and vanishes in the trivial one. It changes through the y-normal edge gap closing at $\phi = \pm \pi/2$ and the x-normal edge gap closing at $\phi = \pm \pi/4$.

If one views φ as the momentum k_z along z, then one obtains a 3D AAH model, which exhibits a 3D higher-order Dirac semimetallic phase. In the phase, there appear Dirac points at both x-normal and y-normal surfaces. As a result, zero-energy hinge arcs arise in this semimetal. Apart from the zero-energy corner modes, the authors also find nonzero-energy corner modes in the energy spectrum as highlighted as red lines in figure 4(a). These modes are chiral with respect to φ, which can be seen as chiral hinge modes for the 3D system with φ acting as k_z . They arise from an effective Chern band localized on the x-normal surfaces, and each boundary band localized at one boundary contributes a Chern number of −1 [118]. They can also be characterized by the Chern number of the Wannier bands calculated from the lowest two degenerate occupied bands. This Wannier-sector Chern number can be interpreted as the spectral flow of the nested Wilson loop as in the case of 3D HOTI with chiral hinge states [38]. The topological invariant is equivalent to the Chern number of boundary bands.

These nonzero-energy corner modes can exist in the continuous bulk spectrum as highlighted by the white circle in figure 4(a). Such states are known as the bound states in the continuum [259–261]. The higher-order zero-energy bound states in the continuum have also been found in a 2D HOTI [262].

In this case, the zero-energy corner modes appear when $\phi \in (-\frac{3\pi}{4}, -\frac{\pi}{4})$ and $(\frac{\pi}{4}, \frac{3\pi}{4})$, which correspond to the regime with zero-energy edge modes in the 1D off-diagonal AAH model as shown in figure 4(c). Similarly, as in the previous case, the corner modes are characterized by the quadrupole moment (see figure 4(d)), which arises from the bulk energy gap closing, rather than the edge energy gap closing.

Similar to the previous case, the authors also find nonzero-energy chiral modes. However, different from the previous case, there appear a group of chiral modes in the bulk gaps which are boundary bands consisting of states localized at the 1D edges. These modes are characterized by the Chern number calculated in the $(k_y,\phi)$ or $(k_x,\phi)$ space for each k_x or k_y . In addition, there exist chiral modes (highlighted by red lines) connecting different boundary bands. These states are localized at the corners. Compared to the previous case, the boundary states do not have a well defined Chern number as these bands are chiral and not well separated. Viewing φ as k_z leads to a 3D model with Dirac points at zero energy in the bulk supporting the coexistence of zero-energy hinge modes, chiral hinge modes and chiral surface modes.

The irrational α_x and α_y lead to a 2D quasicrystal model where the nonzero-energy chiral modes localized at the edges and corners still exist [118]. Similar to the commensurate AAH models, the nonzero-energy corner modes can survive across the continuous bulk spectrum.

The Chern–Simons invariant is given by

$$\begin{equation} \theta = \frac{1}{4\pi} \int \mathrm{d}^3{\boldsymbol{k}} \, \epsilon_{abc} {\text{Tr}}\, \left[ \mathcal{A}_a\partial_b\mathcal{A}_c + \mathrm{i}\frac{2}{3} \mathcal{A}_a\mathcal{A}_b\mathcal{A}_c \right], \end{equation} \tag{ 81 }$$

where $\mathcal{A}_a$ is the Berry gauge field defined as $[\mathcal{A}_{a}]^{n,n^{^{\prime}}} = -\mathrm{i}\langle u_{\boldsymbol{k}}^n|\partial_{k_a}|u_{\boldsymbol{k}}^{n^{^{\prime}}} \rangle$ for $a = x,y,z$, where $|u_{\boldsymbol{k}}^n\rangle$ is the nth occupied eigenstate of the Bloch Hamiltonian. It has been proved that owing to the $C_4^z T$ symmetry, the Chern–Simons invariant is quantized to $\theta = 0,\pi \mod 2\pi$ [38].

The $\mathbb{Z}_2$ topology of 3D chiral HOTIss can be characterized by a Chern–Simons invariant θ, which has been used to characterize the 3D $\mathbb{Z}_2$ TIs protected by TRS [264]. When $\theta = \pi$, the system is in the higher-order topological insulating phase with chiral hinge modes. Moreover, in [62], the authors propose a generalized Pfaffian topological invariant to characterize the 3D HOTIs with $C_4^z T$ symmetry.

One can also use the nested Wilson loop to characterize the chiral higher-order topology [38]. The Wilson loop along x, $\mathcal{W}_{x, {\boldsymbol{k}}}$ is defined in equation (17) in section 2.1.4. For the first-order TIs with TRS, the Wannier spectrum of $\mathcal{W}_{x,(k_y, k_z)}$ is gapless, similar to the energy spectrum under open boundary conditions along x [265, 266]. For the chiral HOTI in equation (73), the Wannier spectrum is fully gapped corresponding to the gapped spectrum on the surfaces perpendicular to x, as shown in figure 6(a).

Zoom In Zoom Out Reset image size

It has been found that for the nontrivial phase of the chiral HOTI, the Wannier Hamiltonian has a Chern number $C = \pm 1$ for the half-filled Wannier bands, which can be seen as a 2D Chern insulator. This can be verified by computing the y-directed nested Wilson loop for a gapped Wannier sector of the Wannier Hamiltonian ${H}_{\mathcal{W}_x}(k_y, k_z)$. For the nontrivial phase, the nested Wilson loop spectrum displays a spectral flow with respect to k_z , indicating the existence of the Chern number for the Wannier bands, as shown in figure 6(b). This nontrivial winding can also be seen in the y-directed Wilson loop spectrum of a slab Hamiltonian with open boundaries along x (see figure 6(c)), illustrating the presence of chiral edge modes in the Wannier spectrum.

It has been proposed that the winding number of the quadrupole moment can serve as a topological invariant to characterize the 3D chiral topology [136, 137, 246]. The quadrupole moment winding with respect to k_z is defined as

$$\begin{equation} W_Q = \int_0^{2 \pi} \textrm{d} k_z \frac{\partial q_{x y}\left(k_z\right)}{\partial k_z}, \end{equation} \tag{ 82 }$$

where $q_{x y}(k_z)$ is the quadrupole moment defined in equation (13) for a 2D Hamiltonian at k_z . For the model in equation (73), $W_Q = 1$ in the nontrivial phase and $W_Q = 0$ in the trivial phase. This invariant can also be applied to many-body or disordered systems if we use a twisted boundary condition with the flux serving as the momentum [137]. Note that this invariant does not need a symmetry to protect its quantization. It can change through either a bulk or surface gap closing. Since this topological invariant is $\mathbb{Z}$-valued, it is possible to observe the phase transition from a trivial phase to a nontrivial phase with $W_Q = 2$.

The authors in [134, 135] prove that chiral symmetry can protect the quantization of the quadrupole moment (also see section 2.1.3 for proof) and find that disorder can drive a topological phase transition with bulk energy gap closing or edge energy gap closing from a trivial insulator to a HOTI, namely, a higher-order topological Anderson insulator. To introduce the disorder effects, we consider the following tight-binding Hamiltonian proposed in [134]

$$\begin{equation} \hat{H} = \sum_\textbf{r} \left[ \hat{c}^\dagger_\textbf{r}h_0\hat{c}_\textbf{r} +\left( \hat{c}^\dagger_\textbf{r}h_x\hat{c}_{\textbf{r}+\textbf{e}_x} +\hat{c}^\dagger_\textbf{r}h_y\hat{c}_{\textbf{r}+\textbf{e}_y}+\mathrm{H.c.} \right) \right]. \end{equation} \tag{ 96 }$$

Here $\hat{c}^\dagger_\textbf{r} = \left( \begin{array}{cccc} \hat{c}^\dagger_{\textbf{r}1} & \hat{c}^\dagger_{\textbf{r}2} & \hat{c}^\dagger_{\textbf{r}3} & \hat{c}^\dagger_{\textbf{r}4} \\ \end{array} \right)$ where $\hat{c}^\dagger_\mathbf{r \nu}$ ($\hat{c}_\mathbf{r \nu}$) are creation (annihilation) operator of a particle at the νth degree of freedom in a unit cell labelled by the position vector $\textbf{r} = (x,y)$ (x and y are integers). $\textbf{e}_x = (1,0)$ and $\textbf{e}_y = (0,1)$ are unit vectors along x and y, respectively.

$$\begin{equation} h_0 = \left( \begin{array}{cccc} 0 & -im_{\textbf{r}}^y & -im_{\textbf{r}}^x & 0 \\ im_{\textbf{r}}^y & 0 & 0 & i\bar{m}_{\textbf{r}}^x \\ im_{\textbf{r}}^x & 0 & 0 & -i\bar{m}_{\textbf{r}}^y \\ 0 & -i\bar{m}_{\textbf{r}}^x & i\bar{m}_{\textbf{r}}^y & 0 \\ \end{array} \right) \end{equation} \tag{ 97 }$$

describes the intra-cell hopping, and

$$\begin{align} \begin{aligned} &h_x = \left( \begin{array}{cccc} 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ t_{\textbf{r}}^x & 0 & 0 & 0 \\ 0 & -\bar{t}_{\textbf{r}}^x & 0 & 0 \\ \end{array} \right)~\text{and}~ h_y = \left( \begin{array}{cccc} 0 & 0 & 0 & 0 \\ t_{\textbf{r}}^y & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & \bar{t}_{\textbf{r}}^y & 0 \\ \end{array} \right) \nonumber\\ \end{aligned} \end{align} \tag{ 98 }$$

describe the inter-cell hopping along x and y, respectively [see figure 8(a) for the hopping parameters]. The inter-cell hopping magnitude are taken to be $t_\textbf{r}^x = \bar{t}_\textbf{r}^x = t_\textbf{r}^y = \bar{t}_\textbf{r}^y = 1$, and the other parameters $m_\textbf{r}^x$, $\bar{m}_\textbf{r}^x$, $m_\textbf{r}^y$, $\bar{m}_\textbf{r}^y$ are real numbers. The system preserves chiral symmetry (sublattice symmetry) given that only the nearest-neighbor hopping is involved in the model. Since the Hamiltonian contains the hopping that takes complex values, the spinless TRS represented by κ does not exist. However, the Hamiltonian can be mapped to the BBH model through a local unitary transformation [134] and thus share similar topological and localization properties with a disordered BBH model studied in [135].

Zoom In Zoom Out Reset image size

The disorder is introduced in the intra-cell hopping, that is,

$$\begin{align} m_\textbf{r}^\nu & = m_\nu+ W V_\textbf{r}^\nu, \end{align} \tag{ 99 }$$

$$\begin{align} \bar{m}_\textbf{r}^\nu & = {m}_\nu+ W \bar{V}_\textbf{r}^\nu \end{align} \tag{ 100 }$$

with $\nu = x,y$, where $m_x = m_y = m_0$ and ${V}_\textbf{r}^\nu$ and $\bar{V}_\textbf{r}^\nu$ obey independent uniform random distribution within the region of $[-0.5,0.5]$. Here, W represents the disorder strength. Due to the random property, the system is characterized by properties averaged over large numbers of random configurations.

In the clean case with $m_\textbf{r}^x = \bar{m}_\textbf{r}^x = m_x$ and $m_\textbf{r}^y = \bar{m}_\textbf{r}^y = m_y$, the system can be described by a Bloch Hamiltonian in momentum space

$$\begin{equation} H_0\left(\textbf{k}\right) = H_x\left(k_x,m_x\right)\otimes \sigma_3+\sigma_0 \otimes H_y\left(k_y,m_y\right), \end{equation} \tag{ 101 }$$

where $H_\nu(k_\nu,m_\nu) = \cos k_\nu \sigma_1+(m_\nu+\sin k_\nu)\sigma_2$ with $\nu = x,y$. This Hamiltonian takes the similar form to the Hamiltonian in equation (45), and thus can support zero-energy corner modes under open boundary conditions. When $m_x = m_y$, it has the C₄ symmetry so that the quadrupole insulator is intrinsic. Similar to the BBH model, the system is topologically nontrivial (trivial) when $|m_{0}|\lt1$ ($|m_0|\gt1$).

The higher-order topology of a QTI can be characterized by its quantized quadrupole moment. In the presence of the chiral symmetry, the system still supports quantized quadrupole moment as a topological invariant in the presence of disorder (see section 2.1.3).

In addition, the authors evaluated the quantized polarization p_x (p_y ) of effective boundary Hamiltonians (see the following introduction) for the y-normal (x-normal) boundary at half filling [134]. For disordered systems without translational symmetries, the polarization can be calculated according to the Resta formula [254, 282]. In the clean case of the model in equation (96), the derived effective boundary Hamiltonian at the y-normal (x-normal) edges is equal to $H_x(k_x,m_x)$ [$H_y(k_y,m_y)$] multiplyed by a nonzero factor as shown in [134]. Evidently, the HOTP arises when these effective boundary Hamiltonians become topological with nontrivial polarizations $p_x = p_y = 1/2$ which can be evaluated by the Resta formula in real space. Thus, the higher-order topology is characterized by the topological invariant

$$\begin{equation} P = 4|p_x p_y|. \end{equation} \tag{ 102 }$$

When P = 1, the system is in a HOTP, and when P = 0, it is in a trivial phase. Compared with the calculated quadruple moment in a system of $80\times 80$, P can be numerically evaluated in a system up to 500.

We now introduce the effective boundary Hamiltonian proposed in [283] and developed in [134] for disordered systems.

In general, we can write the real-space Hamiltonian in a quasi-1D tight-binding form with only nearest-neighor couplings along one direction when the unit cell is chosen to be large enough,

$$\begin{equation} H = \left(\begin{array}{ccccc} h_{1} & V_{1}^{\dagger}\\ V_{1} & h_{2} & V_{2}^{\dagger}\\ & V_{2} & \ddots & \ddots\\ & & \ddots & h_{N-1} & V_{N-1}^{\dagger}\\ & & & V_{N-1} & h_{N} \end{array}\right). \end{equation} \tag{ 103 }$$

Here, h_n is the Hamiltonian for the nth unit cell and V_n describes the hopping between the nth and $(n+1)$th unit cells for $n = 1,\cdots,N$. In disordered systems, the parameters in h_n and V_n can take random values. The full Green's function of the system is defined as $\mathbf{G}(\omega) = (\omega \mathbb{I}-H)^{-1}$. We will focus on the boundary Green's functions defined as the boundary blocks of the full Green's function matrix, $G_{1} = \mathbf{G}_{11}$ and $G_{N} = \mathbf{G}_{NN}$, which can characterize topologically nontrivial gapped boundaries of HOTIs.

The boundary Green's function can be computed recursively through the Dyson equation [283]

$$\begin{equation} \left(g_{N}^{-1} - V_{N-1}G_{N-1}V_{N-1}^{\dagger} \right) G_{N} = \mathbb{I}. \end{equation} \tag{ 104 }$$

Here, G_N ($G_{N-1}$) is the boundary Green's function of a system with N (N − 1) unit cells, and $g_{N}^{-1}(E) = \omega\mathbb{I}- h_{N}$ refers to the bare Green's function of the Nth unit cell. For a translational invariant system, it is expected that when N is sufficiently large, a fixed point G will be approached for the boundary Green's function so that it satisfies the following closed equation [283],

$$\begin{equation} \left(g^{-1} - VGV^{\dagger} \right) G = \mathbb{I}. \end{equation} \tag{ 105 }$$

In numerical computations, the boundary Green's function can be evaluated efficiently using the transfer matrix method. The boundary Green's function encodes information about boundary states of TIs [283]. We can define an effective boundary Hamiltonian from the boundary Green's function,

$$\begin{equation} H_{\mathrm{eff}} = -\left[G\left(E = 0\right)\right]^{-1}. \end{equation} \tag{ 106 }$$

For a QTI, the effective boundary Hamiltonians for both the x-edges and the y-edges have nontrivial topology. Note that if the bulk Hamiltonian respects the chiral symmetry, the effective boundary Hamiltonian also preserves the chiral symmetry so that we can use the 1D winding number or quantized polarization to characterize the edge topology [134]. This approach has also been applied to 3D third-order TIs with disorder [146].

However, in disordered systems unlike the translational invariant systems in [283], there is no fixed-point boundary Green's function in equation (105). Instead, there are fluctuations due to randomness for the effective boundary Hamiltonian $H_{\mathrm{eff}}$ over iterations (104). The authors in [134] use the statistical average over a large number of iteration steps to characterize the topology in disordered systems. In each iteration, the intra-cell hopping parts in h_n and V_n are randomly generated. Specifically, the polarization average along x of the effective boundary Hamiltonian at the y-normal boundary (similarly for x-normal one) is calculated,

$$\begin{equation} p_x = \frac{1}{N_i}\sum_{n = 1}^{N_i} |p_{x,n}|, \end{equation} \tag{ 107 }$$

where

$$\begin{equation} p_{x,n} = \frac{1}{2\pi}\text{Im}\text{log}\langle \Psi_n|e^{2\pi i\hat{x}/L_x}|\Psi_n\rangle \end{equation} \tag{ 108 }$$

[254] (note that the atomic positive charge contribution should be deducted), $\hat{x} = \sum_{x}x\hat{n}_x$ with $\hat{n}_x$ denoting the particle number operator at the site x, and L_x is the system's length along x. The polarization is evaluated for the many-body ground state $|\Psi_n\rangle$ of the boundary Hamiltonian $H_{n} = -G_{2n}(E = 0)^{-1}$ at half filling. Here, $G_{2n}$ is the 2nth boundary Green's function calculated via the Dyson equation. Note that one only considers the Green's functions at even steps since two distinct layers exist in the clean limit. As H_n also preserves chiral symmetry, $p_{x,n}$ is still quantized to either 0 or 0.5 for each iteration.

To see the existence of higher-order TAI phase, the authors consider a topologically trivial phase with $m_x = m_y\gt1$ in the clean limit. The phase diagram with respect to the disorder strength is mapped out in figure 8(b), showing remarkably the disorder-induced phase transition from a trivial insulator to a QTI. As shown in figure 8(b), the topological invariant P experiences a sudden change from 0 to 1 near W ≈ 2.1, implying that a topological phase transition happens there. The induced quadrupole insulator also manifests in the presence of corner states at zero energy. When the disorder strength is further increased, the system becomes gapless while still has quantized topological invariant P, which is named as the gapless HOTAI. Then, the system will enter a regime with nonquantized values of P due to the coexistence of trivial and nontrivial random samples. Finally, when the disorder is sufficiently strong, the topological invariant P vanishes and the system becomes a trivial Anderson insulator.

The results of the quadrupole moment and P (see figure 8(b)) show qualitative agreement with each other. Yet, there exists apparent discrepancy, which arises from finite-size effects. Though the quadrupole moment is quantized by chiral symmetry for a single random configuration, the averaged quadrupole moment may not be quantized to 0 or $1/2$ due to fluctuations. A finite-size scaling performed in [134] indicates that q_xy will approach $1/2$ in the thermodynamic limit.

The disorder-induced topological phase transition can be understood based on the renormalization of a Hamiltonian due to disorder [273, 284]. When disorder is weak, the phase boundaries modified by disorder can be evaluated by the SCBA method. The method has been applied to study disorder-induced quantum spin Hall insulator from a metallic trivial phase [273], disorder-induced phase transition between 3D quantum anomalous Hall state and a Weyl semimetal [277, 278], and disordered 3D TIs [284–286]. It has also been employed to explain the higher-order TAIs and identify their phase boundaries [134, 135]. Figure 8(d) shows the phase boundaries of the disordered model in equation (96) obtained by the SCBA method, which agree well with the numerical phase boundary determined by the topological invariant.

We now introduce the effective medium theory with SCBA [273, 284, 287]. For a Hamiltonian $H = H_0+V(\textbf{r})$ consisting of a clean system H₀ with translational symmetry and a disorder perturbation $V(\textbf{r})$, the disorder effect at the energy E can be described by the self-energy Σ defined by

$$\begin{equation} \left(E-H_{0}-\Sigma \right)^{-1} = \langle\left(E-H\right)^{-1}\rangle, \end{equation} \tag{ 109 }$$

where $\langle \cdot \rangle$ denotes the disorder ensemble average. Upon the disorder average, the system can still be described by a translational invariant Hamiltonian, where the self-energy renormalizes parameters or create new terms. To evaluate the self-energy, we use the SCBA which only takes into account contributions from rainbow diagrams with no crossed lines as shown in figure 9. The approximation neglecting the diagrams with crossed impurity lines usually works well in the regime of weak disorder. In general, the self-energy Σ as the summation of rainbow diagrams can be expressed in a recursive equation like

$$\begin{equation} \Sigma = W^2 U G\left(E\right) U, \end{equation} \tag{ 110 }$$

where $G = \left(E-H_{0}-\Sigma \right)^{-1}$ represents the full Green's function renormalized by disorder. G(E) is related to the bare Green's function $G_0(E) = (E-H_{0})^{-1}$ through Dyson series yielding $G = G_0+G_0 \Sigma G$ (see double solid lines in figure 9). In the recursive equation, W represents a factor proportional to the disorder strength which is contributed from the ensemble average for a disorder term in the random Hamiltonian (see dotted lines in figure 9). U represents the scattering matrix over the internal degrees of freedom due to the disorder term (see solid vertex in figure 9). The self-consistent equation (110) for a single type of disorder term can be easily generalized to cases for multiple disorder terms in the Hamiltonian.

Zoom In Zoom Out Reset image size

We remark that the effective Hamiltonian obtained from the SCBA is a result of the ensemble average of the Green's function. It thus does not correspond to the Hamiltonian for any single disorder configuration. Nevertheless, in the deep region of a topological or trivial phase at weak disorder, it is expected that the majority of random configurations have the similar topological properties except for a small portion of samples due to random fluctuations. Therefore, in these cases, the disorder averaged effective Hamiltonian can correctly describe the bulk topology of most random configurations as well as the presence or absence of their boundary states.

Specifically, we consider the disordered model in equation (96), where the random intra-cell hopping term at each unit cell r reads

$$\begin{align} V\left({\boldsymbol{r}}\right)& = W\left[ \begin{array}{cccc} 0 & -iV^y\left({\boldsymbol{r}}\right) & -iV^x\left({\boldsymbol{r}}\right) & 0 \\ iV^y\left({\boldsymbol{r}}\right) & 0 & 0 & i\bar{V}^x\left({\boldsymbol{r}}\right) \\ iV^x\left({\boldsymbol{r}}\right) & 0 & 0 & -i\bar{V}^y\left({\boldsymbol{r}}\right) \\ 0 & -i\bar{V}^x\left({\boldsymbol{r}}\right) & i\bar{V}^y\left({\boldsymbol{r}}\right) & 0 \\ \end{array} \right] \nonumber\\ & = W\left[\sum_{i = 1}^{4} V_{i}^{^{\prime}}\left({\boldsymbol{r}}\right) U_{i}\right]. \end{align} \tag{ 111 }$$

Here, $V_1^{^{\prime}} = (V^x+\bar{V}^x)/2$, $V_2^{^{\prime}} = (V^x-\bar{V}^x)/2$, $V_3^{^{\prime}} = (V^y+\bar{V}^y)/2$, $V_4^{^{\prime}} = (V^y-\bar{V}^y)/2$, and $U_1 = \sigma_y\otimes\sigma_z$, $U_2 = \sigma_y\otimes\sigma_0$, $U_3 = \sigma_0\otimes\sigma_y$, $U_4 = \sigma_z\otimes\sigma_y$. For the independent random distributions considered, we require

$$\begin{eqnarray} \langle V_i^{^{\prime}}\left({\boldsymbol{r}}\right) \rangle& = &0 \end{eqnarray} \tag{ 112 }$$

$$\begin{eqnarray} \langle V_i^{^{\prime}}\left({\boldsymbol{r}}_1\right) V_j^{^{\prime}}\left({\boldsymbol{r}}_2\right) \rangle & = &\frac{1}{24} \delta_{ij}\delta_{{\boldsymbol{r}}_1{\boldsymbol{r}}_2} \end{eqnarray} \tag{ 113 }$$

for $i,j = 1,2,3,4$, where $\langle \cdots \rangle$ represents the average over disorder configurations.

One thus can write the effective Hamiltonian at E = 0 as

$$\begin{equation} H_{\textrm{eff}}\left({\boldsymbol{k}}\right) = H_0\left({\boldsymbol{k}}\right)+\Sigma\left(E = 0\right). \end{equation} \tag{ 114 }$$

Here, with disorder, one can evaluate the self-energy Σ based on the following self-consistent equation [134]

$$\begin{equation} \Sigma\left(E\right) = \frac{W^2}{96\pi^2}\int_{\textrm{BZ} } \textrm{d}^2{\boldsymbol{k}} \sum_{n = 1}^{4} U_n G\left(E\right) U_n, \end{equation} \tag{ 115 }$$

where the renormalized Green's function $G(E) = [(E+i0^+)I-H_0({\boldsymbol{k}})-\Sigma(E)]^{-1}$. The integral is over the first Brillouin zone. This equation is equivalent to the general equation (110) represented in momentum space.

The self-energy is independent of momentum and only renormalize the onsite terms. At energy E = 0, numerical results suggest the existence of only three terms in the self energy so that it can be written as [134]

$$\begin{equation} \Sigma = i\Sigma_0 I + \Sigma_x \sigma_y\otimes\sigma_z + \Sigma_y \sigma_0\otimes\sigma_y, \end{equation} \tag{ 116 }$$

where $\Sigma_0,\Sigma_x,\Sigma_y$ take real values. These new terms arising from disorder thus renormalize the mass terms m_x and m_y as

$$\begin{align} m_x^{\prime} & = m_x + \Sigma_x, \end{align} \tag{ 117 }$$

$$\begin{align} m_y^{\prime} & = m_y + \Sigma_y. \end{align} \tag{ 118 }$$

$\Sigma_x$ and $\Sigma_y$ thus cause topological phase transitions, since the phase boundaries are determined by $m_x^{\prime} = 1$ and $m_y^{\prime} = 1$.

When disorder is weak, one can approximate the self-energy Σ by taking $\Sigma = 0$ on the right-hand side of equation (115),

$$\begin{equation} \Sigma_\nu = -\frac{W^2}{48\pi^2}\iint_{\textrm{BZ} } \textrm{d}{\boldsymbol{k}} \frac{m_\nu+\sin\left(k_\nu\right)}{F\left({\boldsymbol{k}}\right)}, \end{equation} \tag{ 119 }$$

where

$$\begin{equation} F\left({\boldsymbol{k}}\right) = 2+\sum_{\nu = x,y}\left[m_\nu^2+2m_\nu\sin\left(k_\nu\right)\right] \end{equation} \tag{ 120 }$$

with $\nu = x,y$. For $m_x\gt1$ and $m_y\gt1$, the integrands are positive so that both $\Sigma_x$ and $\Sigma_y$ are negative. The consequence is that for a sufficiently large disorder W, the renormalized parameters will enter a topological nontrivial phase region with $m_x^{\prime}\lt1$ and $m_y^{\prime}\lt1$ (see figure 8(d)).

We now review some concepts and methods used in the study of localization transitions in disordered systems.

Disorder plays an important role in the behavior of materials, impacting not only their topology, as previously discussed, but also their localization properties [269]. Sufficiently strong disorder can result in the absence of diffusion in a 3D lattice with random onsite energies and short-range coupling, as demonstrated by Anderson's pioneering work [268]. This phenomenon, known as Anderson localization, may arise from the mismatch between the random potential and the hopping energy between electrons at neighboring lattice sites [268].

A significant contribution to the understanding of localization is the one-parameter scaling theory introduced by Abrahams et al [288]. The key concept of this theory is the scaling function defined as

$$\begin{equation} \beta\left(\ln g\right) = \mathrm{d} \ln g/ \mathrm{d} \ln L, \end{equation} \tag{ 121 }$$

where g is the dimensionless conductance and L denotes the system size. In this context, $g \gg 1$ and $g \ll$ 1 correspond to extended and localized states, respectively. By examining the asymptotic behavior of β for large and small g, one can establish an approximate expression of β that effectively covers all values of g. In 3D systems, it is anticipated that a critical point g_c exists, separating regimes where β is positive and negative for $g \gt g_c$ and $g \lt g_c$, respectively. When the system size L is increased, the dimensionless conductance g tends to flow towards g = 0 ($g = \infty$) if it initially resides in the regime $g \lt g_c$ ($g \gt g_c$). As g is inversely correlated to the strength of disorder, the result implies the existence of a localization transition at finite disorder strength. In 1D and 2Ds, however, it is conjectured that β is always negative, indicating that Anderson localization occurs for arbitrarily weak disorder.

Despite the broad applicability of the scaling theory, exceptional cases exist in one and two-dimensional systems where states can be delocalized in the presence of disorder. Many of these cases are closely related to the topology of the system. Notable examples include the boundaries of two or three-dimensional topological insulators [16, 17, 289]. In such cases, metallic edge or surface states persist, protected by bulk topology, and remain immune to disorder as long as the bulk gap does not close. Another example arises in 2D systems precisely tuned to the critical point between two integer quantum Hall plateaus. At this transition point, the states become critical and display universal conductance and multifractal behavior, even in the presence of disorder [269].

In fact, the nontrivial bulk topology of a TI is intimately connected to the delocalization of the boundary states against disorder respecting symmetry. For example, the gapless surface states of a 3D TI or the chiral edge states of a Chern insulator are robust to scattering by random impurities. The analysis of Anderson localization at the boundary provides an alternative approach to obtain the Altland-Zirnbauer classification for TIs and superconductors [16, 17]. In brief, to study the Anderson localization problem, we consider noninteracting fermions in a random potential which can be described by an action of the nonlinear sigma model (NLσM) using replica trick in field theory. Then, the boundary states can evade Anderson localization when a topological term is allowed in the action of NLσM at the boundary, which is determined by the homotopy group of the target space in NLσM and directly related to the classification of bulk topology [16, 17].

Such intriguing localization properties can also manifest in disordered HOTIs. However, before delving into these examples, it is necessary to introduce several quantities used to characterize the localization properties of a system. These include the localization length, level-spacing ratio (LSR), inverse participation ratio (IPR), and fractal dimension.

The localization length serves as a measure of the extent of localization for the eigenstates. In the case of Anderson localization, the eigenstates are exponential localized as

$$\begin{equation} \psi\left(r\right) \sim e^{-r/\xi}, \end{equation} \tag{ 122 }$$

where ξ represents the localization length. The localization length can be computed by partitioning the Hamiltonian into layers and employing the formula [290]

$$\begin{equation} 2/\xi = - \lim_{N\rightarrow \infty} N^{-1}\ln \mathrm{Tr} \left(G_{1N} G_{1N}^\dagger\right). \end{equation} \tag{ 123 }$$

Here, $G_{1N}$ is a submatrix of the Green's function $\mathbf{G}(\omega) = (\omega \mathbb{I} - H)^{-1}$ for a Hamiltonian H with N layers (equation (103)), which can be calculated by the iteration $G_{1N} = G_{1, N-1} V_{N-1}^\dagger G_{N}$ and $G_{N} = (\omega \mathbb{I} - h_{N} - V_{N-1} G_{N-1} V_{N-1}^\dagger)^{-1}$. For two (or three) dimensions, the normalized localization length $\Lambda = \xi/L$ is commonly used, with L (or L²) denoting the size of the cross-section. For extended (localized) states, Λ increases (decreases) with L, while critical states exhibit a constant value of Λ with respect to L.

One can also use the LSR to identify the localization properties defined as [291]

$$\begin{equation} r\left(E\right) = \frac{1}{N_E - 2} \sum_{i = 1}^{N_E - 2} \mathrm{min}\left(\delta_i, \delta_{i+1}\right)/\mathrm{max}\left(\delta_i, \delta_{i+1}\right), \end{equation} \tag{ 124 }$$

where $\delta_i = E_{i+1} - E_{i}$ represents the level spacing between adjacent energy levels and N_E denotes the number of selected energy levels around energy E. For a random Hamiltonian belonging to the Gaussian orthogonal ensemble, the LSR for extended states is given by r ≈ 0.53 [291, 292]. In the case of strong disorder where the states become localized, the LSR is approximately r ≈ 0.386, which arises from the Poisson statistics of the level spacing [291, 292].

Another quantity to examine the localization property is the IPR defined by

$$\begin{equation} I\left(E\right) = \frac{1}{N_E} \sum_{i = 1}^{N_E} \sum_{\boldsymbol{r}} \left( \sum_\nu |\Psi_{E_i, {\boldsymbol{r}}\nu}|^2 \right)^2, \end{equation} \tag{ 125 }$$

where ν represents the internal degrees of freedom. The IPR can be used to determine whether the states around energy E are localized or extended. For an extended state, the IPR follows the scaling $I \propto 1/L^d$ where d is the dimension of the system. For example, a state that is evenly distributed in the full lattice has an IPR of $I = 1/L^d$. While for a localized state, the IPR tends to be a constant as the system size increases since the distribution of the state remains unchanged. In the limiting case where the state is localized at a single site, the IPR is I = 1. Based on the IPR, one can further define the fractal dimension D₂ such that $I \propto 1/L^{D_2}$, with D₂ = 0 and $D_2 = d$ corresponding to localized and extended states, respectively. When $0 \lt D_2 \lt d$, the states are neither exponentially localized nor fully extended, but exhibit multifractal behaviors [293].

When $m_x = m_y$ in the model in equation (96), the disorder-induced higher-order topological phase transition occurs through the bulk energy gap closure, as shown in figure 8(c). In the ensemble of disordered systems, there in fact emerges an average C₄ symmetry, that is, a Hamiltonian $\hat{H}$ for a certain disorder configuration occurs with the same probability as its symmetry conjugate partner $\hat{C}_4\hat{H} \hat{C}_4^{-1}$. The average symmetry enforce the bulk energy gap closure (see section 4.1 for detailed discussion). When $m_x \neq m_y$, the symmetry is broken so that the topological phase transition occurs via the edge energy gap closing [135].

At the critical point, the normalized localization length $\Lambda_x = \lambda_x/L_y$ (similarly for $\lambda_y/L_x$) at zero energy becomes scale free, indicating the divergence localization length in the thermodynamic limit, as shown in figure 10(a). This is associated with the delocalized states occurring at the phase transition. On the other hand, at both sides of the critical point, the normalized localization length decreases as the transverse size increases, indicating the localization of states around zero energy.

Zoom In Zoom Out Reset image size

The localization properties can also be confirmed by looking at the energy level statistics from the LSR and the real-space distribution from the IPR (see section 3.1.2). The LSR of zero-energy states at the criticality is close to the value 0.53, corresponding to extended states (see figure 10(c)). The fractal dimension D₂ from the scaling of the IPR is close to the bulk dimension 2, indicating the states are extended in the bulk (see figure 10(d)). Therefore, the topological phase transition accompanies a localization-delocalization-localization (LDL) transition with the bulk gap closing.

For the model studied in [135], the disorder-driven phase transition between topologically trivial and nontrivial phase is accompanied by a LDL transition with gap closing at open boundaries while the bulk remains insulating, corresponding to the case with $m_x \neq m_y$ in the model in equation (96). Li et al performed the finite-size scaling of the localization length at zero energy for a quasi-one-dimensional ribbon with open boundaries in the transverse dimension. As shown in figure 11, around the transition point, the localization length increases monotonically with respect to the ribbon width, which signifies the divergence in the thermodynamic limit. The divergence of the localization length under open boundary conditions indicates the occurrence of zero-energy delocalized states at the critical point due to gap closing along the corresponding boundaries parallel to the longitudinal direction. In contrast, with a periodic boundary condition along transverse direction, the normalized localization length decreases as the width increases, which indicates the absence of delocalized bulk states at the phase transition.

Zoom In Zoom Out Reset image size

For sufficiently strong disorder, the system becomes gapless, leading to gapless higher-order TAI, trivial-II phase and the Griffiths regime (see figure 8). In the former two phases, the states near zero energy are localized, as revealed by the decrease of the localization length with the transverse size in figure 10(a). The localized nature of the states around zero energy can also be evidenced by their LSR approaching 0.386, and the relatively larger IPR than the regime with weak disorder (see figure 10(c)). More precisely, the localization behavior should be characterized by the fractal dimension D₂ obtained from the scaling of the IPR versus the system size, which provides more information than the magnitude of IPR itself. As shown in figure 10(d), the fractal dimension D₂ is around zero reflecting the states near zero energy is localized in the gapless higher-order topological regime. Further finite-size analysis indicates that D₂ approaches zero in the thermodynamic limit [134]. The authors further found that the states of higher energy are all localized via the LSR and localization length calculations, which indicates that the higher-order topology can be carried by the localized bulk states. For other models, the existence of the mobility edge in the gapless HOTP has also been found [143]

The Griffiths regime (or critical regime) has more interesting property that the zero-energy normalized localization length displays scale invariant behavior as the transverse size increases, suggesting a multifractal phase therein (see figures 10(a) and (b)). This multifractal regime appears as a critical region residing between two localized phases, which differs from the conventional case with a critical point between delocalized and localized phase. In the Griffiths regime, other states at nonzero energy remain localized (see figure 10(b)). The multifractal behavior at zero energy can be verified by the LSR close to 0.53 (see figure 10(c)) and the fractal dimension D₂ between 0 and 2 (see figure 10(d)). Similar localization properties have also been observed in the case with $m_{x} = m_{y}\lt1$.

The tight-binding Hamiltonian $\hat{H}_c$ studied in [137] includes the onsite mass term $T_0 = M\tau_z\sigma_0$ and the hopping matrix

$$\begin{align} T_c\left({\boldsymbol{d}}\right) = f\left(|{\boldsymbol{d}}|\right) \left[ t_0\tau_z\sigma_0+it_1\tau_x\left(\boldsymbol{\hat{d}\cdot\sigma}\right)+t_2\left(\hat{d}_x^2-\hat{d}_y^2\right)\tau_y \sigma_0 \right] \end{align} \tag{ 131 }$$

from site r to $\boldsymbol{r+d}$ with ${\boldsymbol{\hat{d}}} = {\boldsymbol{d}}/|{\boldsymbol{d}}| = (\hat{d}_x, \hat{d}_y,\hat{d}_z)$ and $f(|{\boldsymbol{d}}|) = \Theta (R-|{\boldsymbol{d}}|)\exp[-\lambda(|{\boldsymbol{d}}|/d_0-1)]/2$. For regular cubic lattices with only nearest-neighbor hoppings, the Hamiltonian reduces to the paradigmatic 3D SOTI model in equation (73). The SOTI phase is protected by the winding number of the quadrupole moment about the momentum k_z , which requires no spatial symmetry, and thus it is possible for SOTIs to persist in completely random lattice sites. In amorphous systems, there is no well-defined momentum k_z due to the lack of translational symmetry. Therefore, the authors insert a flux $\Phi_z$ through a twisted the boundary condition along z in place of k_z and calculate the winding number of the quadrupole moment with respect to $\Phi_z$ [137], defined as

$$\begin{align} W_Q = \int_{0}^{2\pi} \textrm{d}\Phi_z \frac{\partial Q_{xy}\left(\Phi_z\right)}{\partial \Phi_z}, \end{align} \tag{ 132 }$$

where $Q_{xy}(\Phi_z)$ is the quadrupole moment in the (x, y) plane, which is a function of the inserted flux $\Phi_z$. The flux can also be added by introducing a phase factor $\exp(i\Phi_zd_z/L_z)$ to the hopping term. In a SOTI phase, a nontrivial winding number of $Q_{xy}(\Phi_z)$ arises when one tunes $\Phi_z$ from 0 to 2π, as illustrated in figure 17(b) for a disorder configuration.

The study finds the existence of a second-order topological region with almost quantized winding number of one and quantized two-terminal conductance of $G = 2e^2/h$ along the z direction arising from chiral hinge states (see figure 17(c)).

In the presence of TRS, the authors construct a tight-binding Hamiltonian [137]

$$\begin{align} \hat{H}_h = \sum_{\boldsymbol{r}} \left[\hat{c}_{\boldsymbol{r}}^\dagger T_0\hat{c}_{\boldsymbol{r}}+\sum_{\boldsymbol{d}} \hat{c}_{{\boldsymbol{r}}+{\boldsymbol{d}}}^\dagger T_h\left({\boldsymbol{d}}\right)\hat{c}_{\boldsymbol{r}} \right], \end{align} \tag{ 133 }$$

where $\hat{c}_\textbf{r}^\dagger = (\hat{c}_\textbf{r,1}^\dagger,\cdots,\hat{c}_\textbf{r,8}^\dagger)$, the onsite mass term $T_0 = M\tau_zs_0\sigma_0$ and the hopping matrix

$$\begin{align} T_h\left({\boldsymbol{d}}\right) & = f\left(|{\boldsymbol{d}}|\right) \left[t_0\tau_zs_0\sigma_0+it_1\tau_xs_0\left(\boldsymbol{\hat{d}\cdot\sigma}\right) \right. \nonumber \\ & \quad \left.+t_2\left(\hat{d}_x^2-\hat{d}_y^2\right) \tau_ys_y \sigma_0+it_3\hat{d}_z\tau_ys_x\sigma_0 \right] \end{align} \tag{ 134 }$$

with $\{s_\nu\}$ denoting another set of Pauli matrices. The TRS represented by $\hat{T}$ is respected, i.e. $\hat{T}\hat{H}_h\hat{T}^{-1} = \hat{H}_h$, where $\hat{T}i\hat{T}^{-1} = -i$ and $\hat{T}\hat{c}_{\boldsymbol{r}}\hat{T}^{-1} = U_T\hat{c}_{\boldsymbol{r}}$ with $U_T = i\tau_0s_0\sigma_y$. On a cubic lattice, the Hamiltonian reduces to the model in equation (84) with an extra term proportional to $t_3\sin k_z \tau_y s_x$.

When t₃ = 0, the pseudospin s_y is a conservative quantity and the topological properties can be identified by the winding number of the spin quadrupole moment about a flux $\Phi_z$,

$$\begin{align} W_Q^{\left(s\right)} = \int_{0}^{2\pi} \textrm{d}\Phi_z \frac{\partial Q_{xy}^{\left(s\right)}\left(\Phi_z\right)}{\partial \Phi_z}, \end{align} \tag{ 135 }$$

where $Q_{xy}^{(s)}(\Phi_z)$ is the spin quadrupole moment as defined in section 2.4.2 but with the momentum k_z replaced by the flux $\Phi_z$. As a result, the HOTP can be characterized by the spin winding number defined in equation (86). It is also shown that for small t₃, the spin winding number is still applicable even though s_y is not a conservative quantity. Besides, a $\mathbb{Z}_2$ invariant is defined to characterize the topological phase [see equation (91) in section 2.4.2 for its definition. Here, the momentum k_z should be replaced by a flux $\Phi_z$].

Numerical calculations reveal a range of M with sample averaged $\mathbb{Z}_2$ topological invariant $\overline{\nu}_Q$ and spin winding number $\overline{W}_{QS}$ quantized at 1 up to finite-size effects, as shown in figure 18(b), indicating the existence of an amorphous SOTI phase with TRS. The topological helical hinge states manifest in the two-terminal conductance quantized at $4e^2/h$.

Zoom In Zoom Out Reset image size

Structural disorder refers to position randomness of lattice sites, serving as a bridge between a crystalline system and an amorphous system. The authors in [137] also demonstrate that structural disorder can drive a topological phase transition from a trivial phase to a SOTI. The induced transition can be seen from the jump of the conductance and the winding number of the quadrupole moment from zero and the bulk energy gap closure as shown in figure 18(a). The driven SOTI also exhibits the hinge states as revealed by the local DOS.

In figures 17(d) and 18(a), we see that the higher-order topological phase transition occurs both through the closure of the bulk energy gap, although the $\hat{C}_4^z\hat{T}$ symmetry is broken in each sample. Here, we will follow [137] to show that such a bulk energy gap closure happens due to the $\hat{C}_4^z\hat{T}$ symmetry on average. One can also apply similar argument to the 2D case in section 3.1.

To demonstrate that the Hamiltonian $\hat{H}_c$ respects the average $\hat{C}_4^z\hat{T}$ symmetry, one can calculate the symmetry counterpart of the Hamiltonian,

$$\begin{align} & \hat{C}_4^z \hat{T}\hat{H}_c\left(\hat{C}_4^z \hat{T}\right)^{-1}\nonumber\\ & \quad = \sum_{\textbf{r}\in S} \left[ M\hat{c}_{D_{\hat{C}_4}\textbf{r}}^\dagger \tau_z\sigma_0 \hat{c}_{D_{\hat{C}_4}\textbf{r}} +\sum_{\substack{\textbf{d} = \textbf{r}_1-\textbf{r}\\ \textbf{r}_1 \in S, \textbf{r}_1\neq \textbf{r} }} \hat{c}_{D_{\hat{C}_4}\left(\textbf{r}+\textbf{d}\right)}^\dagger T_c\left(D_{\hat{C}_4}\hat{\mathbf{d}}\right) \hat{c}_{D_{\hat{C}_4}\textbf{r}} \right] \end{align} \tag{ 136 }$$

$$\begin{align} & \quad = \sum_{\textbf{r}^{\prime} \in S^{\prime} } \left[ M\hat{c}_{\textbf{r}^{\prime}} \dagger \tau_z\sigma_0 \hat{c}_{\textbf{r}^{\prime}} +\sum_{\substack{\textbf{d}^{\prime} = \textbf{r}_1^{\prime}-\textbf{r}^{\prime}\\ \textbf{r}_1^{\prime} \in S^{\prime}, \textbf{r}_1^{\prime} \neq \textbf{r}^{\prime} }} \hat{c}_{\textbf{r}^{\prime}+\textbf{d}^{\prime}} \dagger T_c\left(\hat{\mathbf{d}}^{\prime}\right) \hat{c}_{\textbf{r}^{\prime}} \right]. \end{align} \tag{ 137 }$$

In the derivation, the results are used: $\hat{C}_4^z \hat{T} \hat{c}_\textbf{r}(\hat{C}_4^z \hat{T})^{-1} = i\tau_0\sigma_y e^{i\frac{\pi}{4}\sigma_z}\hat{c}_{D_{\hat{C}_4}\textbf{r}}$ and $\hat{C}_4^z \hat{T}i(\hat{C}_4^z \hat{T})^{-1} = -i$. We use S to represent a set of all site positions in a disorder sample and then rotate all site positions in S in a counterclockwise direction about z by 90 degrees to obtain a new set $S^{\prime} = D_{\hat{C}_4} S\equiv \{D_{\hat{C}_4}\textbf{r} = (-y,x,z): \textbf{r}\in S\}$. If we consider cubic lattices, then $S^{\prime} = S$ so that $\hat{C}_4^z \hat{T}\hat{H}_c(\hat{C}_4^z \hat{T})^{-1} = \hat{H}_c$, showing that the $\hat{C}_4^z \hat{T}$ symmetry is obeyed. However, for a typical individual disorder configuration, $S\neq S^{\prime}$, and we usually have $\hat{C}_4^z \hat{T}\hat{H}_c(\hat{C}_4^z \hat{T})^{-1} \neq \hat{H}_c$, indicating that the symmetry is absent.

We now consider a statistical ensemble consisting of all systems. The ensemble is defined to respect an average symmetry if a Hamiltonian and its symmetry conjugate partner exist in the ensemble with the same probability [279, 280]. For instance, a configuration S and its rotational part $S^{\prime} = D_{\hat{C}_4} S$ clearly appear with the same probability since we consider the structural disorder without any correlation. As a consequence, the probability of occurrence of $\hat{H}_c$ is equal to that of its $\hat{C}_4^z \hat{T}$ conjugate partner $\hat{C}_4\hat{T}\hat{H}_c(\hat{C}_4\hat{T})^{-1}$ so that the ensemble respects the $\hat{C}_4\hat{T}$ symmetry on average.

We now argue that this average symmetry enforces the bulk energy gap closure at the phase transition point for a single system [137]. To demonstrate it, we consider a subsystem $\hat{H}_1$ whose z coordinates are restricted in the interval $[z_0,z_0+\Delta z]$ inside a large system (see figure 19). As a single system with a spatial configuration S₁, the subsystem does not respect the $\hat{C}_4^z \hat{T}$ symmetry so that the topological phase transition happens through a surface energy gap closing, e.g. the energy gap on the x-normal surface. For a very large system (e.g. infinitely long along z), it is reasonable to expect that there exists the spatial configuration $S_2 = D_{\hat{C}_4} S_1+\alpha \textbf{e}_z$, which is obtained by rotating S₁ about z by 90 degrees and shifting the coordinate along z by α. This leads to another subsystem $\hat{H}_2$ with the spatial configuration S₂ (see figure 19). It is clear to see that $\hat{H}_2$ is the $\hat{C}_4^z \hat{T}$ conjugate partner of $\hat{H}_1$. Therefore, when the x-normal surface energy gap of $\hat{H}_1$ vanishes, the y-normal surface energy gap of $\hat{H}_2$ must close. It indicates that for the entire system including both $\hat{H}_1$ and $\hat{H}_2$, the closure of the energy gap must simultaneously occur on both surfaces, suggesting that the bulk energy gap should close.

Zoom In Zoom Out Reset image size