PHYSICAL REVIEW RESEARCH 2, 013229 (2020) Computation of topological phase diagram of disordered Pb1−x Snx Te using the kernel polynomial method Dániel Varjas,1,2,* Michel Fruchart,3,4 Anton R. Akhmerov,2 and Pablo M. Perez-Piskunow 2,5,† 1 QuTech, Delft University of Technology, P.O. Box 4056, 2600 GA Delft, The Netherlands 2 Kavli Institute of Nanoscience, Delft University of Technology, P.O. Box 4056, 2600 GA Delft, The Netherlands 3 Instituut-Lorentz, Universiteit Leiden, P.O. Box 9506, 2300 RA Leiden, The Netherlands 4 James Franck Institute and Department of Physics, University of Chicago, Chicago, Illinois 60637, USA 5 Catalan Institute of Nanoscience and Nanotechnology (ICN2), CSIC and BIST, Campus UAB, Bellaterra, 08193 Barcelona, Spain (Received 5 September 2019; accepted 6 February 2020; published 28 February 2020) We present an algorithm to determine topological invariants of inhomogeneous systems, such as alloys, disordered crystals, or amorphous systems. Based on the kernel polynomial method, our algorithm allows us to study samples with more than 107 degrees of freedom. Our method enables the study of large complex compounds, where disorder is inherent to the system. We use it to analyze Pb1−x Snx Te and tighten the critical concentration for the phase transition. Moreover, we obtain the topological phase diagram for related alloys in the family of three-dimensional mirror Chern insulators. DOI: 10.1103/PhysRevResearch.2.013229 I. INTRODUCTION cost with ξ is the main distinction between different numerical approaches. Topological materials have attracted continuing interest Available methods to compute topological invariants either from both the fundamental physics and the material science apply the momentum space Berry curvature formalism to communities for the last decade [1,2]. The program to the- periodic systems with a disordered supercell or use a real- oretically classify noninteracting crystalline insulators has space formulation on a large finite sample [19–27]. However, been completed, tabulating possible topological phases in all these methods involve solving at least one eigenvalue equation space groups [3–5]. Recent efforts focus on automated high- with size equal to the number of degrees of freedom, resulting throughput methods to discover and classify new topological in the complexity ξ 3d in d dimensions. This restricts the materials [6–8], culminating in the production of compre- applicability of such methods to small system sizes, especially hensive databases of topological insulators and semimetals in three dimensions. To our knowledge, the most efficient [8–12]. method in three dimensions is the scattering matrix approach However, not all topological insulators are compounds [28], with a complexity scaling as ξ 3(d−1) , allowing for maxi- with perfect stoichiometry. The first three-dimensional (3D) mum sizes of 5 × 105 degrees of freedom. topological insulator to be predicted [13] and experimentally We present an algorithm to efficiently identify topologi- realized [14] was an alloy—Bix Sb1−x . These systems are cal phases of strongly disordered systems using the kernel usually studied using the virtual crystal approximation (VCA) polynomial method (KPM) [29–31], an approximation based or the coherent potential approximation, which approximate on a polynomial expansion of the quantities of interest. All an alloy by a perfect crystal [15,16]. This approach ignores topological properties of a noninteracting system of electrons the intrinsic disorder in alloys, and it is insufficient to explain are encoded in the projector on the occupied states (spectral topological transitions that appear at strong disorder [17,18], projector), which is efficiently approximated using KPM. Our or accurately find critical concentrations. algorithm builds on the method of topological markers [22] to The topological invariant converges to its bulk value in construct a topological invariant as the trace of an operator. samples larger than the localization length ξ . This is the main Since topological invariants are integers, it is sufficient to limitation in resolving topological phase transitions as ξ di- reduce the statistical uncertainty below 1/2 to obtain the exact verges. Therefore, the asymptotic scaling of the computational value. In particular, the stochastic evaluation of traces [29] from a small number of random vectors, combined with KPM, is suited for this task. With ξ d+1 scaling of the computational * dvarjas@gmail.com effort, it is the most efficient method in three dimensions. † pablo.perez.piskunow@gmail.com As a concrete example, we apply our method to lead-tin- telluride Pb1−x Snx Te alloys—three-dimensional topological Published by the American Physical Society under the terms of the crystalline insulators characterized by a mirror Chern number. Creative Commons Attribution 4.0 International license. Further Thanks to the efficiency of our algorithm, we are able to ana- distribution of this work must maintain attribution to the author(s) lyze the topological properties of large systems that are hardly and the published article’s title, journal citation, and DOI. accessible with previously existing techniques. We study 3D 2643-1564/2020/2(1)/013229(11) 013229-1 Published by the American Physical Society DÁNIEL VARJAS et al. PHYSICAL REVIEW RESEARCH 2, 013229 (2020) the number of degrees of freedom N, so the computational cost of such methods is order ξ 3d , restricting them to small system sizes in three dimensions. The scattering invariant formalism [44] avoids full diago- nalization and only requires the knowledge of the scattering matrix at the Fermi level. The most efficient known algorithm for computing the scattering matrix is based on the nested dissection method [45] and scales as ξ 3d−3 for d > 1. An alternative way to obtain the Chern number in two- dimensional (2D) systems is to compute the Hall conductivity. An efficient method to compute the Kubo-Bastin conductivity is given by García et al. [30] for large disordered systems. This approach uses a two-dimensional KPM expansion to construct a matrix kernel that holds the information of the Hall conductivity for all energies, which is integrated to compute the total Hall conductivity for any Fermi energy. The time and memory requirements of constructing the matrix kernel are large, of order NM 2 for a system with N orbitals and a KPM expansion of order M. Using our method, however, we FIG. 1. Top panels: Surface spectra of a 20 × 80 × 80 sample of compute the Chern number for a single value of the Fermi Pb1−x Snx Te in the trivial (a) and topological (b) phase. The presence energy with time scaling NM and memory requirement of of a gapless surface Dirac cone indicates the mirror Chern insulator order N. phase. (c) Transition between trivial and mirror Chern phase when varying x for Pb1−x Snx Te calculated using our method with various system sizes. Inset: Finite-size collapse of the curves with xc  III. GENERAL STRATEGY 0.28(3) and ν  0.9(6). Computing the exact spectral projector  systems with linear sizes of over 100 lattice constants (L > P̂ = θ (EF − Ĥ ) = |n n| (1) n:En ξ . To ensure that the V. TIGHT-BINDING MODEL OF Pb1−x Snx Te invariant ν obtained with the trace of the local marker con- verges to the momentum space topological invariant we must Topological crystalline insulators protected by reflection repeat the supercell with the disorder realization over the symmetry [36] were theoretically predicted [51] and exper- whole space. Since the approximation of the local operator imentally observed [34,37–41] in Pb1−x Snx Te alloys. They ν̂ is localized with ξ < L, it is sufficient to repeat two times host metallic surface states on the surfaces that are symmetric the supercell in every spatial direction, and compute the with respect to the mirror plane [2,51,53]. Lead-tin-telluride topological invariant as an average of the topological marker was studied using either the VCA [32,34] or ab initio methods over the central L d volume. For a detailed description, see [33], finding a gap closing and phase transition near x = 0.35 Appendix C. or 0.23, respectively. We use a tight-binding approach that The resulting complexity of the computation depends lin- captures long-range correlations, and find a different critical early on the number of random vectors R used (typically of concentration. order 1), on the number of moments M, and on the number We consider the substitutional disorder of the lead-tin- of sites of the system L d . We use a sparse representation telluride alloy Pb1−x Snx Te resulting from substituting some of the short-ranged Hamiltonian. As a result the memory Pb ions for Sn ions. This disorder is nonmagnetic, and it requirement scales linearly with the system size L d and is preserves the reflection symmetry on average, which is suf- independent of other parameters. Setting all quantities to their ficient to protect the gapless surface states [28,54]. We dis- minimal values (L ∼ ξ and M ∼ ξ ) results in an algorithm regard other types of symmetry-breaking disorder appearing with a computational cost scaling of ξ d+1 . naturally in Pb1−x Snx Te [54], such as ferroelectric structural distortion [51,54,55] and magnetic dopants [51,54,56]. IV. APPLICATION TO MIRROR CHERN NUMBER In our investigation we use two atomistic tight-binding models. The first one includes 18 spinful s, p, and d orbitals on Our method provides better scaling than existing ap- both sublattices, with 36 bands in total. This model accurately proaches in d  3. We apply it to disordered 3D mirror Chern describes the energetics, using tight-binding parameters for insulators. These are a widely studied class of topological both SnTe and PbTe derived from ab initio simulations [32]. crystalline materials with a Z topological classification that To simulate the alloy, we substitute randomly Sn for Te relies on reflection symmetry [51]. Several experimental real- with probability x. We incorporate substitutional disorder by izations are known, including alloys [34,37–41]. using the hopping amplitudes of SnTe for Sn-Te bonds and In a reflection-symmetric system of fermions, all wave PbTe amplitudes for Pb-Te bonds. The onsite parameters functions are eigenstates of the mirror operator M̂z , with of Te atoms are slightly different in SnTe and PbTe; we eigenvalues ±i. The Chern numbers C± for mirror-even and therefore use a weighted average of these depending on the mirror-odd wave functions are local environment. We also use the appropriate onsite terms, C± = 2π i Tr A [x̃± , ỹ± ], (5) including L · S spin-orbit coupling (SOC), depending on the type of the Sn or Pb atom. Further details can be found in where x̃± = M̂± P̂x̂ P̂M̂± and ỹ± = M̂± P̂ŷP̂M̂± are the pro- Appendix D. jected position operators restricted to the mirror-even or When investigating the onsite energy dependent phase mirror-odd subspaces. Here M̂± are the projectors on the diagram of X1−x Snx Te alloys, we use a simplified model that 013229-3 DÁNIEL VARJAS et al. PHYSICAL REVIEW RESEARCH 2, 013229 (2020) FIG. 2. Comparison of the critical concentration obtained with FIG. 3. Mirror Chern number CM (a) and localization length ξ our method, against theoretical approximations (squares) and exper- (b) of the X1−x Snx Te model as a function of onsite energy mX and imental measurements (circles). The shapes are empty for the trivial composition x. In both panels the overlay shows the phase boundary phase, half filled for the estimated transition point (transition range in the VCA and the horizontal dashed line corresponds to X = Pb. for Zhong et al. [41]), and fully filled for the topological phase. ξ was calculated from the scaling of conductivity with sample sizes up to L = 20 with ten disorder realizations. CM was calculated in a system of 40 × 40 × 60 unit cells, 2 304 000 degrees of freedom, and only includes six spinful p orbitals, with 12 bands [51,57,58]. R = 5 random vectors and averaged over four disorder realizations. We include L · S SOC terms and first and second neighbor hoppings, with amplitudes that depend on the sublattices but differs from the VCA method, where the topological index not on the types of the atoms. We restrict the effect of disorder only depends on the average onsite energy of the cations. to different onsite energies on Sn and X sites. For more details, see Appendix E. VII. CONCLUSIONS Our method allows the computation of topological invari- VI. RESULTS ants of realistic 3D alloys. Disorder in the crystalline structure We define the Hamiltonians and perform the KPM expan- is present in naturally found and artificially grown com- sions using the Kwant software package [59]. The code to pounds, and it is inherent to substitutional alloys. However, reproduce the figures in this paper is available in Ref. [60]. a computationally efficient method to analyze topological First, we study the topological phase transition in the realistic properties of realistic disordered materials was missing. 18-orbital model of Pb1−x Snx Te. We build a tight-binding We apply our method to study the critical concentration of model with PBC that preserves reflection symmetry and con- Pb1−x Snx Te, and find an estimate that agrees with the most tains W × L110 × Lz unit cells, with 36 degrees of freedom precise experimental data [63], whereas it is not possible to each. For the largest system size used this means 13 824 000 reconcile any of the other theoretical studies (known to us) degrees of freedom in total. This model accurately reproduces with all the experimental data, especially Ref. [41]. the energetics, resulting in full bandwidth of about 25 eV The scaling of computational time with system size of our and band gap of less than 0.3 eV. In order to resolve the method is better than the scaling of other methods available gap that is multiple orders of magnitude smaller than the in the literature. By using the kernel polynomial method, bandwidth, we use M = 5000 moments in the calculation. We we achieve a computational time scaling of ξ d+1 with the use R = 5 random vectors and 12 disorder realizations each. localization length. Since we do not use eigenvalue solvers, This increases the time cost but not the memory cost of the only matrix-vector multiplication, the memory requirement algorithm. only scales linearly with the sample volume as ξ d . We perform finite-size collapse of the data [61] (see Ap- Beyond Chern numbers, our formalism allows calculation pendix F) and find the transition point at xc = 0.28(3), with of all Z-valued strong and weak topological invariants in all critical exponent ν = 0.9(6) accurately describing the transi- dimensions [50]. This method makes the automated discovery tion (see Fig. 2). This result improves significantly from the of topological alloys feasible, and can guide synthesis of new VCA result [32,34,40] and ab initio [33], but is consistent with alloys in the future. In this paper we use energetically accurate the most precise experimental data [39,41]. Experimental esti- tight-binding models obtained from ab initio calculations per- mates are obtained with only two values of the concentration formed on pure materials as input. Using these tight-binding x for works [35,40,62] focused on characterizing the trivial amplitudes for the atoms and bonds that appear in the alloy, and topological phases rather than determining the transition we generate large disordered samples with various concentra- point. tions. We are able to probe the topology, something that would To study a larger parameter space that includes other not be accessible with other methods. Simulation of small possible alloys of the X1−x Snx Te family that manifest the clusters with various disorder configurations is feasible using mirror Chern phase, we use the simplified six-orbital model. ab initio methods [33]; tight-binding parameters obtained for Besides the composition x, we also vary the onsite energy of all local environments would serve as a more accurate input the dopant cation X , approximating compounds with lighter for disordered models [47]. or heavier ions and similar electronic structures. Figure 3 This paper opens several directions for future research. shows the phase diagram. We find that the phase boundary Our method is directly applicable to all types of disorder, as 013229-4 COMPUTATION OF TOPOLOGICAL PHASE DIAGRAM OF … PHYSICAL REVIEW RESEARCH 2, 013229 (2020) well as quasicrystalline and amorphous systems in symmetry classes that admit the topological marker formalism [50]. This approach is not restricted to electrons in solids, and can be combined with finite element methods to analyze topol- ogy in disordered classical mechanical and photonic systems [64–67]. We expect our method to perform well in disordered time-reversal breaking Weyl semimetals with nonzero Hall conductivity; further refinements could extend it to the time- reversal invariant case. While we are not aware of a similar formulation of Z2 indices, KPM could be utilized to calculate quantized responses associated with these phases, such as the quantized magnetoelectric effect of 3D strong topological FIG. 4. Convergence of the projector operator, computed via a insulators [68]. A similar approach could also be applied to KPM expansion with M moments, applied to a vector |x0 , y0  located higher-order topological insulators to calculate multipole mo- at the center of a Chern insulator. The KPM expansion converges to the exact projector operator for M = ∞, and finite values of M ments of the charge density [69]. We expect that KPM could yield finite range approximations of the projector. The error of the be used to study a wide variety of related topics in condensed- approximation scales with exp(−ξ /M ) [29]. matter physics, such as probing localization, topological An- derson insulators, or numerical renormalization-group studies of the topological markers. During the preparation of the Concurrently, each order of the expansion increases the paper, our method served as a base to compute higher-order precision in energy of the approximation, therefore setting an topological invariants [70], and as an aide in the identification energy resolution for the expansion. To apply KPM, operators of topological phases of quasicrystals [71]. are rescaled such that the spectrum is in the [−1, 1] range. As a consequence, the order of the expansion required to resolve the mobility gap is proportional to W/, the full bandwidth ACKNOWLEDGMENTS divided by the mobility gap. This can be a large ratio even if ξ is small, for example, in an atomic insulator with a large We are grateful to M. D. Caio for discussions and in- range of on-site energies and vanishing hoppings. Typically sights on local topological markers, M. Hoskam for help W/ and ξ increase together near a phase transition where with implementing the tight-binding models, and B. Nijholt the gap closes, but the bandwidth changes slowly, and does for technical assistance with the numerical calculations. We not affect the scaling of the computational cost with ξ . thank S. Roche and A. Cummings for critical revision of the The kernel polynomial method provides a stable and effi- paper, and P. Brouwer for helpful discussions. This work was cient method to expand the action of any function of an opera- supported by European Research Council Starting Grant No. tor fˆ that depends on the Hamiltonian H and a set of parame- 638760, The Netherlands Organisation for Scientific Research ters λ, on a vector |v [29,30]. The expansion up to order M is (NWO/OCW) as part of the Frontiers of Nanoscience pro- gram, NWO VIDI Grant No. 680-47-53, the ONR, University  M of Chicago MRSEC through a Kadanoff-Rice postdoctoral fˆ(λ, H ) |v = gm μm (λ)Tm (H ) |v fellowship, and the European Union’s Horizon 2020 Research m=0 and Innovation program under Grant No. 824140.  M The project to study topological invariants using KPM = gm μm (λ) |vm  , was initiated by P.M.P.-P. and D.V., and the scope of the m=0 project was later refined using contributions from all authors. where gm are the Jackson kernel coefficients [29]. The P.M.P.-P. and D.V. wrote the code used for the numerical coefficients μm , called moments in the context of KPM calculations and performed calculations. D.V. performed the expansions, are defined as large-scale numerical calculations on the Pb1−x Snx Te tight-  1 ˆ binding models. All authors took part in the analysis of the 2 1 f (λ, E )Tm (E ) μm (λ) = √ dE , method and writing the paper. π 1 + δm,0 −1 1 − E2 and the vectors |vm  satisfy the recursion relation APPENDIX A: SPECTRAL PROJECTOR OPERATOR |v0  = |v , EXPANDED WITH THE KERNEL POLYNOMIAL METHOD |v1  = H |v0  , The band projector is expanded with the KPM [29] and it |vm+1  = 2H |vm  − |vm−1  . is a finite-range approximation, where the range depends lin- early on the number of moments used in the expansion. Each We approximate the projector operator defined as the step moment of the expansion of the projector operator applied to function: a vector is obtained recursively, by applying the Hamiltonian P̂(ε, H ) = θ (ε − H ), to an expanded vector of the previous iteration. The recursive algorithm effectively spreads a local vector to the neighboring  M sites via the hopping terms of the Hamiltonian, as depicted in P̂(ε, H ) = gm μm (ε)Tm (H ), Fig. 4. m=0 013229-5 DÁNIEL VARJAS et al. PHYSICAL REVIEW RESEARCH 2, 013229 (2020) and in this case the coefficients take the form  1 − arccos(ε) m=0 μm (ε) = −2 sin [m πarccos (ε)] . mπ m=0 Equipped with the KPM expanded projector, we proceed to evaluate matrix elements of topological markers. These are finite polynomials of P̂ and other sparse operators such as position and mirror. The matrix elements are evaluated by successive application of these operators to the states. The resulting memory cost scales linearly with the system size (number of degrees of freedom); by cumulatively summing up the expanded vectors for fixed EF , only a small number of sparse matrices and dense vectors are stored at any given FIG. 5. Off-diagonal matrix elements of the Chern operator in time. Most of the time cost comes from sparse matrix-vector two dimensions as a function of real-space distance in the units of multiplications, linear in the system size. The number of the localization length. Here we use a simple continuum model of operations is proportional to the number of moments M. a Chern insulator, discretized on a square lattice with various lattice constants a and fixed ξ . The collapse of the curves verifies the scaling form of the matrix elements that we use. APPENDIX B: SCALING OF STOCHASTIC TRACE We optimize the calculation further by utilizing the are interested in systems the size of which is proportional stochastic trace approximation to evaluate the trace. We take R to the localization length, so we will set L = cξ in the end. independent random-phase vectors |ri  that are only nonzero We introduce a as the lattice constant here for clarity, but it inside the region S, x, l|ri  = δx∈S exp(iφx,l,i ) with φx,l,i ∈ is an arbitrary reference length scale we can define in fully [0, 2π ] independent random phases for all sites and orbitals. disordered (e.g., amorphous) systems as well, for example, The trace of an operator Ô equals the expectation value as the typical spacing of sites. As it cancels from the final   result, this argument does not rely on the assumption of an 1 R underlying regular lattice. Tr Ô = E ri | Ô |ri  = E(Tr st Ô), (B1) We start with the Chern operator in two dimensions: R i=1 2π i where E denotes the expectation value over random vector Ĉ = [P̂x̂ P̂, P̂ŷP̂], (B3) a2 realizations and we introduced the notation Tr st Ô for the random variable giving the stochastic trace of operator Ô. where the a−2 prefactor is included to measure all distances The above equality is proved by using that the random phases in units of a; this way the sum of the diagonal entries of Ĉ are independent, hence E(ei(φx,l −φx ,l ) ) = δx,x δl,l and only the on a site coincides with the Chern number in a clean system. diagonal entries contribute to the expectation value. We numerically verify (see Fig. 5) that the off-diagonal matrix The standard deviation of the stochastic trace of an opera- elements scale as tor scales with the total square magnitude of the off-diagonal  entries which enter in the expectation value with random x−x y−y x, y| Ĉ |x , y  = f , , (B4) phases [29]: ξ ξ  where f is a dimensionless function and we suppressed the 1 σ (Tr st Ô) = |Ôi j |2 , (B2) dependence on the internal degrees of freedom. f is a quickly R i= j decaying function for (x − x )/ξ 1 in insulating systems, as matrix elements of P̂ also decay at the length scale of ξ . The Chern marker averaged over a square region S of size where σ (X ) = E(|X |2 ) − |E(X )|2 is the standard deviation L around the origin is given by and we used that σ (eiφx,l ) = 1. We also use that the standard deviation of the sum of independent random variables obeys    a 2 σ ( i Xi ) = C = Tr st Ĉ , (B5) i σ (Xi ) . 2 L We are concerned with the stochastic trace of topological markers, such as the Chern and mirror Chern operators. In where we still measure length in units of a. Substituting (B4) order to draw conclusions, we need to know the scaling of we find for the standard deviation of C the off-diagonal matrix elements with respect to the relevant   length scales in the problem. There are three length scales, the 1   a 2 r − r 2 lattice constant a, the localization length ξ , and the system σ (C) = √  f  R r=r ∈S  L ξ  size L (this we take to be the linear size of the subsystem where we take the partial trace; the overall system size is     a constant factor larger). We use units of a to measure the 1  a 4 1  r − r 2 f  = √ d 2 rd 2 r   other two distances, and, as explained in the main text, we R L a4 S ξ 013229-6 COMPUTATION OF TOPOLOGICAL PHASE DIAGRAM OF … PHYSICAL REVIEW RESEARCH 2, 013229 (2020) to the Chern number; the mirror Chern number is effectively a 2D invariant that we evaluate on a thick slab. In a clean system every plane parallel to a mirror plane is also a mirror plane, hence the matrix element can only depend on z + z . The mirror Chern marker averaged over a square region of size L is given by    a 2 CM = Tr st ĈM , (B10) L and we find using a similar derivation for the standard devia- tion (setting ξ = ξz = L/c)  c/2 1 σ (CM ) = √ d 3 r̃ d 3 r̃ | f (r̃, r̃ )|2 , (B11) c2 R −c/2 which is only dependent on the ratio of the localization length FIG. 6. Contributions per layer to the mirror Chern number. The and the system size c and the number of random vectors R. shade represents the standard deviation of the stochastic trace per In the numerical calculations we split the stochastic trace in layer, calculated with ten random vectors per layer. two halves, using two sets of random vectors, each localized in one half of the system separated by mirror planes. This   eliminates most of the large off-diagonal entries with z = −z , c/2 1 1 4 resulting in a constant factor reduction in the error. Splitting = √ ξ d 2 r̃ d 2 r̃ | f (r̃ − r̃ )|2 R L4 −c/2 the stochastic trace into more regions (e.g., separate for each  layer parallel to the mirror plane) results in further reduction c/2 1 in the error, at the cost of increased computational effort. The = √ d 2 r̃ d 2 r̃ | f (r̃ − r̃ )|2 . (B6) c2 R −c/2 overall scaling of the computational time with ξ for a fixed standard deviation is the same up to a constant factor for all of In the second line we took the limit of ξ a, so we can these schemes, ξ d+1 . replace sums over sites with integrals as x∈[−L/2,L/2] =  L/2 1/a −L/2 dx. In the third line we changed the integration APPENDIX C: GEOMETRY USED IN THE NUMERICS variables to r̃ = r/ξ . The result is only dependent on the ratio As described in the main text, we build a tight- of the localization length and the system size c = L/ξ and binding model with PBC using translation vectors W [1, 1, 0], the number of random vectors R, but not ξ . As the integral is L110 [1, −1, 0], and Lz [0, 0, 1]. This geometry preserves the proportional to the integration √ area, c2 for c 1, the overall reflection symmetry with [1, 1, 0] normal and contains W × scaling of the error with c is 1/ Rc2 . L110 × Lz unit cells with 36 degrees of freedom each. The We find a similar scaling for other topological markers, averaging region of the stochastic trace extends the full width such as the 3D winding number in chiral classes [48,49]. of the system in the [1, 1, 0] direction and contains half of the In general, the standard deviation of the stochastic trace linear size in the perpendicular directions, as depicted by the evaluation of the invariant depends only on the ratio of the inner box in Fig. 7. Imposing PBC in all directions eliminates system size and the localization length: gapless surface states, and the (mobility) gap guarantees that  the Fermi projector is short ranged.  1 ξ d However, imposing PBC in the direction normal to the σ (ν̂) ∝ . (B7) R L mirror planes results in two mirror invariant planes. As argued in the main text, in a sample with open boundary conditions in For the mirror Chern operator the other directions, this results in a doubling of the interface π modes, around the edges of the mirror invariant planes. On ĈM = M̂z [P̂x̂ P̂, P̂ŷP̂] (B8) the other hand, CM counts helical modes; while the 3D mirror a2 Chern number is given by the number of chiral modes, this the scaling is factor of 1/2 cancels the previous factor of 2. We conclude  that Eq. (6) is applicable to this 3D Mz symmetric geometry a x−x y−y z+z with PBC in z. x, y, z| ĈM |x , y , z  = f , , , (B9) ξz ξ ξ ξz APPENDIX D: 18-ORBITAL TIGHT-BINDING MODEL where ξz is the localization length in the z direction. This form is justified by the fact that the contributions to the We use the 18-orbital model of SnTe and PbTe derived in mirror Chern number are centered on the invariant planes, Ref. [32]. The cubic rock-salt structure has two sublattices A but are spread out on layers in a thickness proportional to ξz and C (referring to the anion and cation nature of the atoms (see Fig. 6). The total for all layers is, however, constant, occupying them), the first occupied by Te and the second by hence the a/ξz prefactor. This is the key difference compared Sn or Pb atoms. Each site hosts spinful s, p, and d orbitals, 18 013229-7 DÁNIEL VARJAS et al. PHYSICAL REVIEW RESEARCH 2, 013229 (2020) TABLE I. Tight-binding parameters in electron volts for SnTe and PbTe from Lent et al. [32]. Note that we use opposite sign conventions for Vspσ and Vpsσ . All other parameters not listed here vanish. SnTe PbTe Esc −6.578 −7.612 Esa −12.067 −11.002 E pc 1.659 3.195 E pa −0.167 −0.237 Edc 8.38 7.73 Eda 7.73 7.73 λ pc 0.592 1.500 FIG. 7. Geometry used for calculating the mirror Chern number. λ pa 0.564 0.428 The upper and lower halves of the slab have mirror image disorder Vssσ −0.510 −0.474 configurations with mirror plane M1 . Because of the PBC in the z Vspσ −0.949 −0.705 direction, there is a second mirror plane M2 . The central pane with Vpsσ 0.198 −0.633 stars is a schematic representation of the disorder, repeated in the Vppσ 2.218 2.066 x and y directions, and mirrored in z. The box in the center of Vppπ −0.446 −0.430 the sample shows the averaging region where the partial trace is Vpdσ −1.11 −1.29 evaluated, containing one period of the disorder configuration in the Vpdπ 0.624 0.835 x and y directions. Vd pσ −1.67 −1.59 Vd pπ 0.766 0.531 Vddσ −1.72 −1.35 degrees of freedom in total, with annihilation operators cl (l = Vddδ 0.618 0.668 0, 1, 2 for s, p, d orbitals, respectively), which is a vector of length 2, 6, or 10 depending on the value of l. The hopping terms are expressed as two-center inte- grals Hll md in the linear combination of atomic orbitals strength λl . The tight-binding Hamiltonian reads method [72], where l and l are the total angular momenta of the orbitals connected on the two sites and m is the   H = Elr c†lr · clr + λlr c†lr (Ll · S)clr angular momentum of the bonding along the bonding axis d l,r l,r (m = 0, 1, 2 for σ , π , δ bonding, respectively). The matrices  Hll md are 2(2l + 1) × 2(2l + 1) and are proportional to the + Vl,l ,m,r,r c†l r Hl lm(r −r) clr . (D1) identity in spin space. The onsite terms contain different onsite l,l ,m,r,r  energies for the various orbitals El and L · S SOC terms with The first term is the onsite energy, and it is the main source of disorder in our simulation. For sites on the C sublattice the type of the site (Sn or Pb) is chosen randomly with probability 1 − x and x. The value of Elr is assigned accordingly to be ElcSnTe and ElcPbTe , respectively. The super- scripts SnTe and PbTe refer to the two sets of parameters for the two pure materials. If r ∈ A, we use a weighted average Elr = [nElaSnTe + (6 − n)ElaPbTe ]/6 where n is the number of nearest-neighbor sites occupied by Sn atoms. The second term is the L · S spin-orbit coupling; Ll is the vector of angular momentum l operators (zero for l = 0). The values of λlr are assigned in the same fashion, depending on the type of atoms. The third term describes nearest-neighbor hopping terms in the [001] and equivalent crystal directions; the sum runs over all nearest-neighbor pairs with r ∈ A and r ∈ C. Depending on the atoms at sites r and r the value of Vl,l ,m,r,r is set to Vl,lSnTe PbTe ,m if one of the sites is Sn or Vl,l ,m if one of the sites is Pb. All of the onsite energies and hopping terms are spin independent; SOC only enters through the onsite SOC terms. We summarize the parameter values used for numerical results in Table I. Because of the identical outer-shell electronic structure of FIG. 8. Density of states as a function of x near the bulk gap in Sn and Pb, the alloy composition x does not affect the doping the 18-band model of Pb1−x Snx Te. The red line shows the placement level, therefore we set the Fermi level EF to ensure half filling of the Fermi level. for all compositions (see Fig. 8). 013229-8 COMPUTATION OF TOPOLOGICAL PHASE DIAGRAM OF … PHYSICAL REVIEW RESEARCH 2, 013229 (2020) TABLE II. Tight-binding parameters in electron volts for SnTe used in the six-orbital model. We use the same Hamiltonian as Sessi et al. [74], but the numerical values of the parameters differ due to different normalization and sign conventions. SnTe mTe −1.65 mSn 1.65 taa −0.5 tac = tca 0.9 tcc 0.5 λa −0.3 λc −0.3 APPENDIX E: SIX-ORBITAL TIGHT-BINDING MODEL FIG. 9. Logarithm of the cost function (xc , ν ) as a function of We adopt the six-orbital model of SnTe and PbTe orig- the parameters (xc , ν ). The region D allowed by uncertainties around inally described in Mitchell and Wallis [57] and used in the minimum is delimited by a thick black dashed line. Refs. [51,58,73,74]. Each site hosts spinful p orbitals, six degrees of freedom in total with a vector of annihilation operators c. This Hamiltonian is formally identical to (D1) but APPENDIX F: EVALUATION AND UNCERTAINTY OF THE only includes p orbitals and ppσ hopping, while the hopping FINITE-SIZE SCALING PARAMETERS range is extended to second neighbors. The tight-binding We start from the data CM = fi (x) of the mirror Chern Hamiltonian reads number as a function of the concentration of Pb1−x Snx Te where i = 1, . . . , N correspond to different system sizes Li .   H= mr c†r · cr + λr c†r (L · S)cr These data are obtained by averaging multiple disorder re- r r alizations. We wish to perform a finite-size scaling collapse  of the data. The changes of variables x̃i = (x − xc )Li1/ν define + tr,r c†r [1 − (d̂ r,r · L)2 ]cr . (E1) new functions x̃ → f˜i (x̃), where both the critical concentra- r,r  tion xc and the correlation length exponent ν are parameters to estimate so that the curves for different L collapse onto The first term is the onsite energy (also termed the “mass a universal master curve. Hence, we seek to minimize the term”), and it is the main source of disorder in our simulation. distance between f˜i (x̃) over their common domain of defini- mr takes the value of mTe on the A sublattice, while for the tion. To do so, we interpolate the discrete data to perform the B sublattice a value is chosen between mSn and mX with change of variable, evaluate f˜i (x̃) over a fixed interval, and probability 1 − x and x. The second term is the L · S spin-orbit compute the variance of the data. This defines a cost function coupling; its value depends on the sublattice only (identical (xc , ν) that is minimal at optimal parameters (xc∗ , ν ∗ )  for Sn and X ). The third term is a ppσ type of hopping [72] (0.28, 0.9). The cost function is plotted as a function of that only connects p orbitals oriented along the direction of the parameters (xc , ν) in Fig. 9. To evaluate the uncertainty the bond d̂ r,r . on the optimal parameters, we evaluate the cost function We include first neighbor [001] and second neighbor [110] (xc∗ , ν ∗ ) at the optimal parameters (xc∗ , ν ∗ ) for different hoppings, with amplitudes that depend on the sublattices. We disorder realizations [that give different values of the initial summarize the parameter values used for numerical results in data CM,L (x)]. The standard deviation gives an estimation of Table II. the uncertainty u( ) on the cost function. The domain D = Because of the identical outer-shell electronic structure of {(xc , ν) | (xc , ν)  (xc∗ , ν ∗ ) + u( )} where the cost func- Sn and Pb, the alloy composition x does not affect the doping tion is closer to its minimum than the uncertainty u( ) is ap- level, therefore we set the Fermi level EF to ensure half filling proximately an ellipse (see Fig. 9) from which the uncertain- for all compositions. ties u(xc )  0.03 and u(ν)  0.6 can be directly estimated. [1] X.-L. Qi and S.-C. Zhang, Rev. Mod. Phys. 83, 1057 (2011). [4] H. C. Po, A. Vishwanath, and H. Watanabe, Nat. Commun. 8, [2] C.-K. Chiu, J. C. Y. Teo, A. P. Schnyder, and S. Ryu, Rev. Mod. 50 (2017). Phys. 88, 035005 (2016). [5] B. Bradlyn, L. Elcoro, J. Cano, M. G. Vergniory, Z. Wang, C. [3] J. Kruthoff, J. de Boer, J. van Wezel, C. L. Kane, and R.-J. Felser, M. I. Aroyo, and B. A. Bernevig, Nature (London) 547, Slager, Phys. Rev. X 7, 041069 (2017). 298 (2017). 013229-9 DÁNIEL VARJAS et al. PHYSICAL REVIEW RESEARCH 2, 013229 (2020) [6] K. Yang, W. Setyawan, S. Wang, M. B. Nardelli, and S. [37] S.-Y. Xu, C. Liu, N. Alidoust, M. Neupane, D. Qian, I. Curtarolo, Nat. Mater. 11, 614 (2012). Belopolski, J. D. Denlinger, Y. J. Wang, H. Lin, L. A. Wray, [7] S. Curtarolo, G. L. W. Hart, M. B. Nardelli, N. Mingo, S. G. Landolt, B. Slomski, J. H. Dil, A. Marcinkova, E. Morosan, Sanvito, and O. Levy, Nat. Mater. 12, 191 (2013). Q. Gibson, R. Sankar, F. C. Chou, R. J. Cava, A. Bansil, and [8] Z. Zhang, R.-W. Zhang, X. Li, K. Koepernik, Y. Yao, and H. M. Z. Hasan, Nat. Commun. 3, 1192 (2012). Zhang, J. Phys. Chem. Lett. 9, 6224 (2018). [38] Y. Tanaka, Z. Ren, T. Sato, K. Nakayama, S. Souma, T. [9] F. Tang, H. C. Po, A. Vishwanath, and X. Wan, Nature (London) Takahashi, K. Segawa, and Y. Ando, Nat. Phys. 8, 800 (2012). 566, 486 (2019). [39] Y. Tanaka, T. Sato, K. Nakayama, S. Souma, T. Takahashi, Z. [10] F. Tang, H. C. Po, A. Vishwanath, and X. Wan, Nat. Phys. 15, Ren, M. Novak, K. Segawa, and Y. Ando, Phys. Rev. B 87, 470 (2019). 155105 (2013). [11] T. Zhang, Y. Jiang, Z. Song, H. Huang, Y. He, Z. Fang, H. Weng, [40] C. Yan, J. Liu, Y. Zang, J. Wang, Z. Wang, P. Wang, Z.-D. and C. Fang, Nature (London) 566, 475 (2019). Zhang, L. Wang, X. Ma, S. Ji, K. He, L. Fu, W. Duan, Q.-K. [12] M. G. Vergniory, L. Elcoro, C. Felser, N. Regnault, Xue, and X. Chen, Phys. Rev. Lett. 112, 186801 (2014). B. A. Bernevig, and Z. Wang, Nature (London) 566, 480 [41] R. Zhong, X. He, J. A. Schneeloch, C. Zhang, T. Liu, I. (2019). Pletikosić, T. Yilmaz, B. Sinkovic, Q. Li, W. Ku, T. Valla, J. M. [13] L. Fu and C. L. Kane, Phys. Rev. B 76, 045302 (2007). Tranquada, and G. Gu, Phys. Rev. B 91, 195321 (2015). [14] D. Hsieh, D. Qian, L. Wray, Y. Xia, Y. S. Hor, R. J. Cava, and [42] J. E. Avron and R. Seiler, Phys. Rev. Lett. 54, 259 (1985). M. Z. Hasan, Nature (London) 452, 970 (2008). [43] Q. Niu, D. J. Thouless, and Y.-S. Wu, Phys. Rev. B 31, 3372 [15] S. Chadov, J. Kiss, J. Kübler, and C. Felser, Phys. Status Solidi: (1985). Rapid Research Letters 7, 82 (2013). [44] I. C. Fulga, F. Hassler, and A. R. Akhmerov, Phys. Rev. B 85, [16] D. Sante, P. Barone, E. Plekhanov, S. Ciuchi, and S. Picozzi, 165409 (2012). Sci. Rep. 5, 11285 (2015). [45] A. George, SIAM J. Numer. Anal. 10, 345 (1973). [17] J. Li, R.-L. Chu, J. K. Jain, and S.-Q. Shen, Phys. Rev. Lett. [46] M. Aizenman and G. M. Graf, J. Phys. A 31, 6783 (1998). 102, 136806 (2009). [47] E. Prodan and W. Kohn, Proc. Natl. Acad. Sci. USA 102, 11635 [18] C. W. Groth, M. Wimmer, A. R. Akhmerov, J. Tworzydło, and (2005). C. W. J. Beenakker, Phys. Rev. Lett. 103, 196805 (2009). [48] J. Song and E. Prodan, Phys. Rev. B 89, 224203 (2014). [19] A. Kitaev, Ann. Phys. (NY) 321, 2 (2006). [49] I. Mondragon-Shem, T. L. Hughes, J. Song, and E. Prodan, [20] A. M. Essin and J. E. Moore, Phys. Rev. B 76, 165307 Phys. Rev. Lett. 113, 046802 (2014). (2007). [50] P. Perez-Piskunow, D. Varjas, M. Fruchart, M. Caio, and A. [21] T. A. Loring and M. B. Hastings, Europhys. Lett. 92, 67004 Akhmerov (unpublished). (2010). [51] T. H. Hsieh, H. Lin, J. Liu, W. Duan, A. Bansil, and L. Fu, Nat. [22] R. Bianco and R. Resta, Phys. Rev. B 84, 241106(R) (2011). Commun. 3, 982 (2012). [23] T. A. Loring, Ann. Phys. (NY) 356, 383 (2015). [52] We expect that this will not happen with a regular short-range [24] E. Prodan, Appl. Math. Res. eXpress 2013, 176 (2013). correlated disorder. [25] E. Prodan and H. Schulz-Baldes, Bulk and Boundary Invari- [53] L. Fu and C. L. Kane, Phys. Rev. Lett. 109, 246605 (2012). ants for Complex Topological Insulators (Springer, New York, [54] Y. Ando and L. Fu, Annu. Rev. Condens. Matter Phys. 6, 361 2016). (2015). [26] Y. Akagi, H. Katsura, and T. Koma, J. Phys. Soc. Jpn. 86, [55] Y. Okada, M. Serbyn, H. Lin, D. Walkup, W. Zhou, C. Dhital, 123710 (2017). M. Neupane, S. Xu, Y. J. Wang, R. Sankar, F. Chou, A. Bansil, [27] H. Katsura and T. Koma, J. Math. Phys. 59, 031903 (2018). M. Z. Hasan, S. D. Wilson, L. Fu, and V. Madhavan, Science [28] I. C. Fulga, B. van Heck, J. M. Edge, and A. R. Akhmerov, 341, 1496 (2013). Phys. Rev. B 89, 155424 (2014). [56] C. Fang, M. J. Gilbert, and B. A. Bernevig, Phys. Rev. Lett. 112, [29] A. Weiße, G. Wellein, A. Alvermann, and H. Fehske, Rev. Mod. 046801 (2014). Phys. 78, 275 (2006). [57] D. L. Mitchell and R. F. Wallis, Phys. Rev. 151, 581 (1966). [30] J. H. García, L. Covaci, and T. G. Rappoport, Phys. Rev. Lett. [58] I. C. Fulga, N. Avraham, H. Beidenkopf, and A. Stern, Phys. 114, 116602 (2015). Rev. B 94, 125405 (2016). [31] D. Carvalho, N. A. García-Martínez, J. L. Lado, and J. [59] C. W. Groth, M. Wimmer, A. R. Akhmerov, and X. Waintal, Fernández-Rossier, Phys. Rev. B 97, 115453 (2018). New J. Phys. 16, 063065 (2014). [32] C. S. Lent, M. A. Bowen, J. D. Dow, R. S. Allgaier, O. F. [60] D. Varjas, M. Fruchart, A. R. Akhmerov, and P. Perez-Piskunow Sankey, and E. S. Ho, Superlattices Microstruct. 2, 491 (2019), doi:10.5281/zenodo.2667604. (1986). [61] M. D. Caio, G. Möller, N. R. Cooper, and M. J. Bhaseen, Nat. [33] X. Gao and M. S. Daw, Phys. Rev. B 77, 033103 (2008). Phys. 15, 257 (2019). [34] P. Dziawa, B. J. Kowalski, K. Dybko, R. Buczko, A. [62] T. Phuphachong, B. A. Assaf, V. V. Volobuev, G. Bauer, G. Szczerbakow, M. Szot, E. Łusakowska, T. Balasubramanian, Springholz, L. A. De Vaulchier, and Y. Guldner, J. Phys.: Conf. B. M. Wojek, M. H. Berntsen, O. Tjernberg, and T. Story, Nat. Ser. 864, 012038 (2017). Mater. 11, 1023 (2012). [63] In the case of Ref. [40], the ARPES data for x = 0.3 are not [35] J. O. Dimmock, I. Melngailis, and A. J. Strauss, Phys. Rev. Lett. conclusive, but they are interpreted with the help of VCA theory 16, 1193 (1966). as a trivial phase. We claim that for x = 0.3 the experimental [36] J. C. Y. Teo, L. Fu, and C. L. Kane, Phys. Rev. B 78, 045426 data and the theoretical results presented in Ref. [40] are not (2008). sufficient to ensure that x = 0.3 lies on the trivial phase. 013229-10 COMPUTATION OF TOPOLOGICAL PHASE DIAGRAM OF … PHYSICAL REVIEW RESEARCH 2, 013229 (2020) [64] S. Stützer, Y. Plotnik, Y. Lumer, P. Titum, N. H. Lindner, M. [70] O. Pozo, C. Repellin, and A. G. Grushin, Phys. Rev. Lett. 123, Segev, M. C. Rechtsman, and A. Szameit, Nature (London) 560, 247401 (2019). 461 (2018). [71] D. Varjas, A. Lau, K. Pöyhönen, A. R. Akhmerov, D. I. Pikulin, [65] N. P. Mitchell, L. M. Nash, D. Hexner, A. M. Turner, and and I. C. Fulga, Phys. Rev. Lett. 123, 196401 (2019). W. T. M. Irvine, Nat. Phys. 14, 380 (2018). [72] J. C. Slater and G. F. Koster, Phys. Rev. 94, 1498 (1954). [66] J. Cha, K. Woo Kim, and C. Daraio, Nature (London) 564, 229 [73] F. Schindler, A. M. Cook, M. G. Vergniory, Z. Wang, S. S. P. (2018). Parkin, B. A. Bernevig, and T. Neupert, Sci. Adv. 4, eaat0346 [67] D. Zhou, L. Zhang, and X. Mao, Phys. Rev. Lett. 120, 068003 (2018). (2018). [74] P. Sessi, D. Di Sante, A. Szczerbakow, F. Glott, S. Wilfert, [68] R. Bianco and R. Resta, Phys. Rev. Lett. 110, 087202 (2013). H. Schmidt, T. Bathon, P. Dziawa, M. Greiter, T. Neupert, G. [69] W. A. Benalcazar, B. A. Bernevig, and T. L. Hughes, Science Sangiovanni, T. Story, R. Thomale, and M. Bode, Science 354, 357, 61 (2017). 1269 (2016). 013229-11