Disorder-Induced Quantum Phase Transitions in Three-Dimensional Second-Order Topological Insulators

Disorder effects on three-dimensional second-order topological insulators (3DSOTIs) are investigated numerically and analytically. The study is based on a tight-binding Hamiltonian for non-interacting electrons on a cubic lattice with a reflection symmetry that supports a 3DSOTI in the absence of disorder. Interestingly, unlike the disorder effects on a topological trivial system that can only be either a diffusive metal (DM) or an Anderson insulator (AI), disorders can sequentially induce four phases of 3DSOTIs, three-dimensional first-order topological insulators (3DFOTIs), DMs and AIs. At a weak disorder when the on-site random potential of strength $W$ is below a low critical value $W_{c1}$ at which the gap of surface states closes while the bulk sates are still gapped, the system is a disordered 3DSOTI characterized by a constant density of states and a quantized integer conductance of $e^2/h$ through its chiral hinge states. The gap of the bulk states closes at a higher critical disorder $W_{c2}$, and the system is a disordered 3DFOTI in the lower intermediate disorder between $W_{c1}$ and $W_{c2}$ in which electron conduction is through the topological surface states. The system becomes a DM in a higher intermediate disorder between $W_{c2}$ and $W_{c3}$ above which the states at the Fermi level are localized. It undergoes a normal three-dimension metal-to-insulator transition at $W_{c3}$ and becomes the conventional Anderson insulator for $W>W_{c3}$. Our results can be understood from the self-consistent Born approximation that allows one to reveal how the original Hamiltonian been modified by the on-site disorders.

All those newly discovered topological phases should survive in disorder since they are genuine states of matter although disorder has profound effects on them as demonstrated in the Anderson topological insulators [28][29][30].Currently, the second-order topological insulators, characterized by topological states in (d − 2)-dimensional boundary, are the focus because of its experimental realizations in phononics [9][10][11], photonics [14], and circuitry [31].So far, most of the works are on the constructions of second-order topological insula-tors in crystal with well-defined crystalline symmetries, the fate of such a phase under disorder has been less studied [32][33][34].Thus it should be very interesting to find out how the hinge states in 3DSOTIs are modified by disorders.From the knowledge of Anderson localization for topologically trivial states, it is known that the competition between the variation of random energy and the bandwidth determines the metal-toinsulator transition.For 3DSOTIs, the gap ∆ 1 of bulk states and gap ∆ 2 < ∆ 1 of the surface states should be important because a transition from a topological state in a submanifold of dimension (d − n + 1) to another in its boundary of dimension (d − n) can only happen when the band gap in the (d − n + 1) manifold closes at the transition point.We expect a 3DSOTI undergoing three phase transitions that involve four different phases as disorder strength W increases: The 3DSOTI remains stable up to a critical disorder W c1 at which the gap ∆ 2 of surface states closes (∆ 2 = 0) while bulk gap remains open (∆ 1 0).The system enters into the 3DFOTI from the 3DSOTI.Further increase of disorder to the second critical value of W c2 at which bulk gap closes, the 3DFOTI is replaced by the conventional topologically trivial DM.The third quantum phase transition is expected to occur at a strong critical disorder W c3 at which all states are localized and the system becomes an AI.The route of these expected quantum phase transitions is schematically illustrated in Fig. 1.
In this work, we use a 3DSOTI model of class A with a reflection-symmetry in the ten Altland-Zirnbauer classification [35] to verify the physics in Fig. 1, as well as a more comprehensive phase diagram explicated in Supplemental Materials [36].Using highly-accurate numerical calculations, we show that above three quantum phase transitions occur indeed when disorder strength W increases.The 3DSOTI is featured where c † i (c i ) are the electron creation (annihilation) operators at site i = (n x , n y , n z ).M is the Dirac mass that controls the band inversion, and B is a parameter for controlling gap opening on surfaces.Γ 0 and Γ µ=1,2,3,4,5 are, respectively, the fourby-four identity matrix and the five non-unique Dirac matrices satisfying {Γ µ , Γ ν } = 2δ µ,ν Γ 0 and Γ µν = [Γ µ , Γ ν ]/(2i).Here we choose Γ (1,2,3,4,5) with the Pauli matrices s µ and σ µ acting on spin and orbital spaces, respectively.Hopping energy t = 1 is chosen as the energy unit.v i is a white noise, distributing uniformly in the range of [−W/2, W/2].
Clean case.−In the absence of disorder (W = 0), Hamiltonian (1) was well studied [20,42] and can be block diagonalized in the momentum space, Here For B = 0 and 1 < M < 3, Eq. ( 2) describes a reflection-symmetric strong 3DFOTI with reflection plane on x = 0 [20,37].The Hamiltonian does not change under the reflection symmetry of Γ 54 = s 0 ⊗ σ 1 , According to the bulk-boundary correspondence, the nontrivial bulk topology guarantees the appearance of the gapless surface states on the self-reflected surfaces, e.g., 1 on z = 0, while the gapped surface states on the non-reflected surfaces [38][39][40][41].If two non-reflection surfaces encounter at the reflection plane under a sharp angle, the last term BΓ 31 leads to the band inversion of surface states, as well as the emergence of hinge states at their boundary [20,42].Thus, Hamiltonian (1) supports a 3DSOTI.To visualize the above picture, we plot the energy spectrum E(k 3 ) of a rectangle shape with open boundary conditions (OBCs) on the (110) and (1 10) surfaces and periodic boundary condition (PBC) on the z direction for M = 2, B = 0 (Fig. 2(a)) and M = 2, B = 0.2 (Fig. 2(b)).Clearly, for B = 0, topological surface states exist in the bulk gap, i.e., E ∈ [−∆ 1 /2, ∆ 1 /2], while for B = 0.2 surface states are gapped in E ∈ [−∆ 2 /2, ∆ 2 /2] and the hinge states emerge in the gap, exactly the same as reported results [20,42].This can also be seen from the distribution of wave function of E = 0 shown in Fig. 2(c,d).
. The OBCs on the (110), (1 10), and (001) surfaces are assumed if not specified otherwise.Notice that hinge states in 3DSO-TIs are one-dimensional helical edge channels that give rise to quantized conductances and zero conductance fluctuations [20] if no other conduction channels exist.This feature can be used to test the robustness of 3DSOTI hinge state against weak disorder.Thus, we compute the two-terminal conductances g L of a disordered bar connected by two semi-infinite leads along z direction by the Landauer-Bttiker formula [43].The Fermi energy is fixed at E = 0.02 to focus on the hinge states.Figures 3(a) and (b) show, respectively, the sample-averaged dimensionless conductance g L and the conductance fluctuation δg L = ( g 2 L − g L 2 ) 1/2 as a function of disorder W for various sizes from L = 20 to L = 32.Clearly, all g L are exactly quantized at 1 for W < W c1 2.2 with zero conduc-tance fluctuation, a typical feature of hinge states.Beyond W c1 , g L is not quantized and δg L 0. g L increases with L, an indication of states of E = 0.02 being extended.As demonstrated later, these states are surface states, and the system is a 3DFOTI.The dispersion relation of the hinge states in clean 3DSO-TIs is linear in k 3 , i.e., E hinge = ±ck 3 , see Fig. 2(b).Thus its contribution to the density of states (DOS) is a constant, i.e., ρ(E) = 1/c.Interestingly, the average DOS, defined as ρ(E) = 1/(4L 3 ) q 4 p=1 δ(E − E p,q ) , of disordered 3DSOTI from the kernel polynomial method [45] is not only a constant, but also independent of disorders.Averaged ρ(E) of the disordered bar of L = 66 for various W that are obtained from at least 100 ensembles with 1024 Chebyshev moments is plotted in Fig. 3(c).Apparently, the width of the constant DOS regime becomes smaller as W increases.This is expected since disorders tend to reduce the gap ∆ 2 .For large enough disorders W > W c1 , the constant plateau of ρ(E) disappears when the gap ∆ 2 of surface state vanishes and the system becomes a 3DFOTI.Therefore, constant DOS can be a fingerprint of the disordered 3DSOTI.
3DFOTI-to-DM.−Havingestablished the existence of 3DSOTIs for W < W c1 , we now show that the system is a 3DFOTI for W ∈ [W c1 , W c2 ], where the zero-energy states are the surface states rather than the hinge states, and becomes a DM for W > W c2 .Note that the E = 0 wave functions of both 3DFOTIs and 3DSOTIs are highly localized at system boundaries, either on the surfaces in a 3DFOTI or at hinges in a 3DSOTI, in contrast to extending over the whole systems for DMs.Therefore, we introduce the following quantity to distinguish 3DSOTIs and 3DFOTIs from DMs: with the first summation taking over all the sites of the surface.
It is natural to expect that ζ W,L=∞ is a finite non-zero constant for 3DSOTIs and approaches zero for a DM since the ratio of number of surface sites to that of bulk size goes to zero.For a fixed W, ζ L,W should increase with L for 3DFOTIs and 3DSO-TIs and decreases with L for DM.They intersect at critical disorder W c2 .To substantiate the above assertion, we evaluate the sampleaveraged ζ W,L for various L's ranging from 20 to 50 and various W [46][47][48].The results are plotted in Fig. 4(a).A phase transition, between the boundary states (dζ W,L /dL > 0) and the bulk states (dζ W,L /dL < 0) at W = W c2 4 at which all ζ W,L curves cross, can be clearly seen.Since W c2 > W c1 , 3DFOTI, where wave function of E = 0 is localized on the surface rather than the edge of the bar as illustrated in Supplemental Materials [36], exists between the 3DSOTI at weak disorders W < W c1 and the DM beyond W c2 (W > W c2 ).We also consider Fig. 4(a) as an empirical verification of the existence of 3DFOTI-to-DM quantum phase transitions such that E = 0 states for W c1 < W < W c2 and W > W c2 belong, respectively, to the 3DFOTI and the DMs in the thermodynamic limit of L → ∞.
We expect that the 3DFOTI-DM transition happens when the bulk gap closes.To substantiate it, we calculate ρ bulk (E) for the disordered bar of L = 66 with PBCs along all direction so that all boundary states (surfaces or hinges) are eliminated, and only bulk states can contribute to DOS.The DOS for various disorders are displayed in Fig. 4(b).As expected, ρ bulk (0) = 0 below W c2 and ρ bulk (0) 0 beyond the critical disorder of W c2 4.0, in contrast to non-zero ρ(E) around E = 0 in Fig. 3(c) for W < W c2 .This demonstrates from a different angle that non-zero ρ(E) around E = 0 in Fig. 3(c) is from either hinge states of 3DSOTI or surface states of 3DFOTI.The estimate of the critical disorder strength is consistent with that by finite-size scaling of ζ W,L .
We carry out the self-consistent Born approximation (SCBA) to further understand the disorder affects.The self energy Σ [29,30,49,50] is where N is the total number of sites.For simplicity, the BΓ 31 term (B M in this work) is neglected and Σ can be expressed as Σ = Σ 0 Γ 0 + 4 µ=1 Σ µ Γ µ .For E = 0, Σ 0 = −i(1/τ) is a pure imaginary number, with τ being the lifetime of the zero-energy bulk states, i.e., ρ bulk (E = 0) ∝ (1/τ).After some calculations (see Supplemental Materials [36] for more details), we obtain with Σ 1,3,4 = 0.The solutions of Eq. ( 6) can be either 1/τ = 0 for W ≤ W c2 or 1/τ 0 for W > W c2 .The former corresponds to either 3DFOTIs or 3DSOTIs where ∆ 1 0 and ρ bulk (0) = 0, while the later is for DMs with ∆ 1 = 0 and ρ bulk (0) 0. The critical disorder W c2 is given by the gap equation [50] Numerically, we obtain W c2 3.7 for M = 2 and B = 0.2, consistent with the estimates from ζ W,L and ρ bulk (E).According to the SCBA, M is renormalized by the disorder as M = M + ∆ with such that the phase boundary between 3DSOTIs, 3DFOTIs and topological trivial phase are shifted by disorders, see Supplemental Materials [36].
DM-to-AI.−Theextended states in DMs are eventually localized by strong disorders.To investigate the nature of this Anderson transition and its associated universality class, we compute the PR p 2 (E = 0), which measures how many lattice sites are occupied by the wave function of E = 0 [51][52][53].Near the critical disorder W c3 of the Anderson transitions, p 2 satisfies the one-parameter scaling function [54,55]   Near W c3 , the calculated ln Y L (W) and ln f (x) are displayed in Figs.4(c) and (d), respectively.Indeed, data in Fig. 4(c) give W c3 = 18.73 ± 0.03, and dY L (W)/dL is always positive (negative) for W < W c3 (W > W c3 ), indicating the system is a DM (AI).Following the well-established scenario [55], we find D = 1.7±0.2, and ν = 1.5±0.2(see details in Supplemental Materials [36]).The obtained ν and D, characterizing the universality class of transitions, are consistent with previous estimations for 3D Gaussian unitary ensemble [56].
Phase diagram.−Amore inclusive picture is obtained by exhaustive calculations of different M ranging from M = 1.7 to 2.3 and fixed B = 0.2.A phase diagram in the W − M plane elucidated the occurrence of the 3DSOTI, 3DFOTI, DM, and AI phases is shown in Fig. 5. Seemingly, W c1 and W c2 respectively monotonically decreases and increases with M while W c3 does not depend on M.Although the phase boundaries, e.g., W c1 , W c2 , and W c3 , depend on the details of a model, Fig. 5 is a strong evidence of the generality to the physics shown in Fig. 1.
Conclusion.−In conclusion, a generic route of disorderinduced phase transitions for a crystal 3DSOTI is revealed.As random potential strength increase, the 3DSOTI transforms to a 3DFOTI at a lower weak critical disorder of W c1 , followed by a second transition at a higher critical disorder of W c2 to a DM.The system eventually becomes an AI after a metal-to-insulator transition at critical disorder of W c3 .The 3DSOTI is featured by quantized conductance exactly at e 2 /h and zero conductance fluctuation, as well as the constant density of states, while 3DFOTIs are identified by their dominate occupation on surfaces, negligible occupation in the bulk and on hinges.The DM and AI are confirmed by the scaling analysis of participation ratios, and the corresponding Anderson localization transition is found to belong the conventional 3D Gaussian unitary class.
This work is supported by the National Natural Science Foundation of China (Grants No. 11774296, 11704061 and 11974296) and Hong Kong RGC (Grants No. 16301518 and

FIG. 3 .
FIG. 3. (a,b) g L (a) and δg L (b) as a function of W for E = 0.02 and various L. The dash lines locate W c1 .(c) ρ(E) for E ∈ [−0.15, 0.15], L = 66 and various W. The solid (dotted) lines are for W ≤ W c1 (W > W c1 ).The dash line guides the eyes for a non-zero constant ρ(E) = 1/c.
where f (x) is the unknown scaling function, C and y > 0 are a constant and the exponent of the irrelevant variable, respectively.The correlation length ξ diverges at W c3 as ξ ∝ |W − W c3 | −ν with critical exponent ν.D is the fractal dimension of critical wave functions which occupy a subspace of dimensionality smaller than the embedded space dimension d = 3.By defining Y L (W) = p 2 L −D − CL −y , we use the following criteria to identify the Anderson transitions[55]:(1)

FIG. 5 .
FIG. 5. Phase diagram in the W − M plane displaying the occurrences of 3DSOTIs, 3DFOTIs, DMs, and AIs.Colors map g L=12 (W, M) .The boundary between 3DSOTIs and 3DFOTIs (the red dash line) is determined by g L=12 (W, M) = 0.999 while that between FOTIs and DMs is determined by ρ bulk (0).The boundary between DMs and AIs is given by scaling analysis of p 2 .