{'fraction_non_alphanumeric': 0.041933605125218404, 'fraction_numerical': 0.01408880383836703, 'mean_word_length': 4.22352600318702, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 12, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Higher-order Topological Anderson Insulators

We study disorder effects in a two-dimensional system with chiral symmetry and find that disorder can induce higher-order topological insulators from a topologically trivial phase. Their topological properties manifest in a topological invariant defined based on effective boundary Hamiltonians, the quadrupole moment and zero-energy corner modes. We find gapped and gapless topological phases and a Griffiths regime. In the gapless topological phase, all the states are localized, while in the Griffiths regime, the states at zero energy become multifractal. We further apply the self-consistent Born approximation to show that the induced topological phase arises from disorder renormalized masses. We finally introduce a practical experimental scheme with topolectrical circuits where the predicted topological phenomena can be observed by impedance measurements. Our work opens the door to studying higher-order topological Anderson insulators and their localization properties.

We study disorder effects in a two-dimensional system with chiral symmetry and find that disorder can induce higher-order topological insulators from a topologically trivial phase. Their topological properties manifest in a topological invariant defined based on effective boundary Hamiltonians, the quadrupole moment and zero-energy corner modes. We find gapped and gapless topological phases and a Griffiths regime. In the gapless topological phase, all the states are localized, while in the Griffiths regime, the states at zero energy become multifractal. We further apply the self-consistent Born approximation to show that the induced topological phase arises from disorder renormalized masses. We finally introduce a practical experimental scheme with topolectrical circuits where the predicted topological phenomena can be observed by impedance measurements. Our work opens the door to studying higher-order topological Anderson insulators and their localization properties.
Disorder plays an important role in quantum transport, such as Anderson localization and metal-insulator transitions [18]. In the context of first-order topological phases, it has been shown that they are usually stable against weak symmetry preserving disorder. But disorder is not always detrimental to first-order topological phases. Ref. [19] theoretically predicted that disorder can drive a topological phase transition from a metallic trivial phase to a quantum spin Hall insulator; topological insulators induced by disorder are called topological Anderson insulators (TAIs) [19,20]. Since their discovery, there has been great interest and advancement in the study of TAIs [21][22][23][24][25][26]. In addition, disorder can drive a transition from a Weyl semimetal to a 3D quantum anomalous Hall state [27]. Remarkably, the TAI has been experimentally observed in a photonic waveguide array [28] and disordered cold atomic wire [29].
Disorder, topology and symmetry are closely connected, which can be seen from classification theories. For example, random matrix theories are classified based on three internal symmetries, explaining universal transport properties of disordered physical systems [30][31][32]. Similarly, the classification of topological phases is made according to these internal symmetries [33]. Among these symmetries, chiral symmetry plays an important role in disordered systems and many peculiar properties have been found, such as the divergence of density of states (DOS) and localization length at energy E = 0 [34][35][36][37][38]. In 2D, first-order topological phases are not allowed in a system with only chiral symmetry. Yet, it has been reported that a second-order topological phase can exist in a 2D system with chiral symmetry [39,40] and thus provides an ideal platform to study the interplay between disorder and topology.
Here we study the interplay between disorder and higher-order topology in a 2D system with chiral symmetry. We prove that the quantization of the quadrupole moment is protected by chiral symmetry irrespective of reflection symmetry. This gives us an opportunity to explore the effects of off-diagonal disorder respecting chiral symmetry. We theoretically predict the existence of a disorder induced HOTI [dubbed higher-order topological Anderson insulator (HOTAI)] with zero-energy corner modes, which is explained by the self-consistent Born approximation (SCBA). We find gapped and gapless HO-TAIs and a Griffiths regime. In the gapless regime, all the states are localized, while in the Griffiths regime, the states at zero energy become multifractal. In addition, we study the disorder effects on a HOTI and show the existence of gapped and gapless topological phases and a Griffiths regime. Finally, we propose an experimental scheme using topolectrical circuits to realize and detect the HOTAI.
Model Hamiltonian.-We start by considering the following higher-order Hamiltonian withĉ † rν (ĉ rν ) being a creation (annihilation) operator at the νth site in a unit cell described by r = (x, y) with x and y being integers (suppose that the lattice constants are equal to one), e x = (1, 0) and e y = (0, 1). Here h 0 depicts the intra-cell hopping, and h x and h y describe the inter-cell hopping along x and y, respectively [see Fig. 1 parameters and supplementary material for their matrix forms]. The system parameters m x r ,m x r , m y r ,m y r , t x r ,t x r , t y r andt y r all take real values. For simplicity without loss of generality, we take the inter-cell hopping magnitude as energy units so that t x r =t x r = t y r =t y r = 1. For a clean system with m x r =m x r = m y r =m y r = m, the system respects a generalized C 4 symmetry [41]. Since the system contains only the nearest-neighbor hopping, it respects chiral symmetry, i.e., ΠHΠ −1 = −H, where H is the first-quantization Hamiltonian and Π is a unitary matrix [41]. But this system breaks both the time-reversal and particle-hole symmetries, because h 0 is complex. In contrast, if we generalize the Benalcazar-Bernevig-Hughes (BBH) model [1] to the disordered case, it still respects the time-reversal, particle-hole and chiral symmetries. However, these two models are connected through a local transformation and thus have similar topological and localization properties [41].
In the clean case, we show that this system supports zero-energy corner modes [41]. In fact, we also prove the quantization of the quadrupole moment [42,43] protected by chiral symmetry [41], indicating that the quadrupole moment can be used as a topological invariant.
To study the disorder effects, we consider the disorder in the intra-cell hopping, that is, m ν r = m ν + W ν V ν r andm ν r = m ν +W νV ν r with ν = x, y, where V ν r andV ν r are uniformly randomly distributed in [−0.5, 0.5] without correlation. Here W ν andW ν represent the disorder strength. For simplicity, we take W Because of the random character, we perform the average over 200-2000 sample configurations for  numerical calculation. Higher-order topological Anderson insulators.-We map out the phase diagram in Fig. 1(b), showing remarkably the presence of HOTAIs. To characterize the HOTAI, we evaluate the polarization p x (p y ) of effective boundary Hamiltonians at half filling. The higher-order topology can be characterized by a topological invariant defined as P = 4|p x p y |. When P = 1, the system is in a higher-order topologically nontrivial phase, and when P = 0, it is in a trivial phase [41].
We now generalize it to the disordered case. Specifically, we evaluate the average polarization of the effective boundary Hamiltonian at the y-normal boundary (similarly for x-normal one) by p x = 1 Ni Ni n=1 |p x,n |,where p x,n = Imlog Ψ n |e 2πix/Lx |Ψ n /(2π) [44] withx = x xn x ,n x being the particle number operator at the site x, and L x being the length of the system along x (we also deduct the atomic positive charge contribution). Here |Ψ n is the ground state at half filling of the boundary Hamiltonian H n = −G 2n (E = 0) −1 with G 2n being the 2nth boundary Green's function obtained by [45,46] where h n is the Hamiltonian for the nth layer and V n−1 is the coupling between the (n − 1)th and nth layer. We note that p x,n is quantized to be either 0 or 0.5 for each iteration since H n also preserves chiral symmetry. The polarization is evaluated at even steps of Green's function given that there are two different layers in the clean limit. In the disordered case, the intra-cell hopping parts in h n and V n are randomly generated for each iteration [41]. The topological invariant P is finally determined.
In Fig. 1(b), we plot the topological invariant P as the disorder magnitude increases. We see that P suddenly jumps to 1 when W ≈ 2.1, indicating the occurrence of a topological phase transition. P remains quantized to be 1 until W > 3.5, where it begins decreasing continuously. This regime corresponds to the Griffiths phase where topologically trivial and nontrivial sample configurations coexist [41]. When W > 6, P vanishes, showing that the system reenters into a trivial phase. To further identify that the induced topological phase is a higherorder topological phase, we calculate the quadrupole moment, which can be used as a topological invariant since its quantization is protected by chiral symmetry [41]. Figure 1(b) shows that the quadrupole moment qualitatively agrees with the results of P . Yet, conspicuous discrepancy can be observed. The quadrupole moment over many samples is not quantized to 0.5 in the regime where P = 1 and the Griffiths regime is much larger. We attribute this to the finite-size effects [47].
The topological phase transition occurs as the bulk energy gap closes and reopens, as shown in the inset of Fig. 1(b). When W is further increased, the energy gap closes again and remains closed due to the strong disorder scattering, leading to the gapless HOTAI. Even in the gapless regime, the topological invariant P can still be quantized as shown in Fig. 1(b). In fact, in this phase, all the states are localized corresponding to an Anderson insulator (see the following discussion).
To further confirm that the TAI is a higher-order topological state, in Fig. 1(c-d), we display the local density of states (LDOS) at E = 0 for two typical values of W corresponding to a gapped and gapless topological phase, respectively, clearly showing the presence of zero-energy states localized at corners. The evidence above definitely suggests the existence of HOTAIs.
Localization properties.-We now study the localization properties of energy bands in different phases by evaluating their localization length, adjacent levelspacing ratio (LSR), inverse participation radio (IPR) and fractal dimensions. The LSR is defined as being the ith eigenenergy sorted in an ascending order and i denotes the sum over an energy bin around the energy E with N E energy levels counted. For localized states, r ≈ 0.386 corresponding to the Poisson statistics and for extended states of symmetric real Hamiltonians, r ≈ 0.53 corresponding to the Gaussian orthogonal ensemble (GOE) [48].
The localization property can also be character-ized by the real space IPR defined as . This quantity evaluates how much a state in an energy bin around energy E is spatially localized. For an extended state in 2D, I ∝ 1/L 2 with L being the size of a system, which goes zero in the thermodynamic limit; for a state localized in a single unit cell, it is one. It is well known that at the critical point between localized and delocalized phases, the state exhibits multifractal behavior with fractal dimensions D 2 defined through I ∝ 1/L D2 [49]. Clearly, D 2 = 2 and D 2 = 0 indicates that a state is extended and localized, respectively, in the thermodynamic limit; intermediate values of D 2 suggests the multifractal state.
In Fig. 2(a), we plot the normalized localization length is the localization length along ν calculated by the transfer matrix method [50]. In the gapless HOTAI and trivial-II phases, we see the decrease of Λ x as L y is increased, suggesting that the states at E = 0 are localized. The decline can also be clearly seen in the inset of Fig. 2 In fact, all states are localized in these two phases [41]. This shows that even in the higher-order case, the topology can be carried by localized bulk states. Being localized for the states in these regimes is also evidenced by their relative large IPR and the LSR approaching 0.386 [see Fig. 2(c)]. In these regimes, the fractal dimension D 2 becomes negative or approaches zero [see Fig. 2(d)], further indicating that the states around zero energy are localized [51]. Figure 2(a) also demonstrates the existence of a regime (corresponding to the Griffiths regime) where Λ x at E = 0 remains almost unchanged as L y increases, suggesting a multifractal phase in this regime. The multifractal phase resides between two localized phases, which is very different from the conventional wisdom that a multifractal phase lives at the critical point between delocalized and localized phases. In fact, only the states at or very near E = 0 become multifracal, and all other states remain localized [see Fig. 2 [41]. The multifractal properties are also evidenced by the fractal dimension of the states around zero energy as shown in Fig. 2(d). In a system with chiral symmetry, the multifractal phase is usually accompanied by the divergence of the DOS [52,53], showing its critical behavior. We also find the development of a very narrow peak of the DOS at E = 0 in this regime [41].
In the gapped regime, there are trivial-I and gapped HOTAI phases. In the trivial phase, the states at the band edge around zero energy exhibit very small IPR and LSR around 0.53 [see Fig. 2(c)], suggesting the delocalized property of these states. Yet, further delicate calculation shows that D 2 has negative values (in the region around W = 1) for the first two energy levels and becomes larger as the absolute values of energies increase corresponding to the increase of the energy band index [see Fig. 2(d)]. This indicates the existence of mobility edges given that the states at the band edges are localized while other states are delocalized [41].
In the gapped HOTAI, Fig. 2(c) and (d) illustrate that the LSR experiences a drop from around 0.53 to 0.386 and D 2 drops from 1.75 to negative values, suggesting that the states at the band edge undergo a phase transition from delocalized to localized ones. In fact, more evidence indicates that in the phase, all states can be localized [41]. But in the regime closer to the phase transition point W ≈ 2.1, the localization properties are not conclusively determined.
Self-consistent Born approximation.-Based on the SCBA [20], the effective Hamiltonian at E = 0 is given by is the momentum space Hamiltonian without disorder and the self-energy Σ is determined by solving a self-consistent equation [41]. For clarity, let us first approximate the self-energy by taking Σ = 0 in the Green's function in the self-consistent equation, yielding where the integral expressions for Σ x and Σ y can be found in supplementary material. The disorder renormalizes the topological masses so that m ν → m ν = m ν + Σ ν (ν = x, y), leading to topological phase transitions when m ν < 1, which happens in our system as Σ x = Σ y < 0. Disorder effects on HOTIs.
-We now study the effects of disorder on HOTIs. Specifically, we consider m x = m y = 0.5 corresponding to a HOTI in the clean limit. We find that the topological phase is stable against weak disorder as evidenced by the quantized topological invariant P in Fig. 3(b). When the disorder strength becomes sufficiently strong, it enters into a Griffiths regime with fractional P and finally becomes a trivial phase. The strong disorder also closes the energy gap when W > 2.6. In the gapless HOTI and trivial-II phases, all states are localized, as evidenced by the normalized localization length, LSR and IPR [see Fig. 3(c) and (d)] [41]. In the Griffiths regime, the states at E = 0 are multifractal and all other states are localized [see Fig. 3(c)] [41]. In the disordered gapped HOTI, we find that for weak disorder, the system exhibits mobility edges with the states around the band edge being localized as shown by the LSR around 0.386 in Fig. 3(d). For larger disorder, all the states become localized in this phase [41].
Experimental realization.-The BBH model has been experimentally realized in several metamaterials, such as microwave, phononic, photonic and topolectrical circuit systems [54][55][56][57]. In fact, some systems, such as silicon ring resonators [57], have demonstrated the robustness of zero-energy modes to certain disorders. The HO-TAI can be easily realized in these systems when the off-diagonal hopping disorder is considered in the experimentally realized BBH model. The BBH model has also been implemented in topolectrical circuits, and zeroenergy corner modes are probed by measuring two-point impedances [56]. This platform can be easily changed to observe the HOTAI by tuning the capacitance of capacitors and inductance of inductors. In supplementary material, we provide a topolectrical circuit scheme to realize our model (1). We also show that the induced zeroenergy corner modes can be measured by the two-point impedance as performed in the experiments [see In summary, we have discovered the HOTAI in a 2D disordered system with chiral symmetry. Specifically, we show that a topologically trivial phase can transition into a HOTI, when disorder is added. We find gapped and gapless HOTAIs and a Griffiths regime. In the gapless HOTAI, all the states are localized, while in the Griffiths regime, the states at zero energy are multifractal and other states are localized. The Griffiths regime corresponds to a critical regime between two localized phases: a gapless HOTAI and trivial-II phase. Our results demonstrate that disorder can induce higher-order topologically nontrivial insulators with peculiar localization properties from a trivial phase and thus opens a new avenue for studying the role of disorder in higher-order topological phases.
We thank Y.-L. Tao and N. Dai for helpful discussion.

SUPPLEMENTAL MATERIAL
In the supplementary material, we will provide the hopping matrices in the Hamiltonian (1) in Section S-1, show the existence of zero-energy corner modes in the clean case in Section S-2, prove the equivalence between the Hamiltonian and the disordered BBH model in topological and localization properties in Section S-3, present the generalized C 4 symmetry in Section S-4, prove the quantization of the quadrupole moment protected by chiral symmetry in Section S-5, give the effective boundary Hamiltonian of our model without disorder in Section S-6, show the polarizations in different iteration steps in Section S-7, discuss the SCBA in Section S-8, display the DOS in Section S-9, present more evidence on localization properties in Section S-10, and finally introduce an experimental proposal to realize and observe the HOTAI in Section S-11.

S-1. HOPPING MATRICES
We here provide the matrix form of the intra-cell hopping h 0 and the inter-cell hopping h x and h y for the Hamiltonian (1) in the main text, We have shown in the main text that the Hamiltonian respects chiral symmetry but breaks time-reversal and particlehole symmetries in the general case because h 0 is complex.

S-2. EXISTENCE OF ZERO-ENERGY CORNER MODES IN THE CLEAN CASE
To show that the Hamiltonian (1) in the main text describes a higher-order phase supporting zero-energy corner modes in the clean case with m x r =m x r = m x and m y r =m y r = m y , we write the Hamiltonian in momentum space aŝ Here where H ν (k ν , m ν ) = cos k ν σ x + (m ν + sin k ν )σ y (ν = x, y) with σ ν (ν = x, y, z) being the Pauli matrices and σ 0 being a 2 × 2 identity matrix. To see the presence of zero-energy corner modes in the system, we recast the Hamiltonian (S4) to a form in continuous real space by replacing sin k ν by −i∂ ν and cos k ν by 1 + ∂ 2 ν /2 (ν = x, y) so that H ν (k ν ) →H ν . Considering semi-infinite boundaries along x and y, if |u x and |u y are zero-energy edge modes ofH x andH y , respectively, |u x ⊗ |u y is a zero-energy mode ofH 0 localized at a corner.

S-3. EQUIVALENCE BETWEEN OUR MODEL AND THE DISORDERED BBH MODEL
We obtain a Hamiltonian respecting time-reversal, particle-hole and chiral symmetries (the BBH Hamiltonian) described byĤ if we replace h 0 with the following real matrix In fact, while the two Hamiltonians have different symmetries, they are closely related by a local transformation Specifically, one can transformĤ in Eq. (1) in the main text toĤ by the transformationĉ r → U rĉr . In other words, if Ψ Ei,rν is a spatial eigenstate of H, thenΨ Ei,rν = (−i) frν Ψ Ei,rν with f r1 = x + y − 1, f r2 = f r3 = x + y, and f r4 = x + y + 1 is an eigenstate ofH corresponding to the same energy E i . Here H andH are the first-quantization Hamiltonians of H andĤ, respectively. Therefore,Ĥ andĤ have the same energy spectrum and density profiles, indicating identical localization properties that they possess. In addition, this local phase transformation does not change the topological property, and thus the two Hamiltonians have the same topology. Under open boundary conditions, the two models are connected by the transformation irrelevant of the system size. Yet, under periodic boundary conditions, the transformation works well only when L x and L y are equal to four multiplied by positive integers. For topological property, the two models should be equivalent irrelevant of a system size given that the topology does not depend on a specific system size. For localization property, we have also calculated the IPR and LSR of the two Hamiltonians with their sizes being odd and find similar results, showing that their localization properties are irrelevant of the parity of a system size.

S-4. GENERALIZED C4 SYMMETRY
In the clean case, when m x = m y and t x = t y , the Hamiltonian (1) in the main text respects a generalized C 4 symmetry, where and g being a C 4 rotation operator such that gr = (−y, x). In momentum space, let us writeĤ = kĉ † . The generalized C 4 symmetry takes the following form where k = (k x , k y − π) and (S11)

S-5. QUANTIZATION OF QUADRUPOLE MOMENTS BY CHIRAL SYMMETRY
In this section, we will prove that the quadrupole moment is protected to be quantized by chiral symmetry and thus can be used as a topological invariant. Note that the quadrupole moment may not characterize the physical quadrupole moment, we here are only interested in the formula as a topological invariant. We consider the quadrupole moment defined by [42,43] andQ xy = nc j=1x jŷj withx j (ŷ j ) denoting the x-position (y-position) operator for electron j with n c = 2L x L y (the number of occupied bands) at half filling, and |Ψ G is the many-body ground state of a system. Let us write the many-body wave function in real space representation as where ψ n represents the nth occupied eigenstate of a first-quantization Hamiltonian. Then, the quadrupole moment can be evaluated through where Ψ G |Ψ G = ψ 1 |ψ 1 ψ 1 |ψ 2 · · · ψ 1 |ψ nc ψ 2 |ψ 1 ψ 2 |ψ 2 · · · ψ 2 |ψ nc . . . . . . . . . . . . ψ nc |ψ 1 ψ nc |ψ 2 · · · ψ nc |ψ 2 , (S15) ψ n (rν) = e i2πxy/(LxLy) ψ n (rν) and ψ m |ψ n = rα ψ * m (rα)ψ n (rα). For a generic Hamiltonian in real space H with chiral (sublattice) symmetry, U Π HU −1 Π = −H, if |Ψ n is also an eigenstate of H corresponding to energy E n , U Π |Ψ n is an eigenstate of H with energy −E n , corresponding to an unoccupied state. The set {U Π |Ψ 1 , U Π |Ψ 2 , · · · , U Π |Ψ nc } therefore constitutes the unoccupied bands. Clearly, the quadrupole moment q u xy of the unoccupied bands is equal to that of the occupied bands, since ψ m |U † Π U Π |ψ n = ψ m |ψ n given that U Π is a unitary diagonal matrix that adds a minus sign to nearest-neighbor sites. Since q xy + n with n being an integer cannot be distinguished from q xy , q xy = n characterizes a trivial phase. When all bands are occupied, the system is trivial, i.e., q xy + q u xy = n, indicating that q xy can only take either 0 or 0.5 up to an integer. The result is consistent with our numerical results where all disorder configurations exhibit quantized quadrupole moments. Similarly, one can prove the quantization of the octupole moment in 3D protected by chiral symmetry.

S-6. EFFECTIVE BOUNDARY HAMILTONIAN
In this section, we follow the transfer matrix method introduced in Ref. [46] to derive the effective boundary Hamiltonian of our system in the clean case. We will show that the effective boundary Hamiltonian at the y-normal (x-normal) edges are proportional to H x (k x , m x ) [H y (k y , m y )] up to a nonzero factor, implying that the higher-order topology can be characterized by the topological invariant P introduced in the main text.
Specifically, let us write the Hamiltonian asĤ =ĉ † Hĉ whereĉ † = (ĉ † 1ĉ † 2 · · ·ĉ † Ly ) with the index j denoting the jth layer consisting of sites along x and H reads with V n denoting the coupling between the nth and (n + 1)th layer. In disordered systems, the parameters in h n and V 2n−1 describing the intra-cell hopping are randomly generated.
In the clean case, the system has the translational invariance of period 2 and thus there are two different layers described by the Hamiltonian h 1 and h 2 , respectively. If we view these two layers as a unit cell, we use V 1 and V 2 to describe the intra-cell and inter-cell layer coupling, respectively. Now H can be simplified as (S17) Considering the periodic boundaries along x, we write h 1 , h 2 , V 1 and V 2 in momentum space as h 1 = −h 2 = H x (k x , m x ), V 1 = im y σ 0 and V 2 = σ 0 . When m y = 0, it is clear to see that the effective boundary Hamiltonian is H x (k x , m x ). When m y = 0, we obtain the following two transfer matrices at energy E where the transfer matrices connect the eigenstate in neighboring layers through with n ≥ 1 and Ψ n is the component in the nth layer of an eigenstate with the energy E. We now define the total transfer matrix at zero energy as where A(k x ) = 1 + m 2 x + m 2 y + 2m x sin k x . This matrix can be reduced to a diagonal block form through an elementary interchange transformation, where Evidently, H 1 and H 2 have the same eigenvalues. Since T is a symplectic matrix, its eigenvalues show up in pairs as (λ, 1/λ * ). Suppose that λ 1 m y (|λ 1 | > 1) is an eigenvalue of H 1 , then where U is made up of eigenvectors of H 1 and H 2 . Then, the fixed-point boundary Green's function is given by where V 0 can be chosen as any invertible matrix. By calculating eigenvectors of H 1 and H 2 , we obtain the effective boundary Hamiltonian along x is not considered as it corresponds to a phase boundary). Let us further prove that g(k x ) < 0 for all k x . Suppose m y > 0, then where we have used whereḡ(k y ) = −(m x λ 2 + 1) andf (k y ) = 1 + m 2 y + 2m y sin k y with λ 2 = − 1 2mx [B(k y ) + 1 + (B(k y ) + 1) 2 − 4m 2 x ] and B(k y ) = 1 + m 2 x + m 2 y + 2m y sin k y . Evidently, the higher-order topological phase arises when these effective boundary Hamiltonians become topological and thus can be characterized by the topological invariant P .

S-7. GRIFFITHS REGIME
In the main text, we have shown the existence of a Griffiths phase where topologically nontrivial and trivial samples coexist, leading to the topological invariant P that is not quantized. In Fig. S1, we plot the polarizations in 200 different iteration steps corresponding to different sample configurations. We see that in the Griffiths regime, some results show the polarization of 0.5 and others zero.

S-8. SELF-CONSISTENT BORN APPROXIMATION
In this section, we provide some analytical results of the SCBA to understand the disorder effects. As in the main text, we consider a disordered system by adding the following random intra-cell hopping terms at each unit cell r where Here we have changed the notation in the main text by Since we are interested in disorder without correlations, we have V i (R) = 0 and V i (R 1 )V j (R 2 ) = 1 24 δ ij δ R1R2 for i, j = 1, 2, 3, 4 with · · · denoting the average over disorder ensembles.
Based on self-consistent Born approximation, the self-energy in the presence of disorder can be calculated through the following self-consistent equation where At energy E = 0, we find numerically that the self-energy can be expanded as with Σ 0 , Σ x , Σ y being real numbers. It is clear to see that the topological masses m x and m y associated with topological properties are renormalized by disorder to new values Based on Eq. (S31), we approximate the self-energy Σ by taking Σ = 0 in the right-hand side of the equation, yielding BZ dk x dk y m y + sin(k y ) m 2 x + 2 + 2m x sin(k x ) + m 2 y + 2m y sin(k y ) .

(S36)
The correction of m x and m y due to disorder can be negative for the case of trivial quadrupole insulator with bare values m x > 1 and m y > 1 because of the positive integrands, which can induce a transition from trivial phase into topological phase with m x < 1 and m y < 1 in the presence of disorder.

S-9. THE DENSITY OF STATES
In this section, we show the DOS at E = 0 with respect to the disorder strength W in Fig. S2(a), illustrating that the DOS rises to the maximum in the multifractal phase and then fall in the trivial-II phase. Specifically, we see the development of a very narrow peak of the DOS at E = 0 in this regime [ Fig. S2(b)].

S-10. MORE EVIDENCE ON LOCALIZATION PROPERTIES
In the main text, we have demonstrated that for strong disorder, all states are localized in the gapless HOTAI and HOTI and in the Griffiths phase, all states are localized except at E = 0 where the states become multifractal. For weak disorder, all the states can be localized in the topological regime. In the trivial-I and gapped HOTI, the mobility edge appears. In this section, we provide more evidence. For clarity, we discuss the two cases with m x = m y = 1.1 and m x = m y = 0.5 separately, similar to the main text.
Specifically, we first consider the former case. Fig. S3 shows the LSR as a function of energy for five different disorder amplitudes. For small W corresponding to the trivial-I phase [ Fig. S3(a) and its inset], the LSR remains around 0.53 except at the lower band edge where it exhibits a sudden drop towards 0.386, indicating that the states at the band edge are localized. As we have already shown in the main text, the first two energy states at the energy edge near zero energy are localized. This tells us the existence of mobility edges in the trivial-I phase. In the gapped HOTAI, while we cannot conclusively determine that all states are localized when W is near the transition point, we show that this happens when W is larger. For instance, when W = 2.8, Fig. S3(b) illustrates that the LSR decreases towards 0.386 with the increase of the system size. Although the decline at E = −0.1 is not conspicuous, its decline with respect to the system size can be clearly seen in the inset. The clear fall of the normalized localization length with respect to L y is observed, further suggesting that the states at E = −0.1 are localized. These indicators show that all the states are localized. Similarly, Fig. S3(c) suggests that all the states are localized in the gapless HOTAI phase. But in the Griffiths regime, all the states are localized except at E = 0 where the LSR remains unchanged as the system size is increased [see Fig. S3(d)].
In the case with m x = m y = 0.5, Fig. S4 plots the LSR with respect to E for four different disorder strength. When the disorder is weak, e.g., W = 1, the LSR shows the existence of mobility edges in the gapped HOTI [ Fig. S4(a)]. Yet, the delocalized states can become localized in the same phase, as the disorder strength is ramped up, reflected by the decline of the LSR towards 0.386 with the increase of the system size. The localized property is also signalled by the decrease of the normalized localization length with increasing L y [ Fig. S4(b)]. Similarly, in the gapless HOTI, all the states are localized as shown in Fig. S4(c). In this case, the LSR at E = −1 decreases towards 0.386 as the system size is increased, providing further evidence for localization. In the Griffiths regime, the LSR becomes smaller for larger system sizes except at E = 0 where it remains unchanged, suggesting that the states at E = 0 are multifractal and all other states are localized. The multifractal property is also reflected by the unchanged property of the normalized localization length as L y is increased.

S-11. EXPERIMENTAL REALIZATION
In this section, we discuss in detail an experimental scheme to realize the Hamiltonian (1) in the main text using topolectrical circuits and show that the HOTAI phase can be observed by measuring two-point impedances. In fact, the quadrupole topological insulator with time-reversal symmetry without disorder has been experimentally observed in electric circuits by measuring two-point impedances [56]. One can involve disorder in the system by tuning the capacitance of capacitors and inductance of inductors to realize the HOTAI, as we have proved that this model is equivalent to our model in topological and localization properties (see Sec. 3). Here we similarly propose a scheme  to simulate the Hamiltonian (1) using topolectrical circuits and demonstrate that the HOTAI phase can be detected by two-point impedance measurements. Let us consider an electric network composed of different nodes and electric element connecting nodes, as shown in Fig. S5(a). We denote the input current and voltage of each node a by I a and V a , respectively. According to Kirchhoff's law, the circuit at a frequency of ω should satisfy the relation where R ab is the admittance of the corresponding electric element between the node a and b, and R a is the admittance of the electric element between the node a and the ground. We can rewrite the above equation into a compact form as where I and V are N -component column vectors with components I a and V a for N nodes, respectively. Here the matrix J is the circuit Laplacian. Then we can simulate our Hamiltonian H with the Laplacian J(ω) at a proper frequency ω through Each node in circuits represents one lattice site in our Hamiltonian, and each electric element linking two nodes represents the corresponding hopping between the sites, which can be either a capacitor, or an inductor, or a negative impedance converter with current inversion (INIC). For two nodes in the circuit, the electric element between them is determined according to the corresponding matrix element H ab between the site a and b in our Hamiltonian, as illustrated in Fig. S5(a). Specifically, for two neighboring sites in our Hamiltonian, if H ab is a positive (negative) real number, the electric element between a and b should be an inductor (a capacitor) with inductance (capacitance) 1 ωH ab (− H ab ω ). For the case that H ab is an imaginary number, we should connect the two nodes using an INIC with resistance 1 |H ab | and proper direction. In addition, we connect every node with the ground by appropriate electric elements to eliminate the extra diagonal terms in the Laplacian.
Similarly to the experimental work [56], we utilize the two-point impedance measurement in the circuit to characterize the zero-energy corner modes in the HOTAI phase. The two-point impedance between node a and b is defined as Z ab = n |ψ n,a − ψ n,b | 2 j n , where ψ n,a is the component for the node a of the nth eigenvector of J with eigenvalue j n . We define the impedance of each unit cell Z(x, y) as the average two-point impedance between nearest-neighbor nodes within each unit cell as Z(x, y) = [Z 12 (x, y) + Z 24 (x, y) + Z 43 (x, y) + Z 31 (x, y)]/4, where Z ij (x, y) denotes the two-point impedance between the ith node and jth node of the unit cell (x, y).