A New Numerical Method for $\mathbb{Z}_2$ Topological Insulators with Strong Disorder

We propose a new method to numerically compute the $\mathbb{Z}_2$ indices for disordered topological insulators in Kitaev's periodic table. All of the $\mathbb{Z}_2$ indices are known to be derived from the index formulae which are expressed in terms of a pair of projections introduced by Avron, Seiler, and Simon. For a given pair of projections, the corresponding index is determined by the spectrum of the difference between the two projections. This difference exhibits remarkable and useful properties, as it is compact and has a supersymmetric structure in the spectrum. These properties make it possible to numerically determine the indices of disordered topological insulators highly efficiently. The method is demonstrated for the Bernevig-Hughes-Zhang and Wilson-Dirac models whose topological phases are characterized by a $\mathbb{Z}_2$ index in two and three dimensions, respectively.

Since the discovery of Z 2 topological insulators by Kane and Mele, 1) many methods have been proposed to compute the Z 2 index. In particular, Fu and Kane found that the calculation of the Z 2 index can be considerably simplified in a system with inversion symmetry. 2) However, for disordered systems, a numerical determination of the Z 2 index is still very challenging. Roughly speaking, three types of numerical methods have been proposed so far for this problem: • The first one was proposed by Kane and Mele themselves, 1) and later some modifications were introduced. Basically, for a given model, the Z 2 index is defined by a certain Pfaffian with twisted boundary conditions. 3) The methods were applied to class AII models [4][5][6] in two and three dimensions with or without a certain inversion symmetry. [7][8][9] • The second one is based on a scattering matrix theory. 10) The Z 2 indices which are defined by the scattering matrices were numerically computed for models in the classes, AII and DIII, in two and three dimensions. 10,11) • The third one was proposed by Loring and Hastings. 12) The Z 2 indices are defined by Bott indices which are introduced as an obstruction to approximating almost commuting matrices. For models in the class AII in two and three dimensions, the Z 2 indices were numerically computed. 12,13) For systems with chiral symmetry, Loring and Schulz-Baldes 14) proposed a numerical method to obtain the values of Bott indices.
In this paper, we propose an alternative method to numerically calculate the Z 2 indices for disordered topological insulators in arbitrary dimensions. The method is based on the index formulae which were derived in Refs. [15,16]. and has the following two advantages: (i) There is no need to take an average of the Z 2 index over random variables in a model. (ii) The Z 2 index is determined by the discrete spectrum of a certain compact operator with a supersymmetric structure. The latter makes it possible to numerically determine the Z 2 index * E-mail address: akagi@cams.phys.s.u-tokyo.ac.jp highly efficiently.
Our method is demonstrated for Bernevig-Hughes-Zhang (BHZ) [17][18][19] and Wilson-Dirac 20-22) models whose topological phases are characterized by a Z 2 index of the class AII in two and three dimensions, respectively. In consequence, the method enables us to determine all of the values of the Z 2 indices of the strong and weak topological insulators, and the normal insulator phases in the phase diagrams. These values of the Z 2 indices completely coincide with the predictions in previous studies using a reliable transfer-matrix method. 18,19,22) To begin with, we introduce two Dirac operators as 15,16) D a (x) := in two dimension (2D), and in three dimensions (3D), where x = (x 1 , · · · , x d ) ∈ Z d is the position operator and a = (a 1 , · · · , a d ) ∈ R d \Z d is a vector for d = 2, 3 [see Fig. 1(a)]. The three-component vector σ is defined by σ = (σ 1 , σ 2 , σ 3 ) whose components are given by Pauli matrices σ i , each of which acts on an auxiliary Hilbert space C 2 . The whole Hilbert space is given by the tensor product of the auxiliary C 2 and the original Hilbert space for the Hamiltonian of tightbinding models which we will consider shortly. Next, we define the Z 2 index for an infinite-volume system which is a tight-binding model on a square Z 2 or cubic lattice Z 3 . 23) Let P F be the projection operator onto the states below the Fermi energy E F . The difference of two projections is defined as 16,24) A := where D * a is the adjoint of the Dirac operator D a . Then, the Z 2 index ν is defined as When the Fermi energy E F lies in a spectral gap or a mobility gap, the Z 2 index is known to be well defined. 16) In the following, we will consider such situations. Now we describe our numerical scheme for calculating the Z 2 index. Let Λ and Ω be two finite regions satisfying Ω ⊂ Λ ⊂ R d and a ∈ Ω as in Fig. 1. First, we approximate the Fermi projection P F in the infinite volume by the corresponding Fermi projection in the finite region Λ. More precisely, the approximate one is given by where |n are eigenstates of the tight-binding Hamiltonian H on Λ which we consider, and we have denoted by E n the corresponding eigenvalues. To avoid the presence of gapless edge/surface states, we impose periodic boundary conditions for the Hamiltonian in practical numerical calculations. For the operator A in Eq. (3), we replace the Fermi projection P F by P (Λ) F , and write A (Λ) for the approximate one. Further, we restrict the operator A (Λ) to the subregion Ω as where χ Ω is the characteristic function of the subregion Ω.
Under the above gap assumption, the operator A in Eq. (3) is compact even in the infinite volume limit. Therefore, the spectrum of A is discrete and has no accumulation point except for zero. This implies that an eigenstate of A is localized if the corresponding eigenvalue is not equal to zero. Let λ = 0 be an eigenvalue of A. From these observations, it is clear that the eigenvalue λ can be approximated by an eigenvalue λ ′ of the approximate operator A (Λ) Ω if the subregion Ω is sufficiently large. In addition to this, the cutoff function χ Ω in Eq. (6) enables us to remove the boundary effects due to the finite size of the region Λ.
As the first demonstration, we compute the Z 2 index of the BHZ model [17][18][19] with disorder on a square lattice Z 2 . The Hamiltonian is written as (7) where c x = [c x1↑ , c x1↓ , c x2↑ , c x2↓ ] T , and c xiα is the annihilation operator of an electron with spin α and orbital i at site x. The hopping matrices, t 1 and t 2 , are given by where τ k and s k (k = 1, 2, 3) are Pauli matrices for the orbital and the spin, respectively. Here, g 1 , g 2 and g 3 are real parameters. The on-site potential ǫ x is given by x are a random potential whose distribution is uniform in the interval [−W/2, W/2] with a positive parameter W . e 1 and e 2 are the unit vectors in the x and y directions, respectively. As is well known, [17][18][19] this Hamiltonian (7) belongs to the symmetry class AII. We set the Fermi energy E F = 0 which is located at the center of the energy gap.
In the following, we write λ i for the i-th eigenvalue of the operator A (Λ) Ω of Eq. (6) in descending order including multiplicity. 25) We note that the spectrum of A Before showing our numerical results, we abbreviate the topological and the ordinary insulator phases as TI and OI, respectively, in the phase diagram, 18,19) and write ν for the value of the Z 2 index. Figure 1(c) and (d) show, respectively, λ 1 and λ 1 − λ 2 as a function of the mass ∆ and the strength W of disorder. In the region TI, our numerical results are satisfactory because λ 1 ≃ 1 and λ 1 − λ 2 = 0. Actually, these imply ν = 1, i.e., the phase is topological as we expected. In the region OI, λ 1 is significantly different from 1, and thus ν = 0. These results show that the Z 2 index enables us to distinguish TI from OI. In OI phase, one notices λ 1 ≃ λ 2 0.2 as seen in Fig. 1(d). This degeneracy is nothing but a consequence of the two symmetries, 15,16) the time-reversal symmetry of the Hamiltonian and the supersymmetric structure of the operator A. This kind of degeneracy is very useful to determine the Z 2 index. In the region W 6, our results are in good agreement with the previous results 19) which were obtained by using a transfer-matrix method. In the region with a sufficiently large W , Anderson localization is expected to occur. For the intermediate critical region between these two regions, our method is not under control because the existence of a significantly nonvanishing spectral or mobility gap cannot be expected. In fact, our numerical results in this region show large fluctuations in both λ 1 and λ 1 − λ 2 .
The second example is the Wilson-Dirac model 22) with disorder on the cubic lattice Z 3 . The Hamiltonian is written as where Here, c x = [c x1↑ , c x1↓ , c x2↑ , c x2↓ ] T , the vector e k is the unit vector in k = x, y, z direction, and α k = τ 1 ⊗ s k , β = τ 3 ⊗ s 0 ; v x is the on-site random potential whose distribution is uniform in the interval [−W/2, W/2] with a positive parameter W . This Hamiltonian (8) belongs to the symmetry class AII for W = 0 or t 0 = 0. In the following, we will treat the case with t 0 = 0. Similarly, we abbreviate the weak, strong topological, the ordinary insulator and the diffusive metal phases as WTI, STI, OI, and M, respectively, 22) and write ν for The values of the parameters used are E F = 0, g 1 = g 2 = 1, g 3 = 0, and the system size is L 2 = 1600. The curves of the phase boundaries with the dots are plotted by using the results in Ref. [19].
the value of the Z 2 index. 26) Figure 2 shows λ 1 and λ 1 − λ 2 as a function of the mass parameter m 0 /m 2 and the strength W/m 2 of disorder. In the region OI, ν = 0 because λ 1 0.8. In the region STI, λ 1 ≃ 1 and λ 1 − λ 2 = 0, and hence ν = 1. In the region WTI, λ 1 = λ 2 ≃ 1 but λ 3 is significantly different from 1 [see Fig. 3(b)], and hence ν = 0 because the multiplicity of the eigenvalue λ ≃ 1 is two. As for the region M, we cannot expect our method to be effective because the spectral or mobility gap is expected to vanish if the metallic character of the spectrum is true. To summarize, as seen in Fig. 2, our numerical results for the Z 2 index are in good agreement with the predictions of Ref. [22]. In particular, the phase boundaries between WTI and STI, and between STI and OI are considerably sharp. Moreover, although our method is expected to be useless in the metallic phase, there do not appear fluctuations like those in the two-dimensional case. Although the Z 2 index vanishes in both the OI and WTI phases, there is a definite difference between them as follows: OI phase yields no eigenvalue λ ≃ 1 while WTI phase yields the eigenvalue λ ≃ 1 which is doubly degenerate. If the eigenvalue λ ≃ 1 is related to surface states, one can expect, from our numerical results, that the multiplicity of λ ≃ 1 coincides with the number of Dirac cones which appear on the surface of the system. In fact, it was numerically observed in a generalized Kane-Mele model on a diamond lattice that WTI (STI) phase exhibits two surface Dirac cones (odd number of Dirac cones). 27)  Ω . According to Ref. [22], STI emerges for W/m 2 7 and for m 0 /m 2 = −1.0. One can see that λ 1 is significantly different from λ 2 , and that λ 2 and λ 3 are degenerate, and similarly, λ 4 and λ 5 are degenerate. As mentioned in the case of two dimensions, this even degeneracy is a consequence of the two symmetries, 15,16) the time-reversal symmetry of the Hamiltonian and the supersymmetric structure of the operator A. For the region W/m 2 7 in Fig. 3(a), the diffusive metallic phase appears. In this region, one can see that the above double degeneracy in the spectrum of A (Λ) Ω is lifted due to the vanishing of the spectral or mobility gap in the metallic phase.
The inset of Fig. 3(a) shows system-size dependences of λ 1 and λ 2 for fixed parameters m 0 /m 2 = −1.0 and W/m 2 = 1.0. As the system size increases, λ 1 and λ 2 converge to 1 and about 0.7, respectively. Thus in the infinite-volume limit, we can clearly conclude that the Z 2 index ν is unity. One can perform a similar analysis and extract the eigenvalues of A in the infinite-volume limit for other values of the parameters.
As seen in Fig. 3(b), the double degeneracy of λ ≃ 1 appears in the region of WTI where the corresponding parameter satisfies m 0 /m 2 −2. Clearly, one can see that λ 3 is separated from λ 1 ≃ λ 2 ≃ 1, and hence the Z 2 index ν is equal to zero. Thus, if the first and the second eigenvalues, λ 1 and λ 2 , are degenerate, information about the third eigenvalue λ 3 is useful to determine the Z 2 index. Summary: We have presented our new method for numerically calculating the Z 2 indices of topological insulators with strong disorder. Our method has the following two advantages: (i) There is no need to take an average of the Z 2 index over random variables in a model. (ii) The Z 2 index is determined by the discrete spectrum of a certain compact operator with a supersymmetric structure. These properties make it possible to numerically determine the Z 2 index highly efficiently. In order to check the effectiveness of our method, we have demonstrated that all of the numerical values of the Z 2 indices completely coincide with the predictions in previous studies using a reliable transfer-matrix method 18,19,22) for the two-dimensional Bernevig-Hughes-Zhang and the threedimensional Wilson-Dirac models. Thus, the strong topological insulator phases can be characterized by the Z 2 index in the index formulae, 15,16) and can be clearly distinguished from other phases. We believe that the good agreement between the two different approaches is one of the steps toward the understanding of the nature of Z 2 topological insulators although we cannot definitely compare our method with other approaches mentioned at the beginning of the present paper. Finally, we remark that the generalization of our method to models in other symmetry classes in arbitrary dimensions is straightforward.