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

In recent years, higher-order topological phases have attracted great interest in various fields of physics. These phases have protected boundary states at lower-dimensional boundaries than the conventional first-order topological phases due to the higher-order bulk-boundary correspondence. In this review, we summarize current research progress on higher-order topological phases in both crystalline and non-crystalline systems. We firstly introduce prototypical models of higher-order topological phases in crystals and their topological characterizations. We then discuss effects of quenched disorder on higher-order topology and demonstrate disorder-induced higher-order topological insulators. We also review the theoretical studies on higher-order topological insulators in amorphous systems without any crystalline symmetry and higher-order topological phases in non-periodic lattices including quasicrystals, hyperbolic lattices, and fractals, which have no crystalline counterparts. We conclude the review by a summary of experimental realizations of higher-order topological phases and discussions on potential directions for future study.

Quantum phases of matter, according to Landau's paradigm, were believed to be classified based on the symmetry of the ground state for a long time.However, in 1980, the discovery of quantum Hall effect showcased the example of topological phases of matter, for which different phases cannot be distinguished by spontaneously symmetry breaking but are characterized by the bulk topology [1,2].In 1988, the famous Haldane arXiv:2309.03688v2[cond-mat.mes-hall]14 Sep 2023 model was proposed to describe a topological insulator (TI) with quantum Hall conductance resulted from chiral edge states, which is characterized by a nonzero Chern number of bulk wave functions [3].The connection between nontrivial bulk topology and topologically protected boundary states is referred to as bulk-boundary correspondence of TIs.Later, inspired by the discovery of time-reversal-invariant TIs [4][5][6][7], the concept of TIs has been extended to symmetry-protected topological states which are protected by either internal symmetries or crystalline symmetries [8][9][10][11][12][13][14][15].The symmetry-protected TIs cannot be adiabatically deformed into trivial insulators without closing the bulk gap if the symmetry is not broken, and they can be classified based on homotopy theory [16,17] and symmetry indicators developed for topological crystalline insulators [18][19][20].Besides TIs, people have also discovered topological superconductors described by gapped mean-field Hamiltonian [21][22][23] and topological semimetals with topological gapless points and protected surface states [24][25][26][27].During the past decades, the topological phases of matter have attracted broad research interest and become a very active field in condensed matter physics [28][29][30][31][32].
The conventional topological phases in d dimensions host (d − 1)-dimensional boundary states according to the bulk-boundary correspondence [9,10].In recent years, a new type of topological phases dubbed higherorder topological insulators (HOTIs) were discovered, for which a d-dimensional nth-order topological phase supports (d − n)-dimensional boundary states with the codimension n > 1 .In 2017, Benalcazar et al. generalized the concept of the dipole moment to the multipole moment and proposed a model for multipole TIs protected by mirror symmetries, now known as the Benalcazar-Bernevig-Hughes (BBH) model, which exhibits quantized multipole moment and in-gap corner modes as the defining signature of higher-order topology [33,34].Specifically, a quadrupole topological insulator is featured by quantized edge polarizations and fractional corner charges resulted from corner modes.The corner modes arise as the manifestation of nontrivial topology of edge Hamiltonians and cannot be removed by symmetry-preserving perturbations unless the bulk or edge energy gap closes [33,34,66,71].Analogously, an octupole topological insulator is a 3D third-order topological insulator with the quantized octupole moment and corner modes [33].In fact, it is not always true that a quadrupole TI always has corner modes.For instance, very recently, a quadrupole TI (with nonzero quadrupole moment) without corner modes in the energy spectrum has been found [76,77].
Soon after the discovery of the multipole TIs, several works discovered 3D HOTIs with gapped surface states but gapless hinge states, in contrast to first-order TIs with gapless surface states [35][36][37].In the absence of time-reversal symmetry (TRS), a typical class of secondorder topological insulators (SOTIs) is featured by chiral hinge states, which can be protected by the prod-uct of time-reversal and rotation symmetry such as the C 4 T symmetry [37].In the presence of TRS, there exist Kramers pairs of hinge states, namely, helical hinge states, which can be protected by additional crystalline symmetry such as mirror or rotation symmetry [35][36][37].These gapless hinge modes give rise to quantized conductance, which is the transport signature of the higherorder topology.For many HOTIs with crystalline symmetry, the anomalous boundary states can be understood as the domain walls between massive Dirac Hamiltonians of adjacent gapped surfaces with opposite signs of masses due to the symmetry constraints [49].In fact, in early studies [78,79], similar in-gap hinge modes were found in 3D TIs, where the surface states are gapped out by the magnetic field or magnetization breaking the TRS.The chiral hinge modes emerge at the interface of adjacent surfaces with opposite signs of magnetic gaps.However, according to higher-order bulk-boundary correspondence proposed in Ref. [49], these higher-order topological phases (HOTPs) depending on the surface orientations or lattice terminations belong to extrinsic HOTPs, in contrast to intrinsic HOTPs originating from the nontrivial bulk topology protected by crystalline symmetry.The concept of higher-order topology has also been extended to higher-order topological superconductors hosting Majorana modes at corners or hinges , and higher-order topological semimetals hosting gapless points and Fermi arcs at hinges [107][108][109][110][111][112][113][114][115][116].
Since the pioneering works reporting HOTIs, it has witnessed a rapid development of both theoretical and experimental research in this field, and many models with higher-order topology were studied .For example, in 2Ds, the SOTIs with corner modes were found in breathing kagome and pyrochlore lattices, which are characterized by the quantized bulk polarizations [39].In Ref. [57], the authors systematically investigated 2D SOTIs with quantized fractional corner charges protected by rotation symmetries.Moreover, several works studied fractional charge responses to topological defects such as dislocations and disclinations, which are related to higher-order bulk topology [117][118][119][120][121]. In 3Ds, it was found that HOTIs with anomalous hinge states can be protected by various crystalline symmetries including inversion symmetry, rotation symmetry, and magnetic symmetry such as C 2 T [43][44][45][46][47][48][49][50].In addition, the space-time inversion symmetry (spinless P T symmetry) can protect 2D fragile topology with corner modes and 3D nodal-line semimetals with hinge Fermi arcs [51,52,108,114].Very recently, a Klein-bottle HOTI with non-chiral hinge modes protected by a pair of momentum-space glide reflection symmetries has been found [122].
In the presence of crystalline symmetry, HOTIs can usually be detected by symmetry indicators [36,47,57].On the other hand, the higher-order topology can also be characterized by some bulk topological invariants.For instance, the 2D quadrupole TIs are characterized by the quadrupole moment as a Z 2 topological invari-ant [123,124].It was proved that the quadrupole moment is quantized by chiral symmetry [125,126] or particle-hole symmetry [126].The C 4 T symmetric SO-TIs with chiral hinge modes can be characterized by the Chern-Simons invariant [37].Later, it was found that these chiral hinge modes can be also characterized by the winding number of the quadrupole moment [127,128].For the SOTIs with TRS, the helical hinge states can be characterized by the mirror Chern number in the presence of mirror symmetries [37].In the absence of mirror symmetries, the helical modes can be characterized by a Z 2 invariant proposed in Ref. [128].Remarkably, either the quadrupole moment winding for chiral modes and the Z 2 invariant for helical modes do not require the presence of crystalline symmetries.Therefore, both the 2D quadrupole TIs and 3D chiral or helical SOTIs can exist in systems without any crystalline symmetry.
People have also studied the effects of disorder on the higher-order topology, and found that the higher-order topological states remain robust against weak disorder [33,69,94,125,126,[129][130][131][132][133][134][135][136][137][138].More remarkably, the disorder is not always detrimental to topological phases.In two works studying quadrupole TIs in the presence of quenched disorder, it was found that disorder can induce a topological phase transition from a trivial insulator to a quadrupole TI with nonzero quadrupole moment quantized by chiral symmetry [125,126].The HOTP induced by disorder is named higher-order topological Anderson insulator, following the topological Anderson insulator found in first-order topological phases [139].In an earlier study [94], the authors studied a second-order topological superconductor with Majorana corner modes and found that an initially trivial phase can enter into a HOTP upon adding disorder.The disorder-driven topological phase transitions can be explained by the renormalization of Hamiltonian parameters due to disorder.A subsequential study [140] considered a modified Haldane model with disorder in the hopping phases and found that the disorder can drive a transition from a Chern insulator to a HOTIs with corner states, which was further observed in experiment.
In Ref. [141], the authors showed that a HOTI with quantized quadrupole moment and zero-energy corner modes can exist in a 2D amorphous system without translation symmetries, for which the higher-order topology is actually protected by chiral symmetry.Later, the authors in Ref. [128] found SOTIs with either chiral or helical hinge states in a 3D amorphous systems, which are characterized by nontrivial topological invariants, quantized longitudinal conductance, and in-gap hinge states.More intriguingly, it was found that the structural disorder can induce a SOTI starting from a trivial insulator on a regular lattice.The topological transition undergoes a bulk gap closing guaranteed by the average C 4 T symmetry.Recently, it has been realized that the average symmetry plays an important role in the classification of symmetry-protected topological phases [142].
Apart from amorphous systems without any spa-tial symmetry, the HOTPs with in-gap corner modes have been found to exist in other non-crystalline systems such as quasicrystalline [143][144][145][146][147][148][149][150][151][152], hyperbolic lattices [153,154], and fractal lattices [155][156][157].Since quasicrystalline and hyperbolic lattices can possess rotation symmetries that do not exist in regular crystals, they can harbor topological phases without crystalline counterparts.For instance, it was found that a HOTP hosting eight corner modes exists in a quasicrystal with eightfold rotation symmetry, which can be characterized by a Z 2 invariant protected by particle-hole symmetry [143].Similar phases were also found in 2D amorphous lattices with average rotation symmetry [158] and hyperbolic lattices with various rotation symmetries [153,154].In Ref. [159], the authors studied a 3D quasicrystalline lattice constructed by stacking 2D quasicrystals and found the existence of a 3D SOTIs with eight hinge states or helical pairs of hinge states protected by the winding number of the generalized quadrupole moment and a Z 2 invariant based on transformed lattice sites, respectively.In fact, it is found that the 2D HOTI on a quasicrystalline lattice can also be characterized by the generalized quadrupole moment [159].
There have been some reviews on the topics of HOTPs in crystalline systems.In Ref. [49,160], the authors introduced the higher-order bulk-boundary correspondence to classify the bulk and boundary topology of HOTPs.Another excellent review [161] introduced celebrated models of HOTPs and mainly focused on experimental realizations in various classical platforms, such as metamaterials and electric circuits.
The review is organized as follows.In Sec.II, we review several prototypical models of HOTPs in crystals and introduce their topological properties.In Sec.III, we review the studies on disorder effects on HOTPs and demonstrate disorder-induced HOTPs.In Sec.IV, we review the studies on HOTPs in non-crystalline systems including amorphous lattices, quasicrystalline lattices and hyperbolic lattices.In Sec.V, we briefly summarize the experimental progress on HOTPs.Finally, in Sec.VI, we give a summary and perspective.
Here we will introduce a class of HOTPs named quantized electric multipole insulators or multipole TIs.They are characterized by quantized multipole moment and zero-dimensional corner states, which are generalizations of 1D TIs with quantized polarization.We mainly focus on a prototypical model of quadrupole TIs known as the BBH model [see Fig. 1(a)] proposed in Ref. [33], and discuss its bulk and boundary properties and corresponding topological invariants.The model for octupole TIs can be constructed analogously [see Fig. 1(b)], which belongs to third-order topological phases [33].

Benalcazar-Bernevig-Hughes model
The BBH model for quadrupole insulators is a fourband model with Bloch Hamiltonian [33] where γ x,y and λ x,y are real parameters and with U Mx = τ 1 σ 3 and U My = τ 1 σ 1 .When γ x = γ y and λ x = λ y , the Hamiltonian (1) respects a four-fold rotation symmetry C 4 with (C 4 ) 4 = −1, where Hamiltonian (1) also has spinless TRS T represented by the complex conjugation κ, particle-hole symmetry P represented by τ 3 σ 0 κ, and chiral symmetry C represented by τ 3 σ 0 [33].We now prove that the TRS and the inversion symmetry U I = σ 2 , which anticommute with each other, guarantee that each energy band is at least twofold degenerate.Let |u n k ⟩ be an eigenstate of the Hamiltonian (1) with the energy E n k at k. Then U I T |u n k ⟩ corresponds to an eigenstate with the same energy at k. Assume that there were no degeneracy at k, which is equivalent to saying that the two states refer to the same one up to a phase θ, i.e.
On the other hand, due to the anticommutation relation of U I and T , ( , which is impossible.We therefore conclude that each energy band is at least two-fold degenerate. Due to the π-flux in each plaquette, the mirror symmetries M x and M y anticommute with each other, namely, {M x , M y } = 0, which leads to the gapped Wannier bands [33].The mirror symmetries or fourfold rotation symmetry can enforce the quantization of the bulk quadrupole moment q xy (see definition in Sec.II A 3).
Owing to the mirror symmetries, the BBH model has zero bulk polarization and zero Chern number so that there are no first-order edge states [33].The absence of the Chern number can also be seen from the existence of the TRS.However, when |γ x /λ x | < 1 and |γ y /λ y | < 1, the model supports zero-energy corner modes localized at four corners in the energy spectrum with full open boundary conditions [see Fig. 1(c)].For clarity, we consider a limiting case with γ x = γ y = 0 to illustrate these nontrivial signatures for the quadrupole topological phase illustrated in Fig. 1(d).It is clear to see the existence of four isolated sites, contributing to four zero-energy corner modes.In addition, we also easily see that there is a 1D Su-Schrieffer-Heeger (SSH) chain with mirror symmetry localized at each edge, which yields a quantized polarization ±e/2 tangent to the edge.The quantized edge polarization comes from the nontrivial topology of the 1D boundary.The edge polarization of the 1D boundary also indicates the existence of boundary states of boundary, which form the corner modes.
In fact, in the entire topological region, the phase has a quantized bulk quadrupole moment q xy = e/2, quantized edge polarizations p edge x = p edge y = e/2 tangent to the boundaries along x and y, and zero-energy corner states [33].The nontrivial phase is dubbed a quadrupole (topological) insulator.When |γ x /λ x | > 1 or |γ y /λ y | > 1, the model is topologically trivial without boundary signatures.The quantization of the quadrupole moment is protected by the mirror symmetries or chiral symmetry (see Sec. II A 3).The edge polarizations are the polarizations localized along the edges when we take a cylinder geometry with open boundary in x (y) and periodic boundary in y (x) (see Sec. II A 5).The integrated polarization near one edge are quantized to ±e/2 by effective mirror symmetries inherited from the bulk.The edge polarization p edge

Corner modes and corner charges
Here we will follow Refs.[33,34] to give the analytic solution of the corner modes for the BBH model (1).For simplicity, we will solve a continuum model for the lattice Hamiltonian (1) by expanding the Hamiltonian near γ x = γ y = −1 and (k x , k y ) = (0, 0) where the bulk gap closes up to the first order, where we take λ x = λ y = 1 and assume that m x = γ x + 1 and m y = γ y + 1.The parameters m x,y are small and are positive (negative) for the topological (trivial) phase of the BBH model.To solve the corner states under open boundary conditions, we introduce the edges normal to the x or y direction, namely, the x-edge or y-edge forming domain walls and change m x and m y from positive to negative values across the domain walls.Note that for positive m x and m y , the system is in a topological phase while for negative ones, it is in a trivial phase.We first use the continuum model (4) to solve the states localized at the x-edge (x = 0) in the absence of y-edges.The wave function localized at the x-edge takes the form Ψ(x, k y ) = f (x)Φ x (k y ) where f (x) is a scalar function and Φ x (k y ) is a spinor.We require the wave function to satisfy the Schrodinger equation, where we replace k x with −i∂ x and ϵ is the eigenenergy.This equation has a solution f (x) = exp( x 0 m x (x ′ )dx ′ ) and the spinor Φ x should satisfy the equation (Γ 4 − iΓ 3 )Φ x = 0, i.e. (I − τ z ⊗ σ z )Φ x = 0. We have two solutions Φ x1 = (1, 0, 0, 0) T or Φ x2 = (0, 0, 0, 1) T .We can project the Hamiltonian into the subspace of these edge states to obtain the low-energy effective Hamiltonian for the x-edge where µ x,y are Pauli matrices in the basis {Φ x1 , Φ x2 }.We can perform an analogous calculation for the y-edge and have the equation for the spinor (I−τ 0 ⊗σ z )Φ y = 0, which has two solutions Φ y1 = (1, 0, 0, 0) T or Φ y2 = (0, 0, 1, 0) T .Similarly, the y-edge Hamiltonian is obtained by projecting the Hamiltonian into the basis {Φ y1 , Φ y2 }, where γ x,y are Pauli matrices in the basis.
Both of these edge Hamiltonians take the form of massive 1D Dirac models for a 1D topological insulator with the 1D mirror symmetry.Now we consider a corner that is the intersection of the right x-edge and the upper yedge.The zero-energy corner mode can be seen as either a y-domain wall of the x-edge or a x-domain wall of the y-edge, which have the identical solution as Therefore, we have a single corner mode at each corner of a square geometry, which is a simultaneous eigenstate of both edge Hamiltonians.
In general, the boundary states of a HOTI can be understood from the view of domain walls of Dirac mass protected by crystalline symmetries [35][36][37]49].Here we illustrate this picture using a Dirac Hamiltonian form of the 2D BBH model with fourfold rotation symmetry.In the presence of the C 4 symmetry (γ x,y = γ and λ x,y = λ), the BBH Hamiltonian (1) can be rewritten in an alternative form of Dirac Hamiltonian [36], by performing the transformations and Γ 1 → τ 1 σ 2 ; these new matrices constitute a new set of Dirac matrices anticommmuting with each other.Meanwhile, we define the parameters m = √ 2γ and Similarly to the BBH model, this Hamiltonian has a HOTP with zero-energy corner modes protected by C 4 T symmetry when |m/t| < 2 and ∆ ̸ = 0 [36].
Without the mass term, i.e. ∆ = 0, the Hamiltonian (9) has a TRS represented by iσ y κ, belonging to the class AII (time-reversal invariant topological insulator with spin).When |m/t| < 2, gapless helical edge state arise.The τ 2 σ 0 term breaks the TRS as well as the C 4 rotational symmetry represented by e −i π 4 σz , but preserves the symmetry of their combination.The breaking of the TRS leads to the opening of the gap of the helical modes, resulting in an effective mass for the edge Dirac Hamiltonian.Due to the C 4 T symmetry, the Dirac mass of the edge Hamiltonian changes sign under a C 4 rotation.Thus, there exists a domain wall at a corner between two edges related by the C 4 rotation, which gives rise to a corner state.There are a total of four corner states for an open boundary compatible with the C 4 symmetry.
The corner states lead to filling anomaly [71] or fractional charges [57,162] localized at the corners.Specifically, if we require the system to be gapped and symmetric, it is inevitable to have the four degenerate corner states to be either totally occupied or totally unoccupied, resulting in excess fractional charges localized at each corner, which is the manifestation of filling anomaly.On the other hand, if we require the constraints of neutrality rather than symmetry, we can add an infinitesimal symmetry-breaking term in the Hamiltonian to break the degeneracy so that two diagonal corner states are occupied.As a result, fractional corner charges of either e/2 or −e/2 occur near the corners.
In practice, the corner charge for one corner (similarly for the other three corners) is numerically calculated by performing the integration of the charge density over a quadrant of the system [34], Here, is the charge density.The first and second terms arise from the atomic positive charges and the electron distribution, respectively.The latter is calculated by the occupied state |u n ⟩ of the Hamiltonian under open boundary conditions ([u n ] R,α represents the αth component at the site R of the nth occupied state).For the corner charge computation, an infinitesimal term should be added to break the symmetry to lift the fourfold degeneracy of corner modes.In other words, two corner states have positve energy and the other two have negative energy.One thus can fill two negative corner states with electrons at half filling, leading to fractional corner charges of ±e/2 (we consider the existenc of positive charges of 2e in each unit cell).
The topological corner modes of the quadrupole TIs can also be probed by the in-gap modes of the nested entanglement spectrum [37,66] or the entanglement spectrum between a quarter part and its complement [76].In other words, for most models, both the energy spectrum of the Hamiltonian and the entanglement spectrum support in-gap corner modes in the topological region.However, Tao and coworkers propose a model Hamiltonian with a staggered Z 2 gauge field and find that this model leads to a quadrupole insulator without corner modes in the energy spectrum but with in-gap corner modes in the entanglement spectrum [76].In fact, the entanglement spectrum has a one-to-one correspondence with the energy spectrum of the flattened Hamiltonian rather than that of the original Hamiltonian.When we change the flattened Hamiltonian to the original Hamiltonian by varying the eigenenergies, the edge energy gap closes and reopens, leading to the disappearance of the corner modes [76].During this process, the bulk states remain unchanged so that the quadruple moment does not change.

Quadrupole moment as a topological invariant
Here we review a real-space formula of the electric quadrupole moment and show that its value should be quantized in the presence of some crystalline symmetries or chiral symmetry.The formula of the quadrupole moment defined here may not correctly evaluate the physical quadrupole moment for a real system.There have been some discussions on the subtlety of proper definition of the bulk quadrupole moment [163][164][165][166]. Nevertheless, it can be used as an indicator for the second-order topology of quadrupole insulators like the BBH model.
The real-space quadrupole moment is defined by a generalization of the Resta formula for the electric polarization [167].Proposed in Refs.[123,124], the quadrupole moment q xy (in units of e) for a many-electron system is defined as where Qxy = r=(x,y) xy LxLy (n r − n) with nr denoting the electron number operator at site r = (x, y) ∈ (0, L x ] × (0, L y ] in a system of size L x by L y .n denotes the average electron filling per unit cell counting the atomic positive charge contribution that we need to subtract.|Ψ G ⟩ is the many-body ground state of electrons in a periodic system.Similar to the polarization from the Resta formula, the quadrupole moment q xy is only defined in [0, 1) because of the periodicity of complex phases.Note that to have a well-defined quadrupole moment (invariance under translations), we require that the bulk polarizations vanish [33].The octupole moment can be defined analogously by replacing Qxy with Ôxyz = r=(x,y,z) xyz LxLyLz (n r − n) [123,124].For non-interacting fermion systems, the many-body wave function can be written as a Slater determinant of occupied single-particle eigenstates.Then, it can be shown that the quadrupole moment is evaluated by [58,125,126] where we define an n a × n c matrix representing the occupied states of an n a × n a single-particle Hamiltonian under periodic boundary conditions.D = diag{e i2πxj yj /(LxLy) } na j=1 is an n a × n a diagonal matrix [(x j , y j ) is the real-space coordinate of the jth degree of freedom].Here, n c is the total number of occupied states, and n a is the total dimension of the single-particle Hamiltonian.In the calculation, one also needs to deduct the contribution from background positive charge distribution, that is, x j y j /(L x L y ) where (x j , y j ) is the position of the jth positive charge.
We can apply the formula (12) to evaluate the quadrupole moment of typical models of quadrupole insulators.For the BBH model (1) with chiral symmetry, we have q xy = e/2 for the quadrupole topological phase when γ x /λ, γ y /λ are both in the interval (−1, 1).For the SOTI with C 4 T symmetry in Eq. ( 9), q xy = e/2 when |m/t| < 2 and |∆| > 0, which correctly characterizes the topologically nontrivial phase with corner modes.Note that the quadrupole moment as a topological invariant can predict the topological phase transitions between the topological and trivial phases of quadrupole insulators when either the bulk energy gap or edge energy gap closes.In fact, the quadrupole moment characterizes the topology of the flattened Hamiltonian, and thus its topology manifests in the existence of in-gap corner modes in the flattened Hamiltonian or the entanglement spectrum [76].For most models, the original Hamiltonian also supports corner states because it is connected to the flattened version without experiencing any edge energy gap closure.However, there exists a model where quadrupole insulators do not have corner states in the energy spectrum [76].
a. Quantization of q xy protected by crystalline symmetries Under mirror symmetry, M x : (x, y) → (−x, y) or M y : (x, y) → (x, −y), we see that the operator Qxy in Eq. ( 11) transforms into − Qxy .Thus, we have q xy → −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) → (−y, x) also flips the sign of q xy .Thus, the C 4 symmetry enforces that q xy = −q xy mod 1, indicating that q xy is quantized in the presence of the C 4 symmetry.
b. Quantization of q xy protected by chiral 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 (12) protected by chiral symmetry by following Refs.[125,126].Let H be a generic single-particle Hamiltonian in real space preserving chiral (sublattice) symmetry, that is, there exists a unitury operator Π such that ΠHΠ −1 = −H.If we consider an occupied state |ψ n ⟩ with energy E n , then we can obtain an unoccupied state by applying Π to the occupied one, that is, Π|ψ n ⟩ with energy −E n .Since the energy spectrum is symmetric about the zero energy, we consider the half-filling case with n a = 2n c .Then, we have q 4π Im log det D mod 1.The set {Π|ψ 1 ⟩, Π|ψ 2 ⟩, • • • , Π|ψ nc ⟩} therefore constitutes the unoccupied states.These unoccupied states can be represented by The quadrupole moment (12) can be written as Because the unitary matrix D commutes with the chiral (sublattice) symmetry transformation Π with Π † Π = 1, i.e.D, Π = 0, Next, we will prove that 2q xy = 0 mod 1.
Proof.We define an n a × n a unitary matrix It can be easily seen that det D † = det U † t D † U t .Then, we have the following relations In the derivation, the following orthonormal properties are used, 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 th octupole moment is quantized in 3D in the presence of chiral symmetry.Similar conclusion can be generalized to systems with particle-hole symmetry [126].Since the chiral symmetry alone can protect the quantization of the quadrupole moment, we can have quadrupole TIs in disordered systems lacking crystalline symmetry [125,126,141].

Wannier bands and topological invariants
Here we review the topological properties of Wannier bands for quadrupole insulators protected by mirror symmetries.
a. Wilson loop and Wannier band In general, we consider a Bloch Hamiltonian H(k) which has N orb degrees of freedom and N occ occupied bands with Bloch functions |u n k ⟩ for n = 1, • • • , N occ .The Wilson loop is closely related to the position operator projected into the occupied space [167].In the following, we will focus on the Wilson loop along k x , W x,k , and the Wilson loop W y,k along k y can be evaluated similarly.The Wilson loop W x,k starting from the base point k = (k x , k y ) of the loop in Brillouin zone, is defined under full periodic boundary conditions as follows where [F x,k ] mn = u m k+δkx u n k , for δk x = (2π/N x , 0).Since the Wilson loop is unitary in the thermodynamic limit, it can be defined as an exponential of a Hermitian matrix Then, we refer to H Wx (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 [168].
The Wilson loop operator can be diagonalized as where |ν j x,k ⟩ is the eigenvectors with components [ν j x,k ] n for n = 1, • • • , N occ .The Wannier Hamiltonian H Wx (k) has eigenvalues 2πν j x (k y ), j = 1, • • • , N occ , which only depend on k y of the base point k.The eigenvalues ν 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 ν j=± x (k y ) of the Wilson loop W x,k over two occupied bands form a pair of opposite values ν − x (k y ) = −ν + x (k y ) mod 1 due to the mirror symmetry M x , leading to the vanishing of the bulk polarizations.The two Wannier bands ν j=± x (k y ) have a Wannier gap at both ν x = 0 and ν x = 1/2 over the Brillouin zone.Therefore, we can have a 1D gapped Wannier Hamiltonian H Wx (k y ) for k y ∈ [−π, π], which may host nontrivial band topology.
b. Nested Wilson loop and Wannier-sector polarization When the Wannier bands ν j x (k y ) are gapped for k y ∈ (−π, π], 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 ν j x (k y ) are gapped for k y ∈ (−π, π] at both ν x = 0 and ν x = 1/2.Then, we can define two Wannier sectors ν − x ∈ (0, 1/2) and ν + x ∈ (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 ν l x (k y ) and Wilson loop eigenstates |ν l x,k ⟩, the nested Wilson loop is a Wilson loop of one subspace of Wannier bands along k y .Specifically, for the two Wannier sectors ν ± x split by the Wannier gap, we define the Wannier band basis in the Wannier sector ν − x as for l ∈ 1 . . .N W . N W is the number of the Wannier bands in the sector ν − x .This basis obeys w l x,k w l ′ x,k = δ ll ′ .Then, the nested Wilson loop along y for the Wannier sector ν x over the Wannier band basis is defined as where [F νx y,k ] rs = w r x,k+δky w s x,k with r, s ∈ 1 . . .N W over all Wannier bands in the Wannier sector ν x , and δk y = (0, 2π/N y ).The total polarization of the Wannier sector ν x is Im log det W± y,kx , which, in the thermodynamic limit, becomes an integral where Ãνx y,k is the Berry connection of Wannier bands ν x having components where l, l ′ ∈ 1 . . .N W run over the Wannier bands in Wannier sector ν x [33,34].Under the mirror reflection M y and M x , the Wanniersector polarizations satisfy p νx y ≡ −p νx y mod 1 and p νy x ≡ −p νy x mod 1, respectively [33,34].Hence, p νx y and p νy x are quantized to 0 or 1/2.For the topological nontrivial phase of the BBH model (1), we have (p νx y , p ν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 Wx (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 (19) at mirror invariant momenta [33,34].
c. Symmetry constraints of Wannier Hamiltonian Consider a generic Bloch Hamiltonian H(k) with a symmorphic lattice symmetry g: r → D g r, described by gH(k)g † = H(D g 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 where [F k ] mn = ⟨u m k+δk |u n k ⟩ with m, n labeling the occupied bands and δk = (k f − k i )/N .We take N → ∞ in the thermodynamic limit.The Wilson line operator transforms under the symmetry as [33,34] Here a unitary sewing matrix B mn g,k = ⟨u m Dgk |g|u n k ⟩ is defined.The matrix reflects the connection between eigenstates at k with eigenstates at D g k.In particular, a Wilson loop along the contour C at the base point k satisfies Now we focus on the topology of the Wannier band ν x (k y ) from the Wilson loop W x,k in mirror-symmetric systems.The results for ν y (k y ) can be obtained analogously.Under the mirror symmetry M x , we have the transformation for the Wilson loop as where x and −x label the direction that a Wilson loop is obtained.We then define a Wannier Hamiltonian for the Wilson loop W x,(kx=0,ky) .The Wannier Hamiltonian in Eq. ( 17) is a multivalued function due to the logarithm, and thus we need to redefine it relative to a branch cut ϵ [66], where we take log ϵ e iϕ = iϕ, for ϵ ≤ ϕ < ϵ + 2π.With the symmetries, it is proved that the Wannier Hamiltonian should satisfy the following constraints [66], d. Wannier-gap winding number The Wannier bands obey the periodicity from the relation ν x ≡ ν x mod 1, reminiscent of the effective Hamiltonian of Floquet TIs [169,170].It is thus possible that the Wannier bands under open boundary exhibit in-gap modes at both ν x = 0 and ν 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 [42].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 [169].
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 [66].
Specifically, a Wilson line with respect to ϵ is introduced, where W x,kx←0 (k y ) ≡ W (kx,ky)←(0,ky) .It is shown that W ϵ x,2π←0 (k y ) = I Nocc , and one arrives at the fact that W ϵ x,kx←0 (k y ) is periodic versus k x .The symmetry further imposes constraints on the Wilson line with respect to ϵ for k x = π and ϵ = 0, π such that [66] B mx,(π,ky) W ϵ=0 x,π←0 For the quadrupole insulator respecting two anticommuting mirror symmetries with N occ = 2, it can be shown that there exists a basis in which the sewing matrices along the mirror invariant lines take the form [66] As a result, the Wilson loops are 2 × 2 matrices taking forms as The topological invariant for the Wannier gaps at ϵ = 0, π can thus be defined as a winding number 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 Ref. [66].
It has been shown that the change of the winding number can reflect the Wannier gap closing at either ν = 0 or ν = ±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.) describes the edge polarization per unit length along the x (y) direction at the boundaries perpendicular to y (x).The Greek letters β = ±y (α = ±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 Refs.[33,34].Specifically, for a 2D insulator with N x × 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 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 W x for the pseudo-1D system with N occ N y occupied bands.
With the obtained Wilson loop of the pseudo-1D system, we can construct the hybrid Wannier functions as [34] where ν j kx n is the nth component of the jth eigenstate |ν j kx ⟩ for the Wilson loop with the Wannier center ν j x for j = 1, • • • , N occ N y , and |n, k x ⟩ is the nth eigenstates of the pseudo-1D Hamiltonian at k x .One can calculate the density of the hybrid Wannier functions (36) to extract the spatially resolved polarization [34].The density is given by where [u n kx ] Ry,α 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 [34] In the calculation, we need to break the mirror symmetries infinitesimally to split the degeneracy of two Wannier centers at ±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.
= −p edge +y It is quantized to 0 or 0.5 up to an integer.One can similarly calculate the edge polarization along y. b.Edge polarization and Wannier band topology Because of the mirror symmetry, the nontrivial Wannier bands can have mid-gap edge modes at 0 or ±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 ±1/2, these modes give rise to the edgelocalized polarizations along x.Thus, the winding number in Eq. ( 35) 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 Z 2 quantity) with respect to some system parameter γ as [66] Here one uses ∆N q,x to count the number of times that the change of the quadrupole moment occurs from the ynormal edge energy gap closing, as one varies the system parameter γ from γ 0 to γ 1 .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 ∆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 = π for a state localized at a y-normal boundary.If the system is in a trivial phase with p edge x (γ 0 ) = 0 at γ 0 , the formula in Eq. ( 40) can be reduced to by choosing a gauge such that W ϵ=π νx (γ 0 ) = 0 in the topologically trivial phase.It has been shown that the edge polarization calculated based on the formula (41) coincides with the Wannier spectrum and the edge polarization computed using the hybrid Wannier functions in a cylinder geometry [66].

Type-II quadrupole topological insulators
The authors in Ref. [33] derive a general relation among the quadrupole moment q xy , corner charge Q corner ±x,±y and edge polarizations p edge ±x y or p edge ±y x for a 2D square classical system by the multipole expansion of an electric potential in classical electromagnetic theory.The relation reads For the BBH model, the relation reduces to (see Fig. 2).In a quantum system, the relation is enforced by the correspondence that the Wannier band and the edge energy spectrum close their gaps simultaneously [71].The correspondence can be illustrated by transforming the BBH model to the following form [66] where In this model, the mirror symmetries are represented by At the highsymmetric momenta for k x (similarly for k y ), one can write the Hamiltonian as the following direct sum form relative to a basis β consisting of eigenvectors of U Mx .For example, one can choose 3) corresponding to eigenvalues 1 and −1, respectively.In this basis, Suppose that |γ y | < 1 so that the SSH model is in a topologically nontrivial phase.We now decrease γ x across 1.
It is clear to see that the second term in both H 1 and H 2 vanishes at γ x = 1 and k * = π so that both H 1 and H 2 become the SSH model.If we impose open boundaries along y, then there exist two pairs of zero-energy edge states, implying the edge energy gap closing at the y-normal boundary.At the same time, there appear the Berry phase of π for the Hamiltonian H 1 and H 2 , which is equivalent to saying that the Wannier bands ν y (k x ) close their gap at ν y = 1/2 when k x = π.However, the correspondence between the Wannier band and the edge energy spectrum may break down for systems with either particle-hole or chiral symmetry [66] or systems without these symmetries [46].Such a breakdown leads to a distinct type of quadrupole topological insulator (dubbed type-II) with q xy = |Q corner ±x,±y | = |p edge ±y x | = e/2, but p edge ±x y = 0, violating the classical relation (see Fig. 2) [66].A similar discrepancy between the Wannier bands and edge energy spectrum has also been observed in a generalized BBH model [69].
We now follow Ref. [66] to show how the breakdown occurs.Consider a Hamiltonian (the Hamiltonian H II in Ref. [66]) where 2. Schematics of edge polarizations and corner charges.(a) Edge dipole moments exist at all boundaries in type-I QTI and (b) they exist only at the boundaries perpendicular to y in type-II QTIs.Corner charges Q corner = ±e/2 (marked by different colors) appear in both phases.From Ref. [66]. and Here, γ, γ c and b 2 are system parameters with b 2 specifying the strength of the next-nearest-neighbor intercell hopping.This Hamiltonian respects the reflection symmetries U Mx = σ 1 ⊗ σ 3 and U My = σ 1 ⊗ σ 1 and the particlehole symmetry P = σ 3 ⊗ σ 0 κ.
Due to the mirror symmetry, we can write the Hamiltonian at k x = π as When γ = γ c , H NNN (k y ) exhibits zero-energy edge modes under open boundary conditions, indicating that H II closes its edge energy gap at y-normal boundaries.However, the corresponding Berry phase is not equal to π, due to the existence of b 2 .As a result, the correspondence breaks down.The breakdown not only leads to the type-II quadrupole insulator, but also a new insulator with only edge polarizations but without corner modes [66].
Figure 3 shows the phase diagram (reproduced from Ref. [66]) of the Hamiltonian with respect to γ When γ decreases cross γ = 0.13, while the y-normal edge energy gap vanishes, the Wannier band gap remains open at ν y = 1/2, resulting in a type-II quadrupole insulator with nonzero quadrupole moment and nonzero edge polarizations along only one pairs of edges.The specific configurations of edge polarizations can be seen from the Wannier spectra ν x and ν y displayed in the insets: In the type-II phase, only ν x = ±1/2 occurs but not for ν y , while in the type-I phase, both ν x and ν y exhibit isolated eigenvalues at ν x = ±1/2.In the phase diagram, one can also observe a nontrivial phase with only nonzero p edge y but without the quadrupole moment and corner modes, which arises from the Wannier gap closing alone from the trivial phase.Such phases are protected by the topology of Wannier bands rather than energy bands.

B. Chiral-symmetric higher-order topological phases
We have seen that chiral symmetry can protect the quantized quadrupole moment of quadrupole insulators.In the following, we will review a systematic method to construct HOTPs with chiral symmetry and their topological characterizations.

Construction of higher-order topological phases with chiral symmetry
Here we follow Ref. [110] to present how to systematically construct the lattice models for HOTPs protected only by chiral symmetry.We start from the following 2D Hamiltonian in momentum space, where Here, Π 2 i = 1 i , and 1 i is the identity matrix with the same size as H i (k i ) for i = x, y.The Hamiltonian H C (k) thus has chiral symmetry Π = Π x ⊗ Π y .Due to chiral symmetry, the energy spectrum is symmetric about zero energy.When H(k) is gapped at zero energy, it can be characterized by a Z-valued topological invariant given by ν 2D = w x w y , where w i=x,y are 1D winding numbers for H i (k i ) protected by chiral symmetry.The winding numbers characterize the first-order topology of the 1D chiral symmetric Hamiltonian H i=x,y .The topological invariant ν 2D can change when the energy band gap closes at zero energy in the bulk or the edges [110].
The zero-energy corner modes arise when both H x (k x ) and H y (k y ) are in a topological phase with nonzero 1D winding number.Specifically, due to the nontrivial property of The above construction can be generalized to construct 3D Dirac semimetals with chiral symmetry by stacking the 2D Hamiltonian along the z direction.The 3D bulk Hamiltonian can be written as This Hamiltonian respects the same chiral symmetry as the 2D case.If we regard k z as a parameter of the 2D Hamiltonian, then we can have the topological phase transition characterized by the change of ν 2D (k z ), which corresponds to gap-closing points in 3D momentum space.Therefore, we can realize a 3D chiralsymmetric second-order topological semimetal which has gapless points and hinge states at zero energy [110].

Higher-order topological AAH models with chiral symmetry
Now we consider a 2D Aubry-André-Harper (AAH) model with chiral symmetry constructed from 1D offdiagonal AAH models [67,68].The 1D off-diagonal AAH model only contains nearest-neighbor hopping with realspace Hamiltonian where |j⟩ denotes the state at the site j, and the hopping amplitude is modulated as where λ is the amplitude of the modulation whose period is determined by α.The modulation is commensurate when the parameters α = p/q with p and q being mutually prime positive integers.When q is an even number, the system respects a sublattice symmetry as chiral symmetry and can have zero-energy edge modes protected by chiral symmetry [171].In momentum space, the chiral symmetry operator is represented by Π = I q 2 ⊗ σ z , and the bulk Hamiltonian H 1D (k) is a q × q matrix given by where m, n = 1, • • • , q and the momentum k ∈ [0, 2π).For the 1D off-diagonal AAH model with q > 2, in addition to topological zero-energy edge modes, we can obtain nonzero-energy chiral edge modes in the energy spectrum with respect to the parameter ϕ corresponding to Chern bands [171].
According to the previously introduced construction method, the momentum-space Hamiltonian of a 2D lattice model with chiral symmetry constructed from the 1D off-diagonal AAH models in the x and y directions can be written as [68] (60) where and I x is an identity matrix of size q x .Clearly, the 2D Hamiltonian has a chiral symmetry represented by In real space, the above 2D off-diagonal AAH model is defined on a square lattice [68], where |j x , j y ⟩ represents the state at the site (j x , j y ), and the hopping amplitudes along the x and y directions are modulated respectively as where λ x and λ y are the amplitudes of the modulations.The modulation is periodic when the parameters α x = p x /q x and α y = p y /q y with p x and q x (p y and q y ) being mutually prime positive integers, whereas it becomes quasiperiodic when α x or α y is a irrational number.When α x = α y = 1/2, the Hamiltonian can be seen as a 2D lattice model constructed from two 1D SSH models, and can be mapped to the BBH model.As we have seen in the previous section, the 2D off-diagonal AAH model can have zero-energy corner modes if both the 1D off-diagonal Hamiltonian H AAH,i=x,y are topological.In the following, we will follow Ref. [68] to review the results with ϕ x = ϕ y = ϕ for the commensurate case with (α x , α y ) = (1/2, 1/4) and (α x , α y ) = (1/4, 1/4) and the incommensurate case with α x = α y = ( √ 5 − 1)/2.a. (α x , α y ) = (1/2, 1/4) 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 ϕ ∈ (− π 2 , − π 4 ) and ( π 4 , π 2 ).These states are spatially localized at the four corners of the 2D lattice [68].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 AAH,x (k x ) with α x = 1/2 and H AAH,y (k y ) with α y = 1/4, zeroenergy edge states appear when ϕ ), respectively.The regions' intersection thus corresponds to the topological parameter region observed in Fig. 4(a), corresponding to the topological invariant ν 2D = 1 with w x = w y = 1.
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 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 Fig. 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 ϕ = ±π/2 and the xnormal edge gap closing at ϕ = ±π/4.
If one views ϕ as the momentum k z along z, then one obtains a 3D AAH model, which exhibits a 3D higherorder 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 Fig. 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 [68].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 higher-order topological insulator with chiral hinge states [37].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 Fig. 4(a).Such states are known as the bound states in the continuum [172][173][174].The higher-order zeroenergy bound states in the continuum have also been found in a 2D HOTI [175].
b. (α x , α y ) = (1/4, 1/4) In this case, the zeroenergy corner modes appear when ϕ ∈ (− 3π 4 , − π 4 ) and ( π 4 , 3π 4 ), which correspond to the regime with zero-energy edge modes in the 1D off-diagonal AAH model as shown in Fig. 4(c).Similarly, as in the previous case, the corner modes are characterized by the quadrupole moment [see Fig. 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 , ϕ) or (k x , ϕ) 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. c.
The irrational α x and α y lead to a 2D quasicrystal model where the nonzeroenergy chiral modes localized at the edges and corners still exist [68].Similar to the commensurate AAH models, the nonzero-energy corner modes can survive across the continuous bulk spectrum.

The Z-valued multipole invariant protected by chiral symmetry
It is well known that chiral symmetry can protect the first-order topological phases with Z-valued topological invariants in odd dimensions according to the tenfold classification [30].For example, 1D TIs with chiral symmetry can be characterized by a winding number which gives the number of zero-energy edge modes.In fact, from the previous sections, we see that the systematic construction of a 2D system with chiral symmetry can give rise to a higher-order phase that is Z classified.However, the quadrupole moment is a Z 2 topological invariant.The authors in Ref. [72] propose a Z-valued multipole topological invariant to characterize the chiralsymmetric HOTIs.In the following, we will review the integer-valued topological invariant, which extends the Z 2 classification for chiral-symmetric quadrupole TIs.
Let us consider a Hamiltonian H respecting chiral symmetry Π (sublattice symmetry), that is, ΠHΠ = −H.Since Π 2 = 1, there are two eigenvalues of ±1 for Π corresponding to eigenvectors {|ϕ ± n ⟩}, respectively.The chiral symmetry operator is then represented by σ z ⊗I N S in the basis {|ϕ + n ⟩} ∪ {|ϕ + n ⟩}.Under this basis, the Hamiltonian H has the following form where h and h † describe couplings between the states in the two subspaces with opposite chiral charges (eigenvalues of the chiral operator).The chiral-symmetric Hamiltonian (63) can be easily solved by the singular value decomposition (SVD) of the matrix h: where U S for S = ± is an N S ×N S unitary matrix written as and Σ is a diagonal matrix consisting of singular values {ϵ n }.Here, {ψ + n } and {ψ − n } are normalized vectors defined in the + and − subspaces, respectively.The eigenstates of the Hamiltonian can be constructed as which is the eigenstate of H with the eigenvalue ϵ n .Because of the chiral symmetry, for every eigenstate |ψ n ⟩, we have as the eigenstate with the opposite eigenvalue −ϵ n .
For a 1D topological insulator with the bulk Hamiltonian H(k x ), we can define a unitary matrix q = U + U † − from the SVD (64) for each k x in the Brillouin zone.The winding number is thus evaluated by This topological invariant can also be defined in an equivalent real-space form [72], where P S x = U † S P S x U S and P S x is the sublattice dipole operator associated with the Resta formula for polarization, where the 1D system has L unit cells, and |R, α⟩ denotes the orbital α of a unit cell R.
The 1D real-space winding number can be generalized to define the multipole chiral numbers for 2D and 3D systems with chiral symmetry.For example, for a 2D system with chiral symmetry, one can define a quadrupole chiral invariant by where QS xy = U † S Q S xy U S for S = ± with the U S given by the SVD in Eq. (64).Here, Q S xy for S = ± are the quadrupole operators restricted to one sublattice: where |R, α⟩ denotes the orbital α of a unit cell R = (x, y) in a 2D system of size L x × L y .This operator is related to the definition of the Z 2 quadrupole moment in Sec.II A 3. The octupole chiral invariant can be defined for 3D chiral-symmetric systems analogously [72].
As a topological invariant, these multipole chiral numbers are proved to be integer-valued [72].Note that the real-space formula also works for systems without translational symmetry, similar to the traditional quadrupole moment.
Reference [72] illustrates that the quadrupole chiral invariant (69) can characterize the different phases of 2D second-order TIs with chiral symmetry and indicates the number of topologically protected corner states therein.For a model supporting quadrupole insulators, if only the nearest-neighbor hoppings are nonzero, we have a topological phase with the quadrupole chiral invariant N xy = 1 in consistent with the quantized quadrupole moment q xy = 1/2.When the long-range hoppings are increased, the system will undergo a bulk gap closing and transition into another nontrivial phase with N xy = 4 but with q xy = 0.With open boundaries, this new phase possesses four zero-energy states localized at each corner within a single sublattice (see Ref. [72] for details).Therefore, the Z-valued quadrupole chiral invariant can reveal the higher-order topology protected by chiral symmetry beyond the Z 2 -valued quadrupole moment invariant.As shown in Ref. [72], the topological phase transitions between phases of different multipole invariants accompany the gap closing at either the bulk or the boundary.

C. Higher-order topological insulators with chiral hinge states
There exist various types of 3D HOTIs protected by crystalline symmetries [43,49].Here, we review a type of 3D HOTIs with chiral hinge modes introduced in Ref. [37].

Model Hamiltonian
We consider the following 3D Bloch Hamiltonian [37] where τ i , σ i with i = x, y, z, are Pauli matrices describing the orbital and spin degree of freedoms, respectively.When ∆ 2 = 0, the Hamiltonian (71) respects TRS represented by iσ y κ and C z 4 (four-fold) rotational symmetry represented by C z 4 ≡ τ 0 e −i π 4 σz , giving a 3D topological insulating phase when 1 < |M | < 3.In the presence of the ∆ 2 term, both TRS and rotational symmetry are broken individually, but the Hamiltonian respects their combination C z 4 T symmetry, that is, The C z 4 T symmetry is an anti-unitary symmetry with ( Ĉz 4 T ) 4 = −1.
Figure 5 shows the phase diagram of Hamiltonian (71).The system is in a 3D chiral higher-order topological insulating phase when 1 < |M/t| < 3 and ∆ 1 , ∆ 2 ̸ = 0, which hosts the gapless chiral hinge modes in the bulk gap under open boundary conditions along x and y and periodic boundary conditions along z, as shown in Fig. 5 (a) and (c).The chiral hinge states as 1D dissipationless conduction channels will give rise to a quantized longitudinal conductance G = 2e 2 /h along z, which provides the transport signature of the higher-order topology [37].The existence of chiral hinge states can be understood through the domain-wall picture as for the 2D case, which will be detailed in the following.

Hinge modes as domain walls
The hinge states can be seen as the domain-wall states of 2D Dirac Hamiltonians resulting from the sign change of the Dirac mass.Without the ∆ 2 term, the system is in a time-reversal-invariant topological insulating phase with Dirac points on the surfaces.When we turn on the ∆ 2 term, it will yield mass terms to gap out surface Dirac cones on the (100) and (010) surfaces (the (001) surfaces remain gapless).Because of the C z 4 T symmetry, the Dirac mass must change sign between two surfaces related by a C z 4 rotation, so that there must be a chiral mode as a domain wall state along the hinge between the two surfaces.In fact, before the discovery of HOTIs, it has been shown that similar hinge modes due to domain walls can be induced by face-dependent magnetic fields on the surfaces of TIs [78,79].
We now follow Ref. [36] to show how the mass appears and changes its sign based on the effective low-energy theory.Let us consider the Dirac states on the surface parallel to z which deviates from the (y, z) plane by an angle θ [see Fig. 5(d)].We write the bulk Hamiltonian (71) in the coordinates of (k 1 , k 2 , k 3 ), where k 1 = − sin θk x +cos θk y , k 2 = k z , and k 3 = cos θk x + sin θk y , so that k 1 , k 2 are parallel to the surface while k 3 is normal to the surface.We then expand the Hamiltonian to the first order of k 1 , k 2 and the second order of k 3 for M/t ≳ −3 and set t = ∆ 1 = 1 for simplicity, Then, we can solve the surface modes in a semi-infinite region (x 3 ≥ 0) by neglecting the k 1 , k 2 , ∆ terms firstly, where the Dirac mass m (θ) is proportional to ∆(θ) as Therefore, we have m θ + π 2 = − m (θ) so that the mass changes the sign from a surface to another under a C z 4 rotation.

Topological characterization
The authors in Ref. [37]  direction; two neighboring Chern insulators have opposite Chern numbers, so that the whole system does not break the C z 4 T symmetry.As a result, if the initial system has two chiral modes at each hinge, these modes can be removed by attaching the Chern insulators.However, a Z-valued topological invariant (the winding number of the quadrupole moment) is proposed [127,128], making it possible to observe the phase transition from a trivial phase to a topological phase with more than one chiral modes at each hinge if we consider long-range hopping.It in fact remains unclear whether the attached system is topologically equivalent to a chiral HOTI with some winding number of the quadrupole moment.We also want to note that the gapped surfaces of 3D chiral HO-TIs with one chiral mode at each hinge are distinct from Chern insulators in 2Ds in that a single surface carries an anomalous Chern number of ±1/2 rather than an integer measured from the real-space Chern marker [176].
a. Chern-Simons invariant The Chern-Simons invariant is given by where A a is the Berry gauge field defined as , where |u n k ⟩ is the nth occupied eigenstate of the Bloch Hamiltonian.It has been proved that owing to the C z 4 T symmetry, the Chern-Simons invariant is quantized to θ = 0, π mod 2π [37].
The Z 2 topology of 3D chiral HOTIss can be characterized by a Chern-Simons invariant θ, which has been used to characterize the 3D Z 2 TIs protected by TRS [177].When θ = π, the system is in the higher-order topological insulating phase with chiral hinge modes.Moreover, in Ref. [61], the authors propose a generalized Pfaffian topological invariant to characterize the 3D HOTIs with C z 4 T symmetry.b.Nested Wilson loop One can also use the nested Wilson loop to characterize the chiral higher-order topology [37].The Wilson loop along x, W x,k is defined in Eq. ( 16) in Sec.II A 4. For the first-order TIs with TRS, the Wannier spectrum of W x,(ky,kz) is gapless, similar to the energy spectrum under open boundary conditions along x [178,179].For the chiral HOTI in Eq. ( 71), the Wannier spectrum is fully gapped corresponding to the gapped spectrum on the surfaces perpendicular to x, as shown in Fig. 6(a).
It has been found that for the nontrivial phase of the chiral HOTI, the Wannier Hamiltonian has a Chern number C = ±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 Wx (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 Fig. 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 Fig. 6(c)], illustrating the presence of chiral edge modes in the Wannier spectrum.
c.The winding number of the quadrupole moment It has been proposed that the winding number of the quadrupole moment can serve as a topological invariant to characterize the 3D chiral topology [127,128,159].The quadrupole moment winding with respect to k z is defined as where q xy (k z ) is the quadrupole moment defined in Eq. ( 12) for a 2D Hamiltonian at k z .For the model in Eq. ( 71), 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 [128].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 Z-valued, it is possible to observe the phase transition from a trivial phase to a nontrivial phase with W Q = 2.

D. Higher-order topological insulators with helical hinge states
Helical HOTIs host helical hinge modes [see Fig. 6(d)] protected by TRS and crystalline symmetry, such as mirror, rotation, and inversion symmetry [43,44,49].In the case with mirror symmetry, the system has a Z classification [37].In the following, we focus on a model for helical HOTIs as a generalization of the chiral HOTI model (71), which respects both time-reversal and four-fold rotation symmetry proposed in Ref. [36,37].It has a Z 2 classification of bulk topology.

Model Hamiltonian
The model Hamiltonian of the helical HOTI with C 4 and T symmetry takes the form in momentum space as (80) Here, the Pauli matrices τ i and s i act on the orbital degrees of freedom and σ i acts on the spin degrees of freedom.This Hamiltonian respects the TRS and C 4 rotation symmetry represented by T ≡ iσ y κ and C z 4 ≡ s x e −i π 4 σz .When 1 < |M/t| < 3, there exists the nontrivial HOTI phase with helical hinge modes at each of the four hinges.The four helical pairs of hinge states will contribute a quantized longitudinal conductance G = 4e 2 /h along z [37].Analogous to the chiral HOTI, the helical hinge modes can be understood from the perspective of domain walls [36].Without the ∆ 2 term, the Hamiltonian consists of two copies of a strong TI with two Dirac cones on each surface.This system is thus a trivial first-order TI according to the Z 2 classification, whose Dirac cones can be gapped out by a time-reversal invariant mass, such as the ∆ 2 term.The C z 4 symmetry enforces the Dirac mass to have opposite signs between two surfaces so that there exists a helical mode at the domain wall of mass, which can be derived using the low-energy surface theory introduced in Sec.II C 2.
The authors in Ref. [37] introduce a C z 4 -graded Wilson loop for the helical HOTI so that the nontrivial winding of Wilson loop spectra between (k x , k y ) = (0, 0) and (π, π) defines a Z 2 invariant.In addition, they also employ the Wilson loop formalism to characterize the helical HOTI phase, similarly to the chiral HOTI case.Specifically, one can first compute the x-directed Wilson loop to obtain a gapped Wannier spectrum, and then compute the y-directed nested Wilson loop spectrum to obtain a pair of helical modes as a function of k z , reflecting the existence of gapless helical hinge modes.

A Z2 topological invariant for helical HOTIs
In fact, the Hamiltonian in Eq. ( 80) also obeys a Z 2 symmetry of s y with s 2 y = 1, and we thus can divide the Hilbert space into two subspaces with s y = ±1.In each subspace, the quadrupole moment q (s) xy can be evaluated based on Eq. ( 11) so that a winding number of the quadrupole moment about k z in each subspace is defined as [128] It has been proved that due to TRS, so the spin winding number of the quadrupole moment is defined as One can also include a term of t 3 sin k z τ y s x to break the Z 2 symmetry but still preserves TRS and C z 4 symmetry [128].In this case, the authors in Ref. [128] show that for small t 3 , the spin winding number is still applicable.
For generic cases with TRS, a Z 2 invariant ν Q ∈ {0, 1} is derived to characterize the helical HOTI [128,159].It is defined as Here, the k z dependent matrix is defined as where m, n run over both occupied band indices and momenta (k x , k y ) of the Hamiltonian at k z , the operator D is the same as that in the quadrupole moment formula (12), and Pf[•] denotes the Pfaffian of an antisymmetric matrix.For topologically nontrivial (trivial) phase, we have ν Q = 1 (ν Q = 0).Similar to the winding number of the quadrupole moment, this Z 2 invariant can also be generalized to systems without translational symmetry by replacing the momentum k z with a flux twisting through a boundary along z [128].
We now follow Ref. [128] to briefly review the derivation of the topological invariant based on the quadrupole moment for a noninteracting 3D model H h with TRS represented by Using the periodicity of the wave functions about k z and antiunitary property of the time-reversal operator, we can prove the following relations where [A(k z )] T denotes the transpose of the matrix A(k z ).The matrix A(k z ) is antisymmetric at two timereversal symmetric surfaces k z = 0 and k z = π.
We then construct a 1D Hamiltonian based on A(k z ) as Here, we assume that the matrix A(k z ) is invertible so that H A (k z ) has a gap at zero energy.H A respects the chiral symmetry, the TRS and the particle-hole symmetry, belonging to the class DIII for free electrons, and thus its bulk topology is Z 2 classified according to the Altland-Zirnbauer (AZ) classification [30,180].We thus use the Z 2 topological invariant of H A (k z ) to classify a generic 3D SOTI with TRS.The Z 2 invariant ν Q ∈ {0, 1} is given by where the Pfaffian is well defined for A(0) and A(π). Since Note that it is required that the phase of det[A(Φ z )] changes continuously as k z varies from 0 to π.In numerical calculations, one needs to fix the gauge of the eigenstates such that the sewing matrices for TRS take a specific form (see details in Ref. [128]).In Ref. [128], it was also shown that an equivalent formula can be used for practical calculations.Here, q xy (k z ) is the quadrupole moment evaluated at k z .When the pseudospin s y is conserved, the Z 2 topological invariant ν Q and spin winding number W QS are related by (−1) ν Q = exp(iπW QS ) [128].

E. Higher-order topological semimetals
After the discovery of HOTIs, a number of works have explored higher-order topological semimetals with gapless bulk or surface nodes [68,107,[109][110][111][112][113]116].In contrast to surface Fermi arcs in first-order topological semimetals, hinge Fermi arcs appear in higherorder topological semimetals.In the following, we will briefly review how higher-order Dirac and Weyl semimetals arise.
We first consider a model Hamiltonian for higher-order Dirac semimetals proposed in Ref. [107], where γ ν (k z ) = γ ν +0.5 cos k z with ν = x, y.The Gamma matrices are the same as those defined in the BBH model in Eq. ( 1).The Hamiltonian constitutes a stack of the 2D BBH models with the intracell hoppings modified by the hopping along z.Similar to the BBH model, this model respects the spinless TRS T = κ, particle-hole symmetry P = τ 3 κ, chiral symmetry C = τ 3 , reflection symmetries U Mx = τ 1 σ 3 , U My = τ 1 σ 1 and U Mz = I 4 , and inversion symmetry U I ≡ σ 2 .When γ x = γ y = γ, it also has the C z 4 symmetry.
The BBH model exhibits the nontrivial higher-order topology when |γ x | < 1 and |γ y | < 1 [33,34].In the presence of the C z 4 symmetry, the system undergoes a topological phase transition at |γ| = 1 through the bulk energy gap closure, implying the existence of a Dirac point in the bulk at k D z where |γ(k D z )| = 1.For example, in the case with −1.5 < γ < −0.5, there appear two Dirac points at k D z = ± arccos(−2 − 2γ).The zeroenergy hinge modes arise when |k z | < |k D z |, leading to hinge Fermi arcs.Moreover, in the absence of the rotational symmetry, the Dirac points emerge at x-normal or y-normal surfaces.These semimetals are called higherorder Dirac semimetals [107].
As we have proved in Sec.II A 1, the TRS and inversion symmetry ensures the double degeneracy of each band.Since the model in Eq. ( 89) respects both of these symmetries, any touching points must be four-fold degenerate.Therefore, to create Weyl points in the system, one has to add terms breaking one of these symmetries.In addition, the system also respects chiral symmetry so that the Chern number at each k z cannot exist.Therefore, we can add terms to break the TRS so that the chiral symmetry is broken but the particle-hole symmetry is still preserved.In this case, the quadrupole moment is still quantized and the Chern number may appear.For example, let us add a term m 1 σ 2 in the Hamiltonian in Eq. ( 89) that breaks the TRS and chiral symmetry but preserves the particle-hole symmetry.Indeed, this term can split a Dirac point into two Weyl points [see Fig. 7(a)] [111].Between the two Weyl points, the Chern number C = 1 emerges in each k z slice so that there appear surface Fermi arcs connecting the projections of two Weyl points as revealed by Fig. 7(b).The Weyl point thus still corresponds to the anomalous quantum Hall insulating phase transition point.In fact, two of these Weyl points also correspond to the critical point of a higher-order topological phase transition between a quadrupole insulator with nonzero quadrupole moment and a higher-order trivial phase [see Fig. 7(d)].Thus, a hinge Fermi arc also exists connecting these Weyl points as shown in Fig. 7(c).By applying different perturbation terms breaking other symmetries, one can get different configurations of Weyl nodes and Fermi arcs [111].

III. EFFECTS OF QUENCHED DISORDER ON HIGHER-ORDER TOPOLOGICAL PHASES
Disorder is ubiquitous in nature and plays an important role in quantum transports and topological classifications [16,180].For instance, disorder can localize particles, resulting in an Anderson insulating phase [181,182].In the context of symmetry-protected topological phases, these states are usually robust against weak disorder preserving the corresponding symmetry [139,[183][184][185][186][187][188][189][190][191] (or an average symmetry [192,193]).Yet, for a sufficiently strong disorder, a system will transition into a trivial Anderson insulator.Remarkably, it has also been found that disorder can drive a topological phase transition from a trivial phase to a first-order topological phase [139].These TIs emerging from trivial insulators due to disorder are named topological Anderson insulators (TAIs) [139,186].
In the context of HOTPs, it has been found that these phases are stable against weak disorder [33,69,94,125,126,[129][130][131][132][133][134][135][136][137][138].Interestingly, the concept of TAIs is also extended to HOTIs with chiral symmetry [125,126] and higher-order topological superconductors [94], showing that disorder can induce a higher-order topological phase transition from a trivial one to a nontrivial one.Meanwhile, a modified Haldane model has been shown to exhibit a disorder-induced higher-order topological phase from an anomalous quantum Hall insulator [140], which is experimentally observed in electric circuits.In Ref. [194], the authors showed that random flux can induce a metalband insulator transition from a metallic phase to an extrinsic HOTI with zero-energy corner modes in a 2D SSH model.Recently, the disorder-induced third-order topological insulator with the quantized octupole moment in 3D is also found [137].A recent work studied the disorder-induced phase transitions in chiral SOTIs protected by the C z 4 T symmetry, and found that the chiral HOTI is robust with quantized conductance at weak disorder and then driven into a diffusive metal and finally an Anderson insulator [138].
In this section, we will review the disorder effects on quadrupole insulators in subsection III A, on 3D SOTIs in subsection III B, and on HOTSMs in subsection III C.

A. Higher-order topological Anderson insulators
In this subsection, we will review disorder-induced quadrupole TIs, localization properties, self-consistent Born approximation, and effective boundary Hamiltonian from the Green's function.

Disorder-induced quadrupole topological insulators
a. Model Hamiltonian The authors in Ref. [125,126] prove that chiral symmetry can protect the quantization of the quadrupole moment (also see Sec.II A 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 follow-ing tight-binding Hamiltonian proposed in Ref. [125] where ĉ † rν (ĉ rν ) are creation (annihilation) operator of a particle at the νth degree of fredom in a unit cell labelled by the position vector r = (x, y) (x and y are integers).e x = (1, 0) and e y = (0, 1) are unit vectors along x and y, respectively.
describes the intra-cell hopping, and ) describe the inter-cell hopping along x and y, respectively [see Fig. 8(a) for the hopping parameters].The inter-cell hopping magnitude are taken to be t x r = tx r = t y r = ty r = 1, and the other parameters m x r , mx r , m y r , my r 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 [125] and thus share similar topological and localization properties with a disordered BBH model studied in Ref. [126].
The disorder is introduced in the intra-cell hopping, that is, with ν = x, y, where m x = m y = m 0 and V ν r and V ν r 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 x r = mx r = m x and m y r = my r = m y , the system can be described by a Bloch Hamiltonian in momentum space (96) where H ν (k ν , m ν ) = cos k ν σ 1 + (m ν + sin k ν )σ 2 with ν = x, y.This Hamiltonian takes the similar form to the Hamiltonian in Eq. ( 44), and thus can support zeroenergy corner modes under open boundary conditions.When m x = m y , it has the C 4 symmetry so that the quadrupole insulator is intrinsic.Similar to the BBH model, the system is topologically nontrivial (trivial) when b. Topological characterization The higher-order topology of a quadrupole topological insulator 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 Sec. II A 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 [125].For disordered systems without translational symmetries, the polarization can be calculated according to the Resta formula [167,195].In the clean case of the model in Eq. ( 91), 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 Ref. [125].Evidently, the higher-order topological phase 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 higherorder topology is characterized by the topological invariant 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 × 80, P can be numerically evaluated in a system up to 500.
We now introduce the effective boundary Hamiltonian proposed in Ref. [196] and developed in Ref. [125] 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, 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, • • • , 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 G(ω) = (ω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 = G 11 and G N = G N N , which can characterize topologically nontrivial gapped boundaries of HOTIs.
The boundary Green's function can be computed recursively through the Dyson equation [196] Here, G N (G N −1 ) is the boundary Green's function of a system with N (N − 1) unit cells, and g −1 N (E) = ωI − h N refers to the bare Green's function of the N th 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 [196], 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 [196].We can define an effective boundary Hamiltonian from the boundary Green's function, For a quadrupole topological insulator, the effective boundary Hamiltonians for both the x-edges and the yedges 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 [125].This approach has also been applied to 3D third-order TIs with disorder [137].However, in disordered systems unlike the translational invariant systems in Ref. [196], there is no fixedpoint boundary Green's function in Eq. ( 100).Instead, there are fluctuations due to randomness for the effective boundary Hamiltonian H eff over iterations (99).The authors in Ref. [125] 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, where [167] (note that the atomic positive charge contribution should be deducted), x = x xn x with nx 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 |Ψ n ⟩ 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 interation.
c. Disorder-induced HOTIs To see the existence of higher-order TAI phase, the authors consider a topologically trivial phase with m x = m y > 1 in the clean limit.The phase diagram with respect to the disorder strength is mapped out in Fig. 8(b), showing remarkably the disorder-induced phase transition from a trivial insulator to a quadrupole topological insulator.As shown in Fig. 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 Fig. 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 Ref. [125] indicates that q xy will approach 1/2 in the thermodynamic limit.d.Self-consistent Born approximation (SCBA) The disorder-induced topological phase transition can be understood based on the renormalization of a Hamiltonian due to disorder [186].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 [186], disorder-induced 3D quantum anomalous Hall state from a Weyl semimetal [190], and other disorder-induced phase transitions [191,197].It has also been employed to explain the higher-order TAIs and identify their phase boundaries [125,126].Figure 8(d) shows the phase boundaries of the disordered model in Eq. ( 91) 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 [186,198].For a Hamiltonian H = H 0 + V (r) consisting of a clean system H 0 with translational symmetry and a disorder perturbation V (r), the disorder effect at the energy E can be described by the self-energy Σ defined by where ⟨•⟩ denotes the disorder ensemble average.Upon the disorder average, the system can still be described by a translational invariant Hamiltonian, where the selfenergy renormalizes parameters or create new terms.To evaluate the self-energy, we use SCBA where only rain- bow diagrams in the Green's function are included, leading to a self-consistent equation.Specifically, we consider the disordered model in Eq. (91), where the random intra-cell hopping term at each unit cell r reads , and U 1 = σ y ⊗σ z , U 2 = σ y ⊗σ 0 , U 3 = σ 0 ⊗σ y , U 4 = σ z ⊗σ y .For the independent random distributions considered, we require for i, j = 1, 2, 3, 4, where ⟨• • • ⟩ represents the average over disorder configurations.One thus can write the effective Hamiltonian at E = 0 as Here, with disorder, one can evaluate the self-energy Σ based on the following self-consistent equation [125] Σ where The integral is over the first Brillouin zone.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 [125] where Σ 0 , Σ x , Σ y take real values.These new terms arising from disorder thus renormalize the mass terms m x and m y as Σ x and Σ y thus cause topological phase transitions, since the phase boundaries are determined by m ′ x = 1 and m ′ y = 1.
When disorder is weak, one can approximate the selfenergy Σ by taking Σ = 0 on the right-hand side of Eq. ( 109), where with ν = x, y.For m x > 1 and m y > 1, the integrands are positive so that both Σ x and Σ 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 < 1 and m ′ y < 1 [see Fig. 8(d)].

Band and localization properties
a. Introduction 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 [182].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 [181].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 [181].
A significant contribution to the understanding of localization is the one-parameter scaling theory introduced by Abrahams et al. [199].The key concept of this theory is the scaling function defined as where g is the dimensionless conductance and L denotes the system size.In this context, g ≫ 1 and g ≪ 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 > g c and g < g c , respectively.When the system size L is increased, the dimensionless conductance g tends to flow towards g = 0 (g = ∞) if it initially resides in the regime g < g c (g > 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,200].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 [182].
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, levelspacing 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 where ξ represents the localization length.The localization length can be computed by partitioning the Hamiltonian into layers and employing the formula [201] 2/ξ = − lim Here, G 1N is a submatrix of the Green's function G(ω) = (ωI − H) −1 for a Hamiltonian H with N layers [Eq.( 98)], which can be calculated by the iteration For two (or three) dimensions, the normalized localization length Λ = ξ/L is commonly used, with L (or L 2 ) 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 [202] r (118) where δ 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 [202,203].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 [202,203].
Another quantity to examine the localization property is the IPR defined by where ν represents the internal degree 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 ∝ 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 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 2 such that I ∝ 1/L D2 , with D 2 = 0 and D 2 = d corresponding to localized and extended states, respectively.When 0 < D 2 < d, the states are neither exponentially localized nor fully extended, but exhibit multifractal behaviors [204].b.Localization properties at the topological transition point When m x = m y in the model in Eq. ( 91), the disorder-induced higher-order topological phase transition occurs through the bulk energy gap closure, as shown in Fig. 8(c).In the ensemble of disordered systems, there in fact emerges an average C 4 symmetry, that is, a Hamiltonian Ĥ for a certain disorder configuration occurs with the same probability as its symmetry conjugate partner

Ĉ4 Ĥ Ĉ−1
4 .The average symmetry enforce the bulk energy gap closure (see Sec. IV A for detailed discussion).When m x ̸ = m y , the symmetry is broken so that the topological phase transition occurs via the edge energy gap closing [126].
At the critical point, the normalized localization length Λ x = λ x /L y (similarly for λ y /L x ) at zero energy becomes scale free, indicating the divergence localization length in the thermodynamic limit, as shown in Fig. 9(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.
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 Sec. III A 2).The LSR of zero-energy states at the criticality is close to the value 0.53, corresponding to extended states [see Fig. 9(c)].The fractal dimension D 2 from the scaling of the IPR is close to the bulk dimension 2, indicating the states are extended in the bulk [see Fig. 9(d)].Therefore, the topological phase transition accompanies a localization-delocalization-localization (LDL) transition with the bulk gap closing.
For the model studied in Ref. [126], the disorder-driven phase transition between topologically trivial and non- c.The gapless regime For sufficiently strong disorder, the system becomes gapless, leading to gapless higher-order TAI, trivial-II phase and the Griffiths regime [see Fig. 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 Fig. 9(a).The localized nature of the states around zero energy can also be evidenced by their relatively large IPR, the fractal dimension D 2 around zero and the LSR approaching 0.386 [see Fig. 9(c) and (d)].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 higherorder topological phase has also been found [134] 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 Fig. 9(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 Fig. 9(b)].The multifractal behavior at zero energy can be verified by the LSR close to 0.53 [see Fig. 9(c)] and the fractal dimension D 2 between 0 and 2 [see Fig. 9(d)].Similar localization properties have also been observed in the case with m x = m y < 1.

B. Disordered 3D second-order topological insulators
In 3Ds, SOTIs usually support chiral hinge modes or helical hinge modes [35][36][37].Note that it has also been found that non-chiral hinge modes can appear in a higher-order Klein bottle topological insulator protected by momentum-space glide reflection symmetry [122].
Here, we will present the impact of onsite disorder on chiral hinge modes in a 3D system with reflection symmetry [131] and helical hinge modes in a 3D weak SOTI system [133].We will discuss 3D SOTIs in amorphous lattices and structural disorder induced higher-order topological phase transition in Sec.IV A.
For a 3D SOTI protected by reflection symmetry [35], the authors in Ref. [131] have studied the effects of onsite disorder (white noise) on the topological insulator.Note that such a white noise breaks the reflection symmetry for each sample, but the symmetry should be preserved on average (see the discussion in Sec.IV A).It is found that disorder will drive a series of transitions from a SOTI to a 3D first-order topological insulator, then to a diffusive metal, and finally to an Anderson insulator as shown in Fig. 11.Specifically, at weak disorder, the SOTI has the quantized conductance of G = e 2 /h contributed by robust chiral hinge states, although disorder breaks the reflection symmetry.
Meanwhile, due to the linear dispersion of chiral hinge states, the SOTI phase can be characterized by a constant density of states at weak disorder.As the disorder increases, the surface gap closes at a critical point.As a result, the system transitions into the first-order topological insulating phase characterized by the dominate local density of states on surfaces.At a higher critical disorder, the bulk energy gap closes and the system turns into a diffusive metallic phase with nonzero density of states at the Fermi level.Finally, when the disorder is very strong, the system undergoes a 3D metal-insulator transition to an Anderson insulator.The phase boundaries can be obtained by the SCBA method as well [131].
For a 3D weak SOTI with TRS in the clean limit, a first-order topological invariant of the valley Chern number is used to characterize the helical hinge states [133].With disorder respecting TRS, the valley Chern number is not well defined due to the absence of translational symmetry, and the conductance is no longer quantized due to inter-valley scattering [133].Yet, the authors use the distribution of a state on hinges to characterize the hinge states and find that its value remains almost unchanged for weak disorder, suggesting that the hinge states are robust against weak disorder.As disorder increases, the system transitions from the weak SOTI to a weak first-order topological insulator, then to a diffusive metal, and finally turns into an Anderson insulator at strong disorder.

C. Disordered higher-order topological semimetals
In Sec.II E, we introduce how higher-order Dirac and Weyl semimetals arise from the stacking of the BBH model.Here, we discuss how disorder affects the phase diagram.Consider the Bloch Hamiltonian [40] where m, t x , t y , t z , B, λ and ∆ are system parameters.When B = 0, the Hamiltonian respects chiral symmetry represented by τ x σ z so that a higher-order Dirac semimetal can emerge.The Bτ 0 σ z term breaks the chiral symmetry while preserves the particle-hole symmetry represented by τ y σ y κ, splitting a Dirac points into a pair of Weyl points.The Hamiltonian can also realize a higher-order Weyl semimetal with only a pair of Weyl points.
The authors in Ref. [135] study the effects of onsite disorder on the higher-order Weyl semimetals and map out its phase diagram with respect to the disorder strength W and the parameter m as shown in Fig. 12.Five phases are found, including higher-order Weyl semimetals, Weyl semimetals, diffusive anomalous Hall metals, normal insulators and Anderson insulators, by identifying their normalized localization length calculated by the transfer matrix method.Since the disorder breaks the particle-hole symmetry, the quadrupole moment is fragile to weak disorder.They also show that the hinge charge is unstable against disorder.However, numerical results show that the Hall conductance remains unchanged in the semimetal phase, indicating the stability of Weyl points against disorder.The stability of hinge modes is also demonstrated using the neural network method.
From the phase diagram, we see that disorder can induce a phase transition from a Weyl semimetal to a higher-order Weyl semimetal or from a normal insulator to a Weyl semimetal.The SCBA method is also utilized to determine the phase boundaries at weak disorder.

IV. HIGHER-ORDER TOPOLOGICAL PHASES IN NON-CRYSTALLINE SYSTEMS
We are now in a position to explore HOTPs in noncrystalline lattices, including amorphous lattices, quasicrystalline lattices and hyperbolic lattices.
where ĉ † r = (ĉ † r,1 , ĉ † r,2 , ĉ † r,3 , ĉ † r,4 ) with ĉ † r,α being a creation operator of an electron of the αth component at the site of position r, the onsite term T 0 = τ 3 σ 0 (M + 2t 2 ) and the hopping matrix with |d| and θ denoting the distance between two sites and the polar angle of the separation vector, respectively.
which is equivalent to the Dirac Hamiltonian in Eq. ( 9).
To study the topological phases in amorphous lattices, one considers completely randomly distributed sites in a square.It is found that the Hamiltonian in Eq. ( 121) on amorphous lattices can still host four corner states and its quadrupole moment is quantized at q xy = 0.5, indicating the existence of HOTI on amorphous lattices [141].
Based on this model, density-driven higher-order topological phase transitions are observed in Ref. [205].Specifically, the authors map out the phase diagram in the (ρ, M ) plane (ρ is the density of sites and M is the mass) based on the quadrupole moment q xy [see Fig. 13(a)].The figure clearly illustrates that a system in a trivial phase at low density (e.g.ρ = 0.24) will become a nontrivial one at high density (e.g.ρ = 1).Figure 13(b) and (c) also illustrate that as the density increases, the bulk energy gap closes and reopens at the critical density ρ ≈ 0.41 and the quadrupole moment changes from 0 to 1, resulting in the corner modes.
Recently, the authors in Ref. [158] find that amorphous systems can support higher-order topological phases without crystalline counterparts protected by an average C p T symmetry.For instance, the HOTIs can host eight or twelve corner modes, which cannot exist in crystalline systems.Specifically, they consider the Hamiltonian with T 0 = m z τ 3 σ 0 and  To characterize the topological properties, the authors build a Z 2 topological invariant based on the C p M symmetry by imposing twisted boundary conditions [158].In addition, since the model has the chiral symmetry represented by τ 1 σ 3 , one may ask whether the topological phase can be characterized by the traditional quadrupole moment defined in Eq. ( 12).However, when directly evaluating the quadrupole moment, it consistently yields zero.This may be attributed to the presence of an even number of corner states in each quadrant within a nontrivial phase.To address this issue, the authors propose a solution whereby the positions of sites within a 1/p sector and its opposite counterpart are relocated to the first and third quadrants, respectively.Likewise, the sites in other sectors are relocated to the second and fourth quadrants, respectively [158].
It is found that the quadrupole moment based on the transformed site positions and original eigenstates correctly characterizes the topological phase.In addition, it is shown that the higher-order topology in 2D quasicrystalline lattices can also be characterized by the generalized quadrupole moment [159].In 3Ds, the winding number of the generalized quadrupole moment and a Z 2 invariant based on transformed site positions are defined to characterize the chiral and helical modes without crys-talline counterparts [159] (see Sec. IV B for a detailed discussion).

Three dimensions
In Sec.II C and II D, we present a chiral and helical model in 3D crystalline systems to show the existence of chiral and helical hinge modes.There, we see that the crystalline symmetry such as the C z 4 T symmetry is required to protect the quantization of the Chern-Simons invariant, implying that the topological invariant is not well defined if we consider amorphous lattices, which break the symmetry.Fortunately, the winding number of the quadrupole moment is not dependent on any symmetry, making it possible to realize chiral hinge modes in amorphous lattices.Indeed, the authors in Ref. [128] predict the existence of SOTIs without TRS [see the hinge states in Fig. 14(a)] and with TRS in amorphous lattices and structural disorder induced second-order topological phase transition.
a. Amorphous second-order topological insulator The tight-binding Hamiltonian Ĥc studied in Ref. [128] includes the onsite mass term M τ z σ 0 and the hopping matrix  in Eq. (71).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 Φ 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 Φ z [128], defined as where q xy (Φ z ) is the quadrupole moment in the (x, y) plane, which is a function of the inserted flux Φ z .The flux can also be added by introducing a phase factor exp(iΦ z d z /L z ) to the hopping term.In a SOTI phase, a nontrivial winding number of q xy (Φ z ) arises when one tunes Φ z from 0 to 2π, as illustrated in Fig. 14(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 Fig. 14(c)].
In the presence of TRS, the authors construct a tightbinding Hamiltonian [128]  where ), the onsite mass term T 0 = M τ z s 0 σ 0 and the hopping matrix x − d2 y )τ y s y σ 0 + it 3 dz τ y s x σ 0 (128) with {s ν } denoting another set of Pauli matrices.The TRS represented by T is respected, i.e.T Ĥh T −1 = Ĥh , where T i T −1 = −i and T ĉr T −1 = U T ĉr with U T = iτ 0 s 0 σ y .On a cubic lattice, the Hamiltonian reduces the model in Eq. (80) with an extra term proportional to t 3 sin k z τ y s x .
When t 3 = 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 Φ z , where q (s) xy (Φ z ) is the spin quadrupole moment as defined in Sec.D 2 but with the momentum k z replaced by the flux Φ z .As a result, the HOTP can be characterized by the spin winding number defined in Eq. (82).It is also shown that for small t 3 , the spin winding number is still applicable even though s y is not a conservative quantity.Besides, a Z 2 invariant is defined to characterize the topological phase [see Eq. ( 87) in Sec.II D 2 for its definition.Here, the momentum k z should be replaced by a flux Φ z ].
Numerical calculations reveal a range of M with sample averaged Z 2 topological invariant ν Q and spin winding number W QS quantized at 1 up to finite-size effects, as shown in Fig. 15(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.b.Structural disorder induced second-order topological insulator Structural disorder refers to position randomness of lattice sites, serving as a bridge between a crystalline system and an amorphous system.The authors in Ref. [128] 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 Fig. 15(a).The driven SOTI also exhibits the hinge states as revealed by the local DOS.
c. Average symmetry protected bulk energy gap closure In Fig. 14(d) and Fig. 15(a), we see that the higherorder topological phase transition occurs both through the closure of the bulk energy gap, although the Ĉz 4 T symmetry is broken in each sample.Here, we will follow Ref. [128] to show that such a bulk energy gap closure happens due to the Ĉz 4 T symmetry on average.One can also apply similar argument to the 2D case in Sec.III A.
To demonstrate that the Hamiltonian Ĥc respects the average Ĉz 4 T symmetry, one can calculate the symmetry counterpart of the Hamiltonian, 16. Schematic illustration of a system including two subsystems Ĥ1 and Ĥ2 on the site configurations of S1 and S2, respectively.From Ref. [128].
In the derivation, the results are used: Ĉz 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 [192,193].For instance, a configuration S and its rotational part S ′ = D Ĉ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 Ĥc is equal to that of its Ĉz 4 T conjugate partner Ĉ4 T Ĥc ( Ĉ4 T ) −1 so that the ensemble respects the Ĉ4 T symmetry on average.
We now argue that this average symmetry enforces the bulk energy gap closure at the phase transition point for FIG. 17.The spatial distributions of the corner modes of the higher-order topological superconductor on an AB tiling quasicrystal sample.From Ref. [143].
a single system [128].To demonstrate it, we consider a subsystem Ĥ1 whose z coordinates are restricted in the interval [z 0 , z 0 + ∆z] inside a large system (see Fig. 16).As a single system with a spatial configuration S 1 , the subsystem does not respect the Ĉz 4 T symmetry so that the topological phase transition happens through a surface energy gap closing, e.g. the energy gap on the xnormal 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 Ĉ4 S 1 + αe z , which is obtained by rotating S 1 about z by 90 degrees and shifting the coordinate along z by α.This leads to another subsystem Ĥ2 with the spatial configuration S 2 (see Fig. 16).It is clear to see that Ĥ2 is the Ĉz 4 T conjugate partner of Ĥ1 .Therefore, when the x-normal surface energy gap of Ĥ1 vanishes, the y-normal surface energy gap of Ĥ2 must close.It indicates that for the entire system including both Ĥ1 and Ĥ2 , the closure of the energy gap must simultaneously occur on both surfaces, suggesting that the bulk energy gap should close.
Higher-order topological superconductors with Majorana zero modes localized at eight corners of an octagonal sample can occur on a 2D Ammann-Beenker (AB) tiling quasicrystalline lattice with eight-fold rotational symmetry (see Fig. 17) [143].The authors consider the following mean-field real-space many-body Hamiltonian on a AB tiling quasicrystalline lattice, where Ψ † j = (ĉ † j,↑ , ĉj,↑ , ĉ † j,↓ , ĉj,↓ ) with ĉ † j,ν (ĉ j,ν ) creating (annihilating) an electron with spin ν at site j, the mass term H j = µσ z τ z , and the hopping matrix where ∆ is the p-wave pairing strength, θ jl is the polar angle between site j and l, and ⟨• • • ⟩ denotes a pair of sites connected by an edge on the lattice.This Hamiltonian can be written as Ĥ = Ψ † H BdG Ψ where H BdG is the Bogoliubov-de-Gennes (BdG) Hamiltonian and Ψ † = The system respects the particle-hole symmetry, i.e.PH BdG P −1 = −H BdG where P = τ x σ 0 Iκ, and the chiral symmetry, i.e.ΠH BdG Π −1 = −H BdG with Π = τ 0 σ x I (I is an identity matrix acting on the degrees of freedom of lattice sites).
When V = 0, H BdG is a model consisting of two Chern bands with opposite Chern number, thereby contributing helical Majorana edge modes.In fact, the helical edge modes are protected by a time-reversal like symmetry represented by σ y τ x Iκ.The V term breaks the symmetry and gaps out the edge states by introducing an effective mass on each surface; such masses change sign across each neighboring edges, giving rise to Majorana corner modes, as discussed in Sec.II A 2. While the V term also breaks the C 8 symmetry, the C 8 M symmetry is preserved with M = σ z τ 0 I [143].
To define a topological invariant, the authors construct an effective Hamiltonian with G = lim η→0 (H BdG + iη) −1 being the zero-energy Green's function and |k, n⟩ being the plane-wave state for orbital n.Since the system respects the C 8 M symmetry, one can restrict H eff at the C 8 -invariant momenta Γ = (0, 0, 0, 0) and Π = (π, π, π, π) to the eigenspace of C 8 M with eigenvalues ω ±n = exp(±iπn/8) (n = 1, 3, 5, 7).At each subspace, the restricted Hamiltonian still respects the particle-hole symmetry and thus is classified as the zero-dimensional class D. One thus can calculate the sign (denoted by ν n,k * ) of the Pfaffian of an antisymmetric Hamiltonian transformed from the original Hamiltonian.A Z 2 invariant is thus defined as where ν n = ν n,0 /ν n,Π .It is found that in the nontrivial phase with corner modes, ν = −1 whereas in the trivial one, ν = 1.
HOTIs can also arise in quasicrystals [144].The authors in Ref. [144] study a model Hamiltonian with the onsite term T 0 = (M + 2t 2 )σ 0 τ 3 and the hopping matrix on an AB tiling quasicrystal where f (d) = e 1−d/ξ describes the decay of the hopping with respect to the distance d between two sites.Similarly, when g = 0, the Hamiltonian describes a quantum spin Hall insulator protected by TRS represented by iσ y κ and conserved spin σ z [144], hosting helical edge modes.The g term breaks the TRS so as to open an energy gap at boundaries.With this term, the system still respects the particle-hole symmetry represented by σ 3 τ 1 κ and chiral symmetry represented by σ 2 τ 1 .Approximately, the effective mass at each edge depends on cos(ηθ edge ) where θ edge is the polar angle of a boundary.For η = 2, the mass's sign is positive for θ edge ∈ [−π/4, π/4] ∪ [3π/4, 5π/4] (labeled by red color) while it becomes negative for θ edge ∈ [−3π/4, −π/4] ∪ [π/4, 3π/4] (labeled by blue color), as shown in Fig. 18(a) and  (b).For a quasicrystal lattice with pentagon boundaries, shown in Fig. 18(a), (e) and (i), there appear four zeroenergy states localized at corners where the mass exhibits a sign change.For an octagon, the masses are zero at four edges while the sign of the mass at its two neighboring edges is opposite, leading to edge states extended along the four edges [see Fig. 18(b) and (j)].In fact, these edge states can be understood as corner states as one can continuously deform the boundary into a vertex so that a square is formed and corner states are localized at the corners.Similarly, when η = 4, the system exhibits eight-fold symmetric zero-energy corner states for octagon boundaries [see Fig. 18(k)] and for pentagon boundaries, only two corner states appear [see Fig. 18(l)].The HOTI can be characterized by the Z 2 topological invariant ν [143].
The authors in Ref. [144] also construct a Hamiltonian on the AB tiling quasicrystalline lattice with the onsite term T 0 = γ(Γ 2 + Γ 4 ) and the hopping matrix where Γ 4 = τ 1 s 0 and Γ ν = −τ 2 s ν with ν = 1, 2, 3.The system also has chiral symmetry represented by τ 3 and thus the quadrupole moment q xy is a quantized quantity [125,126].Indeed, it is found that in the nontrivial phase q xy = 1/2, and four zero-energy corner states emerge in a geometry with square boundaries [144].However, for the case with eight or more corner states, this quadrupole moment vanishes, and one needs to calculate the generalized quadurpole moment based on transformed lattice sites [158,159].
We have shown that an even number of corner modes can exist in a quasicrystalline lattice.One may wonder whether five corner modes can appear in a Penrose tiling quasicrystalline lattice as the lattice has the five-fold rotational symmetry.Since the sign change of the mass can only give rise to even number of corner modes, one may consider other mechanisms.Indeed, the authors in Refs.[147,149] find that a Penrose tiling quasicrystal can support five corner modes associated with a fractional charge Q = e/5 as shown in Fig. 19(a).To generate the corner modes, they apply an in-plane Zeeman field to a TI with helical edge states on a quasicrystal to gap out the edge states, leading to an effective Hamiltonian at a boundary where ϕ j = ζ j + θ − π/2 with ζ j being the orientation angle of the jth edge [see Fig. 19(b)], k 1 is the momentum along the edge, θ is the polar angle of the Zeeman field.
It has been shown that even though a mass kink ∆ζ across neighboring boundaries does not lead to the sign change of the mass, it can still result in a corner state associated with a fractional charge of Q = e|∆ζ/2π| [208][209][210].It is found that for a pentagonal Penrose tiling quasicrystal, five corner modes appear associated with a fractional charge Q = e/5 due to the angle change of adjacent edges ∆ζ = 2π/5 [149].While the HOTIs in 2D quasicrystals indicate the existence of new TIs without crystalline counterparts, it remains unclear whether novel chiral and helical HOTIs without crystalline counterparts can occur in 3Ds.Very recently, the authors in Ref. [159] construct a 3D model by stacking 2D models on the AB tiling quasicrystalline lattices and find the existence of eight gapless chiral hinge modes, leading to the two-terminal conductance of 4e 2 /h.The TI is characterized by the winding number of the generalized quadrupole moment based on transformed lattice sites.In the presence of TRS, it is shown that eight helical pairs of hinge modes arise, resulting in the conductance of 8e 2 /h.The phase is characterized by a Z 2 topological invariant as well as the spin winding number of the quadrupole moment based on transformed lattice sites.In addition, higher-order Weyl-like and Diraclike semimetals in quasicrystals supporting both surface Fermi arcs and hinge Fermi arcs have been theoretically predicted [159,207].

C. Hyperbolic lattices
Recently, there has been a surge of interest in hyperbolic lattices [211][212][213][214][215][216][217][218][219][220][221][222][223][224][225][226][227][228][229], which belong to the category of 2D non-crystalline lattices.Compared with a crystalline system, any rotational symmetry is allowed in a hyperbolic lattice, providing another platform to realize HOTPs without crystalline counterparts.The arbitrary rotational symmetry comes from the fact that a hyperbolic lattice with constant negative curvature can be tessellated by regular p-gons with any integer p.The hyperbolic lattice in fact respects the hyperbolic translational symmetry.In other words, the lattice can be obtained by applying translational operations to the first unit cell on the hyperbolic plane.For example, in the hyperbolic {8,3} lattice, the translational operators are generated by four generators γ 1 , γ 2 , γ 3 and γ 4 [see Fig. 20(a)].These translational operations form a non-Abelian group.One can also project the hyperbolic lattice onto a Poincaré disk as shown in Fig. 20.Since the lattice sites on the disk break the traditional translational symmetry, we classify them as non-crystalline lattices.
It is found that hyperbolic lattices can also give rise to HOTPs without crystalline counterparts due to the rotational symmetry absent in crystalline systems [153,154].Here, we follow Ref. [153] to present how the new phases arise.The authors construct two tight-binding models.For the first model, the Hamiltonian within the where |riα⟩ represents the αth orbital at site i in the rth unit cell on a hyperbolic lattice, and T (θ ij ) is given by Eq. (124).The hopping between sites in the first unit cell and sites in neighboring unit cells (the set S 1 ) is given by (140) with modified polar angle θ(ri),(1j) .The entire Hamiltonian can be derived through the application of hyperbolic translational operations, ensuring that respects the hyperbolic translational symmetry.
The momentum space Hamiltonian H(k 1 , k 2 , k 3 , k 4 ) in a 4D Brillouin zone with k j ∈ [0, 2π] (j = 1, 2, 3, 4) is constructed based on the hyperbolic band theory [216].One can keep two momenta kī, kj out of the 4D momentum space (k 1 , k 2 , k 3 , k 4 ) as parameters and apply the Fourier transformation to the other two momenta, yielding a real-space Hamiltonian on a square lattice.In this way, the quadrupole moment q ij (kī, kj) can be calculated and the averaged quadrupole moment over k i , k j is defined as q ij = 1/(2π) 2 dkīdkjq ij (kī, kj).It is proved that Q 12 = Q 23 = Q 34 = Q 41 and Q 13 = Q 24 due to the C 8 M symmetry.It is found that the average quadrupole moment q 12 is quantized to 0.5 when g is not very large, suggesting the existence of HOTP.However, the realspace calculations do not suggest the existence of robust corner modes as the energy spectrum differs significantly for distinct system sizes.
To find a higher-order phase, they relax the requirement of the hyperbolic translational symmetry and study the following Hamiltonian on the hyperbolic lattice [153], where |iα⟩ represents the αth orbital at site i.This system exhibits zero-energy modes in the energy spectrum for a hyperbolic lattice under open boundary conditions [see Fig. 21(a)].The topological properties are characterized by the corner charge which approaches the quantized value of 0.5.Interestingly, in contrast to the quasicrystal case, there are two energy gaps: the minigap determined by the first nonzero eigenenergy and the edge energy gap determined by the energy where the degeneracy suddenly changes from the eight-fold to the double one [see Fig. 21(b)].Inside the edge gap, there are infinite number of eight-fold degenerate corner modes besides the zero-energy ones [see the spatial distribution of the states in Fig. 21(c)].Although the minigap may decrease to zero in the thermodynamic limit, finite-size analysis suggests that the zero-energy corner modes persist.In addition, HOTIs with 12, 16 or 20 corner modes have been found as shown in Fig. 21(d)-(f).

V. EXPERIMENTAL PROGRESS
In this section, we briefly summarize the progress in experimental realizations of HOTPs for both crystalline and non-crystalline systems.
For experimental aspects, HOTPs have only been found in a few solid state materials.Bulk bismuth (Bi) [230,231], bismuth bromide (Bi 4 Br 4 ) [232,233] and layered T d -WTe 2 [234][235][236][237] are experimentally denonstrated to be time-reversal 3D HOTIs with 1D helical hinge states.The higher-order properties of Bi are protected by the C 3 and inversion symmetry and Bi 4 Br 4 are protected by the C 2 symmetry.T d -WTe 2 is non-centrosymmetric and hosts nonsymmetry indicated higher-order topology.Among them, only Bi 4 Br 4 is really gapped.The hinge states can be observed directly by scanning tunnelling microscopy (STM) [230,231,233] and angle-resolved photoemission spectroscopy (ARPES) [232] and their transport properties can be measured by Josephson interferometry [230,[234][235][236].The spin polarization of the helical hinge states is resolved by the magneto-optic Kerr-rotation measurement in T d -WTe 2 [237].Helical hinge states are also observed in Cd 3 As 2 by Josephson interferometry measurements [238], consistent with the hinge Fermi arc calculated theoretically [113], implying Cd 3 As 2 is a higherorder topological semimetal.Spin-momentum locking of the hinge states is evidenced by spin transport measurements [239].There is also evidence that FeTe 0.55 Se 0.45 might be a higher-order topological superconductor with helical hinge zero modes [240].

VI. SUMMARY AND PERSPECTIVES
In this review, we have introduced the recent progress on HOTPs in both crystalline and non-crystalline systems.Different from conventional first-order topological phases, the HOTPs are featured by the topological boundary states of lower dimension.We review several prototypical examples of HOTPs including 2D quadrupole TIs, 3D SOTIs with either chiral or helical hinge states, and higher-order topological semimetals with hinge Fermi arcs, and introduce their topological boundary states and corresponding topological invariant.It has been demonstrated that some HOTIs are protected by the topological invariants which do not require the presence of crystalline symmetries.Therefore, these HOTPs can exist in disordered systems or non-crystalline systems that break crystalline symmetries.Then, we turn to discuss the HOTPs in the presence of quenched disorder.Similar to first-order topological phases, it has been shown that HOTPs are robust to weak disorder.We further demonstrate that disorder can induce a HOTI from a trivial insulator, which is the generalization of topological Anderson insulator to the higher-order case.On the other hand, structural disorder in irregular lattices can also have effects on HOTPs.It has been found that there exist SOTIs with chiral or helical hinge states in amorphous systems and structural disorder can remarkably drive a topological phase transition between a trivial insulator and a SOTI.In addition, we review the progress on HOTPs in other non-crystalline systems such as quasicrystals and hyperbolic lattices, which can host HOTPs beyond crystalline systems.
We remark that this review has only discussed a few models for HOTPs and many other aspects of higherorder topology have not been covered here.In addition to static systems, the HOTPs in nonequilibrium systems have also been extensively investigated.In Refs.[66,242,[343][344][345][346][347][348][349], people have studied the far-fromequilibrium dynamics of various HOTPs after a sudden quench of the Hamiltonian.Interestingly, the authors in Ref. [66] found that a quadrupole TI after the quench will undergo a topological transition from a trivial phase to a quadrupole topological insulating phase during the unitary time evolution.On the other hand, a number of works studied the Floquet higher-order topology in periodically driven systems [343,.Similarly to the first-order cases, people discovered the Floquet HOTPs with anomalous corner or hinge modes and the corresponding topological invariants were also proposed [369].In the review, we have focused on the HOTPs for noninteracting fermions, which can be effectively described by topological band theory.Several works have stud-ied higher-order symmetry-protected topological phases in interacting fermion or boson systems which exhibit gapless corner or hinge modes [371][372][373][374][375][376][377][378][379][380][381][382][383][384][385][386].Besides HOTIs and semimetals we have discussed, higher-order topological superconductors have also gained much attention for the Majorana modes at corners or hinges, which provide alternative platforms to realize non-Abelian braiding for topological quantum computation in the future [387][388][389][390][391][392][393].
On the topological characterization of HOTIs, there exists an open question about the properties of the quadrupole moment as a topological invariant.It has been found that the quadrupole moment can change when either the bulk energy gap or edge energy gap closes while the quantity itself is evaluated from the bulk wave functions without boundaries.Moreover, the current operator-based definition of the quadrupole moment may have some inconsistent results for some cases as pointed out in Ref. [163].Another problem worth studying is about HOTPs in a system with gauge flux leading to momentum-space nonsymmorphic symmetry.Recently, multipole TIs with corner modes and HOTIs with hinge modes have been found in lattice models respecting the glide-reflection symmetry in momentum space [76,77,122].The momentum-space nonsymmorphic symmetries may enrich the family of HOTPs.
τ 1 σ 0 are mutually anticommuting Hermitian matrices.Here τ j , σ j are Pauli matrices representing the degrees of freedom within a unit cell.We take lattice constants a x,y = 1.The Hamiltonian (1) can be realized by a tight-binding model with four orbitals per unit cell, see Fig.1(c).This model describes spinless electrons on a square lattice with π-flux per plaquette and the hoppings are dimerized along x and y directions.The Hamiltonian H BBH (k) has two energy bands and each band is two-fold degenerate enforced by the TRS and inversion symmetry (see the following proof).The model is an insulator at half filling unless the bulk energy gap closes at |γ x /λ x | = |γ y /λ y | = 1.The Hamiltonian (1) satisfies two mirror symmetries M x and M y represented by U Mx and U My , respectively,

x and p edge y jumps to e/ 2
from zero as we change parameters γ x /λ x and γ y /λ y across the phase boundaries at |γ x /λ x | = 1 with |γ y /λ y | < 1 and |γ y /λ y | = 1 with |γ x /λ x | < 1.This happens either through the edge energy gap closing or the Wannier gap closing of the Wannier bands (see Sec. II A 4).The BBH model can be generalized to 3D [see Fig.1(b)], which is a model for octupole TIs with quantized octupole moment[33].The model hosts eight zeroenergy corner modes.

FIG. 1 .
FIG. 1.(a) Schematic illustration of the 2D BBH model with quantized octupole moment where γ and λ represent intracell and intercell hopping, respectively.The dashed lines represent the hopping terms with negative signs.The numbers indicate the basis for Γ matrices [see Eq. (1)].(b) Schematic of the tight-binding model with quantized octupole moment.(c) The energy spectrum of the BBH model with open boundary conditions along x and y as a function of γ/λ with γx,y = γ and λx,y = λ.The red lines describe four degenerate corner states.From Ref. [33].(d) Schematic illustration of the BBH model in a limiting case with γx = γy = 0, showing how corner states and edge polarizations arise.

5 .
Edge polarizations a. Edge polarization from the Wilson loop The edge polarization p edge β x (p edge α y

FIG. 3 .
FIG. 3. (a)The phase diagram of the Hamiltonian (53) versus a system parameter γ.In the phase diagram, we observe the topologically trivial insulator, the type-I quadrupole insulator, the type-II quadrupole insulator with only edge polarizations along x, and a nontrivial phase with nonzero edge polarizations but without the quadrupole moment and zero-energy corner modes in the orange region.Distinct phases are identified by the quadrupole moment (blue line), edge polarizations (p edge x , p edge y ) illustrated by the Wannier spectrum νx (νy) obtained in a cylinder geometry with periodic boundaries along x (y) and open ones along y (x).In the Wannier spectrum, the isolated edge Wannier centers are highlighted by red solid circles.In addition, the Wannier-sector polarization (p νx y , p νy x ) and the number of edge modes in the Wannier Hamiltonian Nν ≡ (N 0 νx , N π νx , N 0 νy , N π νy ) are displayed.The vertical dashed lines mark out the critical points where the closure of the edge energy gap or the Wannier band gap happens.Reproduced from Ref.[66].

FIG. 4 .
FIG. 4. (a) The energy spectrum of the 2D AAH model in Eq. (61) under open boundary conditions along both x and y directions.The zero-energy and nonzero-energy topological corner modes are described by the blue and red lines, respectively.The nonzero-energy corner modes in a continuous bulk band are highlighted by the white circle.(b) The quadrupole moment with respect to ϕ at half filling.In (a) and (b), αx = 1/2, and αy = 1/4.(c) and (d) show the same quantities as (a) and (b) but for αx = αy = 1/4.Here tx = ty = 1, λx = λy = 0.8.Reproduced from Ref. [68].

FIG. 5 .
FIG. 5. (a) Schematic illustration of chiral hinge modes.(b) The phase diagram of the model (71), where NI represents the normal insulator.When ∆2 = 0, the chiral HOTI becomes 3D topological insulator with TRS.(c) The energy spectrum of the chiral HOTI model in Eq. (71) with open boundaries along x and y and periodic boundary along z.The chiral hinge states are highlighted by the red lines.(d) A surface that deviates from the yz plane by an angle θ. (k1, k2, k3) denotes a new coordinate system with k3 normal to the surface and k1 and k2 parallel to the surface.(a)-(c) are reproduced from Ref. [37].

1 −
with the boundary condition ψ (0) = ψ (∞) = 0.It has two solutions ψ n=1,2 as ψ n (x 3 ) = f (x 3 )u n , where f (x 3 ) = e −λ+x3 − e −λ−x3 up to a normalization factor, and u n=1,2 are two spinors that satisfy τ y σ 3 u n = u n and λ ± = 1 ± √ 2m with m > 0 in the topological regime.By projecting the remaining terms into the subspace of the two surface modes ψ n=1,2 , one finally gets the effective surface Hamiltonian in terms of Pauli matrices in the basis of u n=1,2

FIG. 6 .
FIG. 6.(a) The eigenvalues of the Wannier Hamiltonian HW x (ky, kz) with respect to (ky, kz).(b) The eigenvalues of the nested Wilson loop along ky over the occupied subspace of the gapped Wannier Hamiltonian HW x (ky, kz).(c) The Wilson loop spectrum along y in a slab geometry along x.(d) Schematic illustration of the helical hinge modes.Reproduced from Ref. [37].

FIG. 7 .
FIG. 7. (a)-(c) The energy spectrum of the Hamiltonian HHODSM perturbed by a term of m1σ2.(a) The bulk band structure at kx = ky = 0.A Dirac point splits into two Weyl point whose charges are marked out as positive and negative signs.The energy spectrum (b) at ky = 0 (kx = 0) with open boundaries along x (y) and that (c) with open boundaries along both x and y.(d) The quadrupole moment qxy(kz), edge polarizations Px(kz), Py(kz) and corner charge Qc(kz) with respect to kz.Here, γ = −1 and m1 = 0.4.Reproduced from Ref. [111].

FIG. 8 .
FIG. 8. (a) Schematic illustration of the disordered tight-binding model in Eq. (91).(b) The phase diagram with respect to the disorder strength W mapped out by the topological invariant P and the quadrupole moment qxy.There exist distinct phases including gapped/gapless higher-order TAI, Griffiths phase, and trivial I/II phases, separated by vertical dashed lines.(c) The bulk energy gap versus W .In (b)-(c), mx = my = 1.1.(d) The topological invariant P with respect to W and m0 with mx = my = m0.The topological phase boundary calculated by the SCBA method is highlighted by the red dotted line.Reproduced from Ref. [125].

386 FIG. 9 .
FIG. 9. (a) The normalized localization length Λx = λx/Ly at E = 0 as a function of the disorder strength W for several distinct Ly.(b) The scaling of Λx at different energies when W = 4.6, illustrating the multifractal behavior at zero energy and localized behavior at other energies.(c) The LSR and IPR with respect to W for the eigenstates around zero energy in a system with size L = 500.(d) The fractal dimension D2 as a function of W for the eigenstates near zero energy.In (a) and (d), the vertical dashed lines mark out the boundaries of different phases shown in Fig. 8(b).Here, mx = my = 1.1.Reproduced from Ref. [125].

FIG. 10 .
FIG. 10.The localization-delocalization-localization transition at the boundaries: finite-size scaling of the unnormalized localization length at zero energy along the x and y direction.(a) Λx and (b) Λy as a function of the disorder strength W with open boundary condition along transverse direction.From Ref. [126].
trivial 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 ̸ = m y in the model in Eq. (91).Li et al. performed the finitesize scaling of the localization length at zero energy for a quasi-one-dimensional ribbon with open boundaries in the transverse dimension.As shown in Fig. 10, 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.
Here, t(|d|) = Θ(R−|d|) exp(−|d|/d 0 ) with Θ(R−|r|) being the Heaviside step function to cut off hoppings longer than R. The lattice constant d 0 is taken to be d 0 = 1.For a square lattice including only nearest-neighbour hoppings with t(|d N N |) = 1, the Hamiltonian reduces to where p takes the value of a multiple of 4 andf (|d|) = Θ(R − |d|) exp[−λ(|d|/d 0 − 1)]/2.The lattice sites are positioned completely randomly with a hardcore radius r h = 0.2 in a regular p-gon.Although each configuration of lattice sites lacks p-fold rotation (C p ) symmetry

FIG. 13 .
FIG. 13.(a) The phase diagram of the tight-binding model in Eq. (121) on amorphous lattices with respect to the density ρ and the mass M , which is mapped out by the quadrupole moment qxy.(b) The energy spectrum under open boundary conditions (blue circles) and periodic boundary conditions (red circles) at M = 0. (c) The energy gap (red line for open boundaries and yellow line for periodic boundaries) and the quadrupole moment qxy (blue line) with respect to the density ρ at M = 0. From Ref. [205].

FIG. 14 .
FIG. 14.(a) The density distribution of four states near zero energy in the SOTI phase, exhibiting the hinge modes.(b) The quadrupole moment Qxy with respect to the flux Φz, reflecting the nontrivial winding number of the quadrupole moment.(c) Sample averaged two-terminal conductance G [blue (system size L = 30) and green lines (L = 40)] and winding number of the quadrupole moment W Q [red line (L = 16)] of the 3D amorphous model as a function of the mass M .The black line describes the conductance of the Hamiltonian on a cubic lattice.(d) Sample averaged bulk energy gap (green and blue lines) of the 3D amorphous model with respect to M compared with the results (black line) on a cubic lattice The red line depicts the DOS at zero energy.Reproduced from Ref. [128].

FIG. 15
FIG. 15.(a) Structural disorder induced SOTI, illustrated by the sample averaged two-terminal conductance G (blue and green lines), winding number of the quadrupole moment W Q (black line) and bulk energy gap (red line) of the 3D amorphous model without TRS with respect to the disorder strength W .(b) Sample averaged two-terminal conductance G (blue line) and Z2 topological invariant (red line) of the 3D amorphous model with TRS as a function of M .The black and green lines describe the conductance and Z2 topological invariant of the model on a cubic lattice, respectively.Reproduced from Ref. [128].

4 T
ĉr ( Ĉz 4 T ) −1 = iτ 0 σ y e i π 4 σz ĉD Ĉ4 r and Ĉz 4 T i( Ĉz 4 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 ′ = D Ĉ4 S ≡ {D Ĉ4 r = (−y, x, z) : r ∈ S}.If we consider cubic lattices, then S ′ = S so that Ĉz T Ĥc ( Ĉz 4 T ) −1 ̸ = Ĥc , indicating that the symmetry is absent.

FIG. 18
FIG. 18. (a)-(d) Illustrations of the angle dependence of the sign of the effective masses near a corner of pentagonal and octagonal AB tiling quasicrystal samples.The red and blue color represents the effective mass with opposite signs.(e)-(h) The energy spectrum under open boundary conditions.(i)-(l) The corresponding real-space distributions of zero-energy states.For the first and second columns, η = 2 and for the third and fourth columns, η = 4. From Ref. [144].

FIG. 19 .
FIG. 19.(a) The energy spectrum of a pentagonal Penrose tiling quasicrystal sample subject to an in-plane Zeeman field, exhibiting five corner modes at the energy EF .The inset displays the spatial distribution of the corner states.(b) Schematic illustration of an edge whose polar angle measured relative to the positive x axis is ζ.Reproduced from Ref. [149].

4 FIG. 20 .
FIG. 20.Schematic illustration of the {8,3} hyperbolic lattice projected onto the Poincaré disk.(b) The {8,3} hyperbolic lattice on the Poincaré disk with the first unit cell enclosed by the green curve and γ1, γ2, γ3 and γ4 representing four translational operators.The corner states are highlighted by yellow circles.From Ref. [153].

FIG. 21 .
FIG. 21.(a) The energy spectrum of the Hamiltonian in Eq. (141) with respect to the system parameter g on a hyperbolic {8,3} lattice at epoch 5 with 2888 sites, illustrating the presence of the zero-energy states.(b) The eigenenergies when g = 0.6 at epoch 6.The minigap represents the first nonzero positive energy and the edge gap represents the energy where the degeneracy suddenly changes from the eight-fold to the double one.Inside the edge gap, all states are eight-fold degenerate.(c) The spatial distributions of the wave functions in the vicinity of zero energy over the 1/8 sector on the boundary of the Poincaré disk, suggesting that the eight-fold degenerate states are spatially localized around a corner while the states whose energy is slightly above the edge gap are more like the edge state.The zero-energy local DOS of hyperbolic (d) {12, 3} lattice, (e) {16, 3} lattice and (f) {20, 3} lattice.Reproduced from Ref. [153].
y), there exists a zero-energy edge state |ψ zero i ⟩ spatically residing at one edge.Then, it is easy to see that |ψ corner ⟩ = |ψ zero OBC , which is spatically localized at a corner.Because H OBC i=x,y have |w i | zero-energy modes, we will generally have |ν 2D | = |w x w y | topological corner modes at each corner.