Application of Convolutional Neural Network to Quantum Percolation in Topological Insulators

Quantum material phases such as the Anderson insulator, diffusive metal, and Weyl/Dirac semimetal as well as topological insulators show specific wave functions both in real and Fourier spaces. These features are well captured by convolutional neural networks, and the phase diagrams have been obtained, where standard methods are not applicable. One of these examples is the cases of random lattices such as quantum percolation. Here, we study the topological insulators with random vacancies, namely, the quantum percolation in topological insulators, by analyzing the wave functions via a convolutional neural network. The vacancies in topological insulators are especially interesting since peculiar bound states are formed around the vacancies. We show that only a few percent of vacancies are required for a topological phase transition. The results are confirmed by independent calculations of localization length, density of states, and wave packet dynamics.

Introduction-Three-dimensional (3D) topological insulators (TI) [1][2][3] have been attracting considerable research attention. Oe of the many interesting features of these materials is the appearance of surface states.These states are protected topologically, namely, they appear because the topology of quantum states in the bulk of the material and that of the vacuum are different.4) A similar situation is realized when there are spherical or line vacancies in TIs.5,6) In this case, instead of surface states, bound states emerge around vacancies.The question addressed in this paper is how the bound states are connected when we have many lattice vacancies, and how the system changes with the increase in the density of vacancies.
) The main topic studied here is the change in the topological nature of the material by the quantum percolation of bound states.
Electron states on random lattice systems are difficult to study, because the conventional methods of using the transfer matrix are not applicable.The scaling analyses of the energy level statistics 13) are also difficult, if not impossible, 14) owing to a spiky density of states (DoS). 11)One way to avoid the difficulty of the transfer matrix is to introduce small but finite transfers in the disconnected bonds, which is applied to the quantum Hall effect. 15)Another approach to analyze quantum phase transitions in random lattice TI systems 16) is to calculate Chern and Z 2 numbers. 17,18) ][21][22][23][24][25][26][27][28][29] We have shown in refs. 30,31) hat a method based on a convolutional neural network (CNN) is free from the above difficulties and works well in determining the phase diagrams of quantum percolation.In contrast to the study of the DoS where all the information of the energy spectrum is obtained, our approach here is to focus on the eigenstates closest to E = 0, which makes it easier and more efficient to draw a global phase diagram.In this paper, we adopt the idea of using a CNN to draw the phase diagram of percolative 3D topological insulators (3DTIs) in site occupation probability vs gap parameter space (see Fig. 2).The results are confirmed by the calculation of the quasi-one-dimensional (Q1D) localization length via the iterative Green's function method, 32) the DoS via the kernel polynomial method (KPM), 33) as well as wave function dynamics based on the equation of motion method. 34,35) dels and Methods-We consider the following Wilson-Dirac-type tight-binding Hamiltonian on a cubic lattice, 36,37)

H =
x µ=x,y,z where |x denotes a four-component state on a site x = (x, y, z), and e µ is a unit vector in the µ-direction.α µ and β are gamma matrices, where σ µ and τ µ are Pauli matrices.m 0 is the mass parameter, and m 2 and t are hopping parameters.In the rest of this paper, we take m 2 = 1 as the energy unit and set t = 2.The parameter V x+e µ ,x is defined as Namely, V x+e µ ,x = 1 if and only if the nearest-neighbor sites x and x + e µ are connected.We assumed that a site is randomly occupied with probability P and empty with probability 1 − P, so that the Hamiltonian describes the site percolation model in 3DTIs.When all the sites are occupied, P = 1, we can analytically determine the phases: 38,39) the system is in the ordinary insulator (OI) phase for m 0 > 0, the strong topological insulator [STI(000)] phase for 0 > m 0 > −2, the weak topological insulator [WTI(111)] phase for −2 > m 0 > −4, and the STI(111) phase for −4 > m 0 > −6.The indices such as (000) and (111) are the weak indices. 40)ypical eigenfunctions at E ≃ 0 are shown in Fig. 1.We fixed m 0 = −2.5 so that the system is in the WTI(111) phase in the clean limit, P = 1.The fixed boundary condition (FBC) is imposed in the z direction, while periodic boundary conditions (PBCs) are imposed in x and y directions.For P = 0.99, the gapless surface state appears in the x − y plane [Fig.1(a)] as in the clean limit.When we reduced the site occupancy to P = 0.97, the surface states disappear, and the wave function is extended over the whole system [Fig.1(b)].Note that the surface state is destroyed by only 3% site vacancies despite the robustness of topological materials against disorder.When we further reduced the site occupancy to P = 0.90, the surface state reappears [Fig.1(c)].We first derive a rough phase diagram of quantum percolation in 3DTI.We adopt supervised learning with the CNN, which is an efficient machine learning approach to draw phase diagrams of random electron systems. 19,20,30,31) Inhe supervised training, we need a correctly labeled data set (training data) in advance to optimize the weight parameters of the CNN.Without prior knowledge of quantum percolation in 3DTI, however, preparing enough training data for each phase is difficult.We therefore use the weight parameters of the CNN that has learned the features of k-space eigenfunctions of the 3DTI with random on-site potential, 31) where the lattice is regular (V x+e µ ,x ≡ 1), but the random on-site potential is added with V(x) random potentials.Because of the generalization capability of the CNN, we expect it to determine the correct phases for the unlabeled data set (the wave functions of the 3DTI with vacancies). 30,31) n analogy of this approach is that the CNN trained for the Anderson transition with random on-site potential can correctly identify quantum phases in quantum percolation.30) To train the CNN, we need to prepare the k-space eigenfunctions with varying m 0 and random potential.20,31) We consider the 24×24×24 lattice with the FBC in the z direction and the PBC in the other directions, diagonalize the Hamiltonian in real space, calculate the eigenfunctions around E = 0 using the sparse matrix diagonalization Intel MKL/FEAST, 41) and obtain the k-space eigenfunctions through discrete Fourier transform.Random numbers are generated by the Mersenne Twister algorithm.42) We also study the Q1D localization length for the detailed analysis. As ntioned in the introduction, in the quantum percolation, the transfer matrix method, an effective method for the case of a regular lattice, is not applicable since a matrix relating one layer to the other has zero determinant owing to disconnected bonds.31) In this paper, we employ the itera-tive Green's function method 32) and calculate the Q1D localization length in the geometry L × L × L z , L z ≫ L with PBCs in the transverse (x and y) directions.(A similar geometry is employed in the case of wave packet dynamics simulation.)For numerical simulation, the system size in the z direction is truncated at L z = 100, 000.We shifted the energy slightly from the band center, E = 0 to E = 0.001, to avoid the numerical instability of the inverse matrix calculation. Note tht if all the clusters that include sites of the first layer are disconnected at z < L z , the Green's function method also breaks down.In practice, this situation does not occur because the parameter regime studied here is well above the classical percolation threshold P c ≈ 0.312.8,43) In the case of DoS calculation, cubic systems of size L = 160 with PBCs are considered.
Results-In Fig. 2, we show the phase diagram of quantum percolation in the 3DTI obtained by the CNN.The abscissa shows the mass parameter m 0 , while the ordinate is the site occupancy P. From this figure, we see that as the site occupancy P decreases, namely, the disorder increases, the absolute value of the topological mass |m 0 + 3| effectively increases by the renormalization of m 0 .The STI(000) and the STI(111) phases are therefore in contact with each other around m 0 = −3, and the diffusive metal (DM) phase appears around there since topologically different phases cannot be connected continuously.In the case of random on-site potential, in contrast, the renormalized mass |m 0 + 3| decreases with the increase in disorder, resulting in the transitions from the OI phase to the STI phase and from STI to WTI (so-called topological Anderson insulator transition [44][45][46][47] ).The difference between our case and that of the random on-site potential indicates that the role of lattice vacancy disorder is qualitatively different.We also emphasize that the transition from the WTI(111) phase to the STI phase occurs with less than 5% lattice vacancies, despite the common belief that the topological phase is robust against randomness.The CNN outputs the confidence values P OI , P W111 , P S000 , P S111 , and P DM that the wave function belongs to OI, WTI(111), STI(000), STI(111), and DM, respectively, and the intensity 0 × P OI + 1 × P W111 + 2 × P S000 + 3 × P S111 + 4 × P DM is plotted.Sample-to-sample fluctuations are confirmed to be small, so the sample average is not taken.
Since the above phase diagram shows unexpected features, we verify the phase boundaries by calculating the Q1D localization length ξ Q1D via the iterative Green's function method. 32)Figure 3 shows the normalized localization length [48][49][50] Λ = ξ Q1D /L as a function of P for m 0 = −1.5, −2.2, and −2.4. W see that for each m 0 , the peak of the localization length is located at the phase boundary in Fig. 2. At m 0 = −2.4,for example, the localization length exhibits two peaks around P = 0.97 and P = 0.56, corresponding to the phase boundary of WTI(111)-STI(000) and that of STI(000)-OI, respectively (Fig. 2).This increase in Q1D localization length at the phase boundary is consistent with the well-known fact that Dirac semimetals (DSMs) continue to exist at the topological phase boundary even in disordered systems.51,52) To further understand the nature of the topological phase transition due to vacancies, we also study the DoS using the KPM.33) The DoSs in the clean limit, P = 1, at m 0 = −1.5, −2.0, and −2.5 are shown in Figs.4(a)-4(c), respectively.The case m 0 = −2.0(b) corresponds to the phase boundary between the STI(000) and the WTI(111) phases, on which the system is a DSM.The linear energy dispersion around E = 0 leads to a parabolic DoS, as is clearly seen in Fig. 4(b).
The DoSs with 4% lattice vacancies, P = 0.96, with the same m 0 's are shown in Figs.4(d)-4(f).When lattice vacancies are present, the bound states appear around them, and these bound states form minibands inside the original band gap.In Fig. 4(f), the miniband shows a parabolic DoS in the vicinity of E = 0 (see the inset), which is consistent with the phase diagram: the parabolic DoS appears for the parameter set (m 0 , P) = (−2.5,0.96), which is on the phase boundary in Fig. 2. In the following, we show that this parabolic DoS in the miniband comes from DSM states on the phase boundary.
A DSM is characterized by its ballistic transport.To confirm this, we employ the Chebyshev polynomial expansion for the time-evolution operator 34) U(∆t) = exp(−iH∆t), and study the wave packet dynamics in the Q1D geometry.The results for two parameter sets corresponding to the phase boundaries, (a) (m 0 , P) = (−2.0,1.0) and (b) (m 0 , P) = (−2.3,0.975), are shown in Fig. 5.The ballistic transport in the absence of disorder [Fig.5(a)] survives even in the presence of randomness [Fig.5(b)].The decay of the wave packet in Fig. 5(b) is because part of the wave packet is trapped by the bound states around the vacancies during the transport process. 53,54) rom the time dependence of the peak positions, we evaluate the speed v of the ballistic transport and obtain (a) v a = ta/ = 2 × m 2 a/ (a being the lattice constant) and (b) v b = 1.15 × m 2 a/ .The ratio α of the renormalized velocity to the bare one, α = v b /v a , is estimated to be 0.575.We also estimate the curvatures of the DoSs around E = 0 by fitting the DoS with a quadratic polynomial in the vicinity of E = 0, and obtain (a) ρ a ≈ 0.0103 E 2 and (b) ρ b ≈ 0.0536 E 2 .From ρ ∝ E 2 /v 3 , we estimate α ′ = (ρ a /ρ b ) 1/3 = 0.577, in good agreement with the estimate of the ballistic velocity ratio α.We therefore conclude that the parabolic DoSs inside the minibands [e.g., Fig. 4(f)] are formed by a DSM.
Summary and concluding remarks-In this paper, we have shown that a small amount of vacancies induces a topological phase transition.For example, a weak topological insulator undergoes topological phase transition and becomes a strong topological insulator with only a few percent of vacancies.A few percent of vacancies, at first sight, may sound very small.Around a vacancy, however, a bound state on the order of 10 sites is formed.Such bound states therefore fill a significant amount of sites even for a few percent of vacancies, resulting in metallic states that spread over the system [see Fig. 1(b)].We emphasize that the vacancy-driven topological transition is in sharp contrast to the case of on-site random potential, where the topological phase transition does not occur as long as the bulk band gap remains.Note also that contrary to the present problem, in the case of on-site potential, the strong-toweak topological insulator and the ordinary-to-strong topological insulator transitions take place with the increase in randomness.
][57] Furthermore, we have an additional sublattice symmetry for m 0 = −3m 2 , where the last term in Eq. ( 1) vanishes.The latter symmetry is rather artificial, since this term may fluctuate and deviate from 0 in actual materials, and some features of the phase diagram such as symmetry around m 0 = −3m 2 will not be observed in real materials.The sensitivity of topological phases against vacancies, however, remains even in the presence of site-to-site fluctuation of 3m 2 .We have confirmed that changing the last term in Eq. ( 1), (m 0 + 3m 2 ) x |x β x|, to x (m 0 + m(x)/2)|x β x| with m(x)/m 2 given by the number of bonds connected to the site x breaks the symmetry of the phase diagram around m 0 + 3m 2 , but the overall features of the phase diagram remain unchanged.The addition of on-site randomness, Eq. ( 4), with V(x) independently and uniformly distributed in [−W/2, W/2], derives the system to the Wigner-Dyson symplectic class AII, [56][57][58][59][60] but as long as W is small, say W = 0.5, the features of the phase boundaries remain almost unchanged.

Fig. 3 .
Fig. 3. (Color) Normalized localization length Λ as a function of P for m 0 = −1.5, −2.2, and −2.4.The size of the cross section is L = 8.The peak of the localization length indicates the phase boundary where the states are (semi)metallic and show a larger localization length.We slightly shift the energy to E = 0.001 from E = 0. We have confirmed that the peak position is insensitive to the small shift by examining the case for E = 0.0001.