Computation of topological phase diagram of disordered Pb$_{1-x}$Sn$_{x}$Te using the kernel polynomial method

TL;DR

This paper introduces a kernel polynomial method-based algorithm to compute topological invariants in large, disordered systems, enabling detailed phase diagram analysis of complex alloys like Pb$_{1-x}$Sn$_{x}$Te.

Contribution

The paper develops a scalable algorithm for topological invariant calculation in inhomogeneous systems, allowing analysis of large disordered materials with over 10 million degrees of freedom.

Findings

Determined the topological phase diagram of Pb$_{1-x}$Sn$_{x}$Te.

Identified the critical concentration for the topological phase transition.

Enabled analysis of large, disordered systems with high computational efficiency.

Abstract

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 $10^7$ degrees of freedom. Our method enables the study of large complex compounds, where disorder is inherent to the system. We use it to analyse Pb$_{1-x}$Sn$_{x}$Te and tighten the critical concentration for the phase transition.

Tables

Table 1. Table 1: Tight-binding parameters in electronvolts for \ce SnTe and \ce PbTe from Lent  et al. ( 1986 ) . Note that we use opposite sign convention for V s ​ p ​ σ subscript 𝑉 𝑠 𝑝 𝜎 V_{sp\sigma} and V p ​ s ​ σ subscript 𝑉 𝑝 𝑠 𝜎 V_{ps\sigma} . All other parameters not listed here vanish.

	\ceSnTe	\cePbTe
$E_{sc}$	$-6.578$	$-7.612$
$E_{sa}$	$-12.067$	$-11.002$
$E_{pc}$	$1.659$	$3.195$
$E_{pa}$	$-0.167$	$-0.237$
$E_{dc}$	$8.38$	$7.73$
$E_{da}$	$7.73$	$7.73$
${\lambda }_{pc}$	$0.592$	$1.500$
${\lambda }_{pa}$	$0.564$	$0.428$
$V_{ss\sigma }$	$-0.510$	$-0.474$
$V_{sp\sigma }$	$-0.949$	$-0.705$
$V_{ps\sigma }$	$0.198$	$-0.633$
$V_{pp\sigma }$	$2.218$	$2.066$
$V_{pp\pi }$	$-0.446$	$-0.430$
$V_{pd\sigma }$	$-1.11$	$-1.29$
$V_{pd\pi }$	$0.624$	$0.835$
$V_{dp\sigma }$	$-1.67$	$-1.59$
$V_{dp\pi }$	$0.766$	$0.531$
$V_{dd\sigma }$	$-1.72$	$-1.35$
$V_{dd\delta }$	$0.618$	$0.668$

Table 2. Table 2: Tight-binding parameters in electronvolts for \ce SnTe used in the 6-orbital model. We use the same Hamiltonian as Ref.  Sessi  et al. ( 2016 ) , but the numerical values of the parameters differ due to different normalization and sign conventions.

	\ceSnTe
$m_{\text{ce}Te}$	$-1.65$
$m_{\text{ce}Sn}$	$1.65$
$t_{aa}$	$-0.5$
$t_{ac}=t_{ca}$	$0.9$
$t_{cc}$	$0.5$
${\lambda }_a$	$-0.3$
${\lambda }_c$	$-0.3$

Equations

$$\hat{P}=\theta (E_F-\hat{H})=\sum _{n:E_n<E_F}|n⟩⟨n|$$

$$\nu ={\mathrm{Tr}}_S(\hat{\nu })=\frac{1}{|S|}\sum _{\lambda ,\mathbf{x}\in S}⟨\mathbf{x},\lambda |\hat{\nu }|\mathbf{x},\lambda ⟩,$$

$$C={\mathrm{Tr}}_A\hat{C}=2\pi \mathrm{i}{\mathrm{Tr}}_A\left[\hat{P}\hat{x}\hat{P},\hat{P}\hat{y}\hat{P}\right].$$

$${\mathrm{Tr}}_S(\hat{\mathcal{O}})\approx \frac{1}{R|S|}\sum _{i=1}^R⟨r_i|\hat{\mathcal{O}}|r_i⟩,$$

$$C_{\pm }=2\pi \mathrm{i}{\mathrm{Tr}}_A\left[{\tilde{x}}_{\pm },{\tilde{y}}_{\pm }\right],$$

$$C_M=\pi {\mathrm{Tr}}_A\left({\hat{M}}_z\left[\hat{P}\hat{x}\hat{P},\hat{P}\hat{y}\hat{P}\right]\right).$$

$$\hat{f}(\lambda ,H)|v⟩$$

$$=\sum _{m=0}^M{\mu }_m(\lambda )|v_m⟩,$$

$${\mu }_m(\lambda )=\frac{2}{\pi }\frac{1}{1+{\delta }_{m,0}}{\int }_{-1}^1\frac{\hat{f}(\lambda ,E)T_m(E)}{\sqrt{1-E^2}}\mathrm{d}E,$$

$$|v_0⟩$$

$$|v_1⟩$$

$$|v_{m+1}⟩$$

$$\hat{P}(\epsilon ,H)$$

$${\mu }_m(\epsilon )=\{\begin{array}{ll}1-\frac{1}{\arccos (\epsilon )}& m=0\\ \frac{-2\sin \left(m\arccos (\epsilon )\right)}{m\pi }& m\ne 0\end{array}$$

$$\mathrm{Tr}\hat{\mathcal{O}}=\text{mathbbm}E\left(\frac{1}{R}\sum _{i=1}^R⟨r_i|\hat{\mathcal{O}}|r_i⟩\right)=\text{mathbbm}E\left({\mathrm{Tr}}_{\mathrm{st}}\hat{\mathcal{O}}\right),$$

$$\sigma ({\mathrm{Tr}}_{\mathrm{st}}\hat{\mathcal{O}})=\sqrt{\frac{1}{R}\sum _{i\ne j}|{\hat{\mathcal{O}}}_{ij}{|}^2},$$

$$\hat{C}=\frac{2\pi \mathrm{i}}{a^2}\left[\hat{P}\hat{x}\hat{P},\hat{P}\hat{y}\hat{P}\right],$$

$$⟨x,y|\hat{C}|x^{\prime },y^{\prime }⟩=f\left(\frac{x-x^{\prime }}{\xi },\frac{y-y^{\prime }}{\xi }\right),$$

$$C={\mathrm{Tr}}_{\mathrm{st}}\left[{\left(\frac{a}{L}\right)}^2\hat{C}\right],$$

$$\sigma \left(C\right)$$

$$\sigma (\hat{\nu })\propto \sqrt{\frac{1}{R}{\left(\frac{\xi }{L}\right)}^d}.$$

$${\hat{C}}_M=\frac{\pi }{a^2}{\hat{M}}_z\left[\hat{P}\hat{x}\hat{P},\hat{P}\hat{y}\hat{P}\right]$$

$$⟨x,y,z|{\hat{C}}_M|x^{\prime },y^{\prime },z^{\prime }⟩=\frac{a}{{\xi }_z}f\left(\frac{x-x^{\prime }}{\xi },\frac{y-y^{\prime }}{\xi },\frac{z+z^{\prime }}{{\xi }_z}\right),$$

$$C_M={\mathrm{Tr}}_{\mathrm{st}}\left[{\left(\frac{a}{L}\right)}^2{\hat{C}}_M\right],$$

$$\sigma \left(C_M\right)=\frac{1}{c^2\sqrt{R}}\sqrt{{\int }_{-c/2}^{c/2}{d}^3\tilde{\mathbf{r}}{d}^3{\tilde{\mathbf{r}}}^{\prime }{\left|f(\tilde{\mathbf{r}},{\tilde{\mathbf{r}}}^{\prime })\right|}^2},$$

$$\begin{gathered} H=\sum _{l,\mathbf{r}}{E}_{l\mathbf{r}}{\mathbf{c}}_{l\mathbf{r}}^{†}\cdot {\mathbf{c}}_{l\mathbf{r}}+\sum _{l,\mathbf{r}}{\lambda }_{l\mathbf{r}}{\mathbf{c}}_{l\mathbf{r}}^{†}\left({\mathbf{L}}_l\cdot \mathbf{S}\right){\mathbf{c}}_{l\mathbf{r}} \\ +\sum _{l,l^{\prime },m,⟨\mathbf{r},{\mathbf{r}}^{\prime }⟩}{V}_{l,l^{\prime },m,\mathbf{r},{\mathbf{r}}^{\prime }}{\mathbf{c}}_{l^{\prime }{\mathbf{r}}^{\prime }}^{†}{H}_{l^{\prime }lm({\mathbf{r}}^{\prime }-\mathbf{r})}{\mathbf{c}}_{l\mathbf{r}}. \end{gathered}$$

$$\begin{gathered} H=\sum _{\mathbf{r}}{m}_{\mathbf{r}}{\mathbf{c}}_{\mathbf{r}}^{†}\cdot {\mathbf{c}}_{\mathbf{r}}+\sum _{\mathbf{r}}{\lambda }_{\mathbf{r}}{\mathbf{c}}_{\mathbf{r}}^{†}\left(\mathbf{L}\cdot \mathbf{S}\right){\mathbf{c}}_{\mathbf{r}} \\ +\sum _{⟨⟨\mathbf{r},{\mathbf{r}}^{\prime }⟩⟩}{t}_{\mathbf{r},{\mathbf{r}}^{\prime }}{\mathbf{c}}_{{\mathbf{r}}^{\prime }}^{†}\left[\text{mathbbm}1-{\left({\hat{\mathbf{d}}}_{\mathbf{r},{\mathbf{r}}^{\prime }}\cdot \mathbf{L}\right)}^2\right]{\mathbf{c}}_{\mathbf{r}}. \end{gathered}$$

Full text

.0.0.1 People search for topological materials using high throughput algorithms.

Introduction. — Topological materials have attracted continuing interest from both the fundamental physics and the material science communities for the last decade Qi and Zhang (2011); Chiu et al. (2016). The program to theoretically classify non-interacting crystalline insulators has been completed, tabulating possible topological phases in all space groups Kruthoff et al. (2017); Po et al. (2017); Bradlyn et al. (2017). Recent efforts focus on automated high-throughput methods to discover and classify new topological materials Yang et al. (2012); Curtarolo et al. (2013); Zhang et al. (2018), culminating in the production of comprehensive databases of topological insulators and semimetals Tang et al. (2019a, b); Zhang et al. (2019); Vergniory et al. (2019); Zhang et al. (2018).

I.0.0.1 Alloys are important.

However, not all topological insulators are compounds with perfect stoichiometry. The first three-dimensional topological insulator to be predicted Fu and Kane (2007) and experimentally realized Hsieh et al. (2008) was an alloy—\ceBi_xSb_1-x. These systems are usually studied using the virtual crystal approximation or the coherent potential approximation, which approximate an alloy by a perfect crystal Chadov et al. (2013); Sante et al. (2015). This approach ignores the intrinsic disorder in alloys, and it is insufficient to explain topological transitions that appear at strong disorder Li et al. (2009); Groth et al. (2009), or accurately find critical concentrations.

I.0.0.2 Alloys need large supercell, this is expensive.

The topological invariant converges to its bulk value in samples larger than the localization length $\xi$. This is the main limitation in resolving topological phase transitions as $\xi$ diverges. Therefore, the asymptotic scaling of the computational cost with $\xi$ is the main distinction between different numerical approaches.

I.0.0.3 Several alternatives exist, but their performance is unsuitable for 3D.

Available methods to compute topological invariants either apply the momentum-space Berry curvature formalism to periodic systems with a disordered supercell, or use a real space formulation on a large finite sample Kitaev (2006); Essin and Moore (2007); Loring and Hastings (2010); Bianco and Resta (2011); Loring (2015); Prodan (2013); Prodan and Schulz-Baldes (2016); Akagi et al. (2017); Katsura and Koma (2018). However, these methods involve solving at least one eigenvalue equation with size equal to the number of degrees of freedom, resulting in the complexity $\xi^{3d}$ in $d$ dimensions. This restricts the applicability of such methods to small system sizes, especially in three dimensions (3D). To our knowledge, the most efficient method in 3D is the scattering matrix approach Fulga et al. (2014), with a complexity scaling as $\xi^{3(d-1)}$, allowing for maximum sizes of $5\times 10^{5}$ degrees of freedom.

I.0.0.4 We develop a new method based on KPM.

We present an algorithm to efficiently identify topological phases of strongly disordered systems using the kernel polynomial method (KPM) Weiße et al. (2006); García et al. (2015); Carvalho et al. (2018), an approximation based on a polynomial expansion of the quantities of interest. All topological properties of a non-interacting system of electrons are encoded in the projector on the occupied states (spectral projector), which is efficiently approximated using KPM. Our algorithm builds on the method of topological markers Bianco and Resta (2011) to construct a topological invariant as the trace of an operator. Since topological invariants are integers, it is sufficient to reduce the statistical uncertainty below $1/2$ to obtain the exact value. In particular, the stochastic evaluation of traces Weiße et al. (2006) from a small number of random vectors, combined with KPM, is suited for this task. With $\xi^{d+1}$ scaling of the computational effort, it is the most efficient method in three dimensions.

I.0.0.5 We use this method to analyze large samples of real materials.

As a concrete example, we apply our method to lead tin telluride \cePb_1-xSn_xTe alloys—three-dimensional topological crystalline insulators characterized by a mirror Chern number. Thanks to the efficiency of our algorithm, we analyze 3D systems with linear sizes of over a hundred lattice constants ($L>100$), and more than $10^{7}$ degrees of freedom (see Fig. 1). In contrast to previous theoretical estimations for the case of \cePb_1-xSn_xTe Lent et al. (1986); Gao and Daw (2008); Dziawa et al. (2012), we find a critical concentration that matches the one found experimentally Dimmock et al. (1966); Teo et al. (2008); Xu et al. (2012); Tanaka et al. (2012, 2013); Yan et al. (2014); Zhong et al. (2015).

I.0.0.6 Linear size of a disordered sample must be of order $\xi$.

Review of existing algorithms. — We focus on the vicinity of disorder-driven phase transitions, where the localization length $\xi$ diverges. The finite but potentially large value of $\xi$ defines the relevant length scale of our problem. Two regions that are closer than $\xi$ feel each other’s presence. Therefore, with open boundary conditions, the bulk has to be further than $\xi$ to the edge; analogously, with periodic boundary conditions when the linear size of the system is smaller than $\xi$ it hosts states whose extent is larger than the system size. These states span the whole system and overlap with themselves because of the finite size. In order to simulate the bulk, the system needs to have a linear size $L\gtrsim\xi$. Fluctuations of local quantities resulting from disorder scale with $\xi$, so averaging over a larger sample provides a good approximation of the thermodynamic limit.

II.0.0.1 Most methods require full diagonalization and cost $\xi^{3d}$.

The momentum-space formalism of topological invariants applies to disordered systems by studying a periodic system with a large disordered supercell of volume $\xi^{d}$. This is equivalent to taking a finite torus and threading fluxes through its cycles Avron and Seiler (1985); Niu et al. (1985); Essin and Moore (2007). The final formula of the invariant is identical to the momentum space Berry curvature treatment applied to the supercell. Other approaches include the Bott index Loring and Hastings (2010), topological markers Bianco and Resta (2011), pseudospectra Loring (2015), and noncommutative index theorems Prodan (2013); Prodan and Schulz-Baldes (2016); Akagi et al. (2017); Katsura and Koma (2018). All of these methods involve diagonalization of a matrix of size proportional to the volume of the system. Diagonalization scales as $N^{3}$ with the number of degrees of freedom $N$, so the computational cost of such methods is order $\xi^{3d}$, restricting them to small system sizes in three dimensions.

II.0.0.2 Scattering approach is faster and requires $\xi^{3d-3}$.

The scattering invariant formalism Fulga et al. (2012) avoids full diagonalization 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 George (1973) and scales as $\xi^{3d-3}$ for $d>1$.

II.0.0.3 Projector can be approximated using KPM with a finite range of the expansion.

General strategy. — Computing the exact spectral projector

[TABLE]

by full diagonalization of the Hamiltonian is numerically expensive. Instead, we approximate the projector using the kernel polynomial method with the Jackson kernel Weiße et al. (2006), detailed in Appendix A. For a $d$-dimensional system of linear size $L$ the computational cost scales linearly with the number of degrees of freedom $L^{d}$, and with the number of moments $M$—the order of the expansion. The order of the expansion sets a real-space cutoff in the approximate projector, which is an $M$’th order polynomial of the finite-range Hamiltonian. In an insulating system, the projector is a local operator with matrix elements $\bra{\bm{x}}\hat{P}\ket{\bm{x^{\prime}}}\propto\exp\left(-|\bm{x}-\bm{x^{\prime}}|/\xi\right)$ that have a decay length $\xi$. Hence, the error of the approximation scales as $\exp(-M/\xi)$, and the number of moments necessary for fixed precision scales linearly with the localization length as $M\sim\xi$ Aizenman and Graf (1998); Prodan and Kohn (2005). See Appendix A.

III.0.0.1 Topological markers converge to the bulk invariant for a large subsystem.

We use the topological marker formalism introduced by Bianco and Resta (2011). All $\mathbb{Z}$-valued topological markers are a partial trace per unit volume of a local operator $\hat{\nu}$

[TABLE]

where the sum runs over the sites $\bm{x}$ inside the subsystem $S$ with volume $|S|$, and their internal degrees of freedom $\lambda$. The operator $\hat{\nu}$ is a polynomial of the spectral projector, position and symmetry operators, such that $\nu$ is dimensionless and independent of the detailed energetics or the overall length scale of the system. The marker coincides with the momentum-space invariant in periodic systems, and converges to a quantized integer for large $S$ in insulating homogeneous disordered systems Bianco and Resta (2011).

III.0.0.2 We illustrate the general approach to $Z$ invariants by considering Chern number, but we do specify how to do $Z_{2}$.

An example of topological marker is the real space expression for the Chern number Bianco and Resta (2011):

[TABLE]

Here, $A$ is the area of the subsystem, $\hat{x}$ and $\hat{y}$ are the two components of the position operator, and $[\cdot,\cdot]$ is the commutator. Topological markers for all strong and weak $\mathbb{Z}$-valued topological invariants have similar algebraic expressions of the projected position operators Song and Prodan (2014); Mondragon-Shem et al. (2014); Chiu et al. (2016), making a straightforward application of our method to these cases Perez-Piskunow et al. . We are not aware of similar formulations of $\mathbb{Z}_{2}$ topological indices suitable for KPM.

III.0.0.3 To average the volume contribution to the invariant we apply stochastic trace.

To estimate the trace per volume, we use the stochastic trace approximation Weiße et al. (2006)

[TABLE]

where $|r_{i}\rangle$ are random phase vectors localized in the region $S$. The standard error of this approximation scales as $\sqrt{\xi^{d}/(R|S|)}$, meaning that the number of random vectors $R$ required for a given precision is constant if the system size is proportional to the localization length (see Appendix B).

III.0.0.4 Repeating the supercell twice eliminates boundary fluctuations of the invariant.

We build a supercell of size $L>\xi$. To ensure that the invariant $\nu$ obtained with the trace of the local marker converges to the momentum space topological invariant we must repeat the supercell with the disorder realization over the whole space. Since the approximation of the local operator $\hat{\nu}$ is localized with $\xi<L$, it is sufficient to repeat two times the supercell in every spatial direction, and compute the topological invariant as an average of the topological marker over the central $L^{d}$ volume. For a detailed description, see Appendix C.

III.0.0.5 Combining all factors our method scales as $\xi^{d+1}$.

The resulting complexity of the computation depends linearly on the number of random vectors $R$ used (typically of order $1$), on the number of moments $M$, and on the number of sites of the system $L^{d}$. We use a sparse representation of the short-ranged Hamiltonian. As a result the memory requirement scales linearly with the system size $L^{d}$ and is independent of other parameters. Setting all quantities to their minimal values ($L\sim\xi$ and $M\sim\xi$), results in an algorithm with a computational cost scaling of $\xi^{d+1}$.

III.0.0.6 Our method offers an edge in analysis of 3D materials, we therefore apply it to mirror Chern alloys.

Application to mirror Chern number. — Our method provides better scaling than existing approaches in $d\geq 3$. We apply it to disordered 3D mirror Chern insulators. These are a widely studied class of topological crystalline materials with a $\mathbb{Z}$ topological classification that relies on reflection symmetry Hsieh et al. (2012). Several experimental realizations are known, including alloys Xu et al. (2012); Dziawa et al. (2012); Tanaka et al. (2012, 2013); Yan et al. (2014); Zhong et al. (2015).

IV.0.0.1 We rewrite the mirror Chern number in the topological marker formalism.

In a reflection-symmetric system of fermions, all wave functions are eigenstates of the mirror operator $\hat{M}_{z}$, with eigenvalues $\pm\mathrm{i}$. The Chern numbers $C_{\pm}$ for mirror-even and mirror-odd wave functions are

[TABLE]

where ${\tilde{x}_{\pm}=\hat{M}_{\pm}\hat{P}\hat{x}\hat{P}\hat{M}_{\pm}}$, and ${\tilde{y}_{\pm}=\hat{M}_{\pm}\hat{P}\hat{y}\hat{P}\hat{M}_{\pm}}$ are the projected position operators restricted to the mirror-even or mirror-odd subspaces. Here $\hat{M}_{\pm}$ are the projectors on the mirror-even and mirror-odd subspaces and $A$ is the area in the $xy$ plane. The total Chern number $C$ is the sum of the Chern numbers for each subspace $C=C_{+}+C_{-}$, while the mirror Chern number equals to their difference $C_{M}=\left(C_{+}-C_{-}\right)/2$. In the presence of time-reversal invariance, the total Chern number vanishes (${C_{+}=-C_{-}}$), and the mirror Chern number ${C_{M}=C_{+}}$ counts the helical surface modes. Since the mirror operator ${\hat{M}_{z}=\mathrm{i}\left(\hat{M}_{+}-\hat{M}_{-}\right)}$ commutes with the projector $\hat{P}$ and position operators $\hat{x}$ and $\hat{y}$, we express the mirror Chern number as

[TABLE]

IV.0.0.2 In order to compute mirror Chern number we enforce reflection symmetry.

In order to compute the mirror Chern number we consider a system with a disorder configuration that is mirror symmetric across the $M_{1}$ plane at $z=0$. The PBC in the $z$-direction results in another mirror plane $M_{2}$ across the boundary (see Appendix C). The bulk of this system is locally indistinguishable from a sample without reflection symmetry except for the two mirror planes $M_{1}$ and $M_{2}$. Therefore as long as the two mirror planes do not undergo a two dimensional topological transition, the mirror Chern number of the symmetric sample equals that of the bulk system111We expect that this will not happen with a regular short-range correlated disorder..

IV.0.0.3 \cePb_1-xSn_xTe is the canonical mirror Chern insulator.

Tight binding model of $\text{Pb}_{1-x}\text{Sn}_{x}\text{Te}$. — Topological crystalline insulators (TCI) protected by reflection symmetry Teo et al. (2008) were theoretically predicted Hsieh et al. (2012) and experimentally observed Tanaka et al. (2012); Dziawa et al. (2012); Xu et al. (2012); Tanaka et al. (2013); Yan et al. (2014); Zhong et al. (2015) in \cePb_1-xSn_xTe alloys. They host metallic surface states on the surfaces that are symmetric with respect to the mirror plane Hsieh et al. (2012); Fu and Kane (2012); Chiu et al. (2016). Lead tin telluride was studied using either the virtual crystal approximation (VCA) Lent et al. (1986); Dziawa et al. (2012) or ab initio methods Gao and Daw (2008), finding a gap closing and phase transition near $x=0.35$, or $x=0.23$, respectively. We use a tight-binding approach that captures long-range correlations, and find a different critical concentration.

IV.0.0.4 We consider substitutional disorder as appropriate for alloys.

We consider the substitutional disorder of the lead tin telluride alloy \cePb_1-xSn_xTe coming from replacing some \cePb ions for \ceSn ions. This disorder is nonmagnetic, and it preserves the reflection symmetry on average, which is sufficient to protect the gapless surface states Fulga et al. (2014); Ando and Fu (2015). We disregard other types of symmetry breaking disorder appearing naturally in \cePb_1-xSn_xTe Ando and Fu (2015), such as ferroelectric structural distortion Ando and Fu (2015); Hsieh et al. (2012); Okada et al. (2013), and magnetic dopants Ando and Fu (2015); Hsieh et al. (2012); Fang et al. (2014).

IV.0.0.5 We use the full LCAO 18-orbital model and a simplified 6 orbital model.

In our investigation we use two atomistic tight-binding models. The first one includes 18 spinful $s$, $p$ and $d$ orbitals on both sublattices, with 36 bands in total. This model accurately describes the energetics, using tight-binding parameters for both \ceSnTe and \cePbTe derived from ab initio simulations Lent et al. (1986).

To simulate the alloy, we substitute randomly \ceSn for \ceTe with probability $x$. We incorporate substitutional disorder by using the hopping amplitudes of \ceSnTe for \ceSn–\ceTe bonds and \cePbTe amplitudes for \cePb–\ceTe bonds. The onsite parameters of \ceTe atoms are slightly different in \ceSnTe and \cePbTe; we therefore use a weighted average of these depending on the local environment. We also use the appropriate onsite terms, including ${\bm{L}\cdot\bm{S}}$ spin-orbit coupling (SOC), depending on the type of the \ceSn or \cePb atom. Further details can be found in Appendix D.

When investigating the onsite energy dependent phase diagram of \ceX_1-xSn_xTe alloys, we use a simplified model that only includes 6 spinful $p$ orbitals, with 12 bands Mitchell and Wallis (1966); Hsieh et al. (2012); Fulga et al. (2016). We include ${\bm{L}\cdot\bm{S}}$ SOC terms, first and second neighbor hoppings, with amplitudes that depend on the sublattices but not on the types of the atoms. We restrict the effect of disorder to different onsite energies on \ceSn and \ceX sites. For more details, see Appendix E.

IV.0.0.6 We run heavy simulations using Kwant and KPM.

Results. — We define the Hamiltonians and perform the KPM expansions using the Kwant software package Groth et al. (2014). The code to reproduce the figures in this article is available in Ref. Varjas et al. (2019). First, we study the topological phase transition in the realistic 18-orbital model of \cePb_1-xSn_xTe. We build a tight-binding model with PBC that preserves reflection symmetry and contains $W\times L_{1\overline{1}0}\times L_{z}$ unit cells, with 36 degrees of freedom each. For the largest system size used this means $13\,824\,000$ degrees of freedom in total. This model accurately reproduces the energetics, resulting in full bandwidth of about $25\text{\,}\mathrm{eV}$ and band gap of less than $0.3\text{\,}\mathrm{eV}$. In order to resolve the gap that is multiple orders of magnitude smaller than the bandwidth, we use $M=$5000$$ moments in the calculation. We use $R=5$ random vectors and 12 disorder realizations each. This increases the time cost, but not the memory cost of the algorithm.

V.0.0.1 We find an accurate transition point.

We perform finite size collapse of the data Caio et al. (2019) (see Appendix F) and find the transition point at $x_{c}=$0.28(3), with critical exponent $\nu=$0.9(6) accurately describing the transition, see Fig. 2. This result improves significantly from the virtual crystal approximation (VCA) result Lent et al. (1986); Dziawa et al. (2012); Yan et al. (2014), and ab initio Gao and Daw (2008), but is consistent with the most precise experimental data Tanaka et al. (2013); Zhong et al. (2015). Other experimental estimates are made with only two values of the concentration $x$ Dimmock et al. (1966); Yan et al. (2014); Phuphachong et al. (2017).

V.0.0.2 With the 6-orbital model we find the phase diagram of a broader class of materials.

To study a larger parameter space that includes other possible alloys of the \ceX_1-xSn_xTe family that manifest the mirror Chern phase, we use the simplified 6-orbital model. Besides the composition $x$, we also vary the onsite energy of the dopant cation $\ce{X}$, approximating compounds with lighter or heavier ions and similar electronic structures. Figure 3 shows the phase diagram. We find that the phase boundary differs from the VCA method, where the topological index only depends on the average onsite energy of the cations.

V.0.0.3 We implemented the first practical approach to 3D topological invariants for alloys.

Conclusions. — Our method is the first to allow the computation of topological invariants of realistic 3D alloys. Disorder in the crystalline structure is present in naturally found and artificially grown compounds, and it is inherent to substitutional alloys. However, a computationally efficient method to analyze topological properties of realistic disordered materials was missing.

VI.0.0.1 Our method makes accurate prediction on the phase transition concentration.

We apply our method to study the critical concentration of \cePb_1-xSn_xTe, and find an estimate that agrees with the most precise experimental data 222In the case of Ref. Yan et al. (2014), the ARPES data for $x=0.3$ is not conclusive, but it is interpreted with the help of VCA theory as a trivial phase. We claim the that for $x=0.3$, the experimental data and the theoretical results presented in Yan et al. (2014) are not sufficient to ensure that $x=0.3$ lies on the trivial phase.. Whereas it is not possible to reconcile any of the other theoretical studies (known to us) with all the experimental data, specially Ref. Zhong et al. (2015).

VI.0.0.2 Our method has the best scaling.

The scaling of computational time with system size of our method is better than the scaling of other methods available in the literature. By using the kernel polynomial method, we achieve a computational time scaling of $\xi^{d+1}$ with the localization length. Since we do not use eigenvalue solvers, only matrix-vector multiplication, the memory requirement only scales linearly with the sample volume as $\xi^{d}$.

VI.0.0.3 Our method is general and should allow automated search for topological alloys.

Beyond Chern numbers, our formalism allows calculation of all $\mathbb{Z}$ valued strong and weak topological invariants in all dimensions Perez-Piskunow et al. . This method makes the automated discovery of topological alloys feasible, and can guide synthesis of new alloys in the future. In this study we use energetically accurate tight-binding models obtained from ab initio calculations performed on pure materials as input. Using these tight-binding amplitudes for the atoms and bonds that appear in the alloy, we generate large disordered samples with various concentrations. We are able to probe the topology, something that would not be accessible with other methods. Simulation of small clusters with various disorder configurations is feasible using ab initio methods Gao and Daw (2008), tight-binding parameters obtained for all local environments would serve as a more accurate input for disordered models Prodan and Kohn (2005).

VI.0.0.4 Future research in $\mathbb{Z}$ invariants, open questions on how to compute $Z_{2}$, HOTI or topological semimetals.

This work opens several directions for future research. Our method is directly applicable to all types of disorder, as well as quasicrystalline and amorphous systems in symmetry classes that admit the topological marker formalism Perez-Piskunow et al. . This approach is not restricted to electrons in solids, and can be combined with finite element methods to analyze topology in disordered classical mechanical and photonic systems Stützer et al. (2018); Mitchell et al. (2018); Daraio et al. (2018); Zhou et al. (2018). We expect our method to perform well in disordered time-reversal breaking Weyl-semimetals with nonzero Hall conductivity, and further refinements could extend it to the time-reversal invariant case. While we are not aware of a similar formulation of $\mathbb{Z}_{2}$ indices, KPM could be utilized to calculate quantized responses associated with these phases, such as the quantized magnetoelectric effect of 3D strong topological insulators Bianco and Resta (2013). A similar approach could also be applied to higher order topological insulators to calculate multipole moments of the charge density Benalcazar et al. (2017). We expect that KPM could be used to study a wide variety of related topics in condensed matter physics, such as probing localization, topological Anderson insulators, or numerical renormalization group studies of the topological markers.

Acknowledgements.

Author contributions. — The project to study topological invariants using KPM was initiated by P. Perez-Piskunow and D. Varjas, the scope of the project was later refined using contributions from all authors. P. Perez-Piskunow and D. Varjas wrote the code used for the numerical calculations and performed calculations. D. Varjas performed the large-scale numerical calculations on the \cePb_1-xSn_xTe tight-binding models. All authors took part in the analysis of the method and writing the manuscript. Acknowledgments. — We are grateful to M. D. Caio for discussions and insights on local topological markers, M. Hoskam for help with implementing the tight-binding models, and B. Nijholt for technical assistance with the numerical calculations. We thank S. Roche and A. Cummings for critical revision of the manuscript, and P. Brouwer for helpful discussions. This work was supported by ERC Starting Grant 638760, the Netherlands Organisation for Scientific Research (NWO/OCW) as part of the Frontiers of Nanoscience program, NWO VIDI grant 680-47-53, the US Office of Naval Research, University of Chicago MRSEC through a Kadanoff-Rice postdoctoral fellowship, and the European Union’s Horizon 2020 research and innovation programme under grant agreement No 824140.

Appendix A Spectral projector operator expanded with the kernel polynomial method

The band projector is expanded with the kernel polynomial method (KPM) Weiße et al. (2006) and it is a finite range approximation, where the range depends linearly on the number of moments used in the expansion. Each moment of the expansion of the projector operator applied to a vector is obtained recursively, by applying the Hamiltonian to expanded vector of the previous iteration. The recursive algorithm effectively spreads a local vector to the neighboring sites via the hopping terms of the Hamiltonian, as depicted in Fig. 4.

Concurrently, each order of the expansion increases the precision in energy of the approximation, therefore, setting an energy resolution for the expansion. To apply KPM, operators 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/\Delta$, the full bandwidth divided by the mobility gap. This can be a large ratio even if $\xi$ is small, for example in an atomic insulator with a large range of on-site energies and vanishing hoppings. Typically $W/\Delta$ and $\xi$ increase together near a phase transition where the gap closes, but the bandwidth changes slowly, and does not affect the scaling of the computational cost with $\xi$.

The kernel polynomial method provides a stable and efficient method to expand the action of any function of an operator $\hat{f}$ that depends on the Hamiltonian $H$ and a set of parameters $\lambda$, on a vector $\ket{v}$ Weiße et al. (2006); García et al. (2015). The expansion up to order $M$ is

[TABLE]

The coefficients, called moments in the context of KPM expansions, are defined as

[TABLE]

and the vectors $\ket{v_{m}}$ satisfy the recursion relation

[TABLE]

We approximate the projector operator defined as the step function

[TABLE]

and in this case, the coefficients take the form

[TABLE]

Equipped with the KPM expanded projector, we proceed to evaluate matrix elements of topological markers. These are finite polynomials of $\hat{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 $E_{F}$, only a small number of sparse matrices and dense vectors are stored at any given time. Most of the time cost comes from sparse matrix-vector multiplications, linear in the system size. The number of operations is proportional to the number of moments $M$.

Appendix B Scaling of stochastic trace

We optimize the calculation further by utilizing the stochastic trace approximation to evaluate the trace. We take $R$ independent random phase vectors $\ket{r_{i}}$ that are only nonzero inside the region $S$, $\left\langle\bm{x},l|r_{i}\right\rangle=\delta_{\bm{x}\in S}\exp(\mathrm{i}\phi_{\bm{x},l,i})$ with $\phi_{\bm{x},l,i}\in[0,2\pi]$ independent random phases for all sites and orbitals. The trace of an operator $\hat{\mathcal{O}}$ equals the expectation value

[TABLE]

where $\mathbbm{E}$ denotes the expectation value over random vector realizations and we introduced the notation $\operatorname{Tr}_{\rm st}\hat{\mathcal{O}}$ for the random variable giving the stochastic trace of operator $\hat{\mathcal{O}}$. The above equality is proved by using that the random phases are independent, hence $\mathbbm{E}\left(e^{\mathrm{i}(\phi_{\bm{x},l}-\phi_{\bm{x^{\prime}},l^{\prime}})}\right)=\delta_{\bm{x},\bm{x}^{\prime}}\delta_{l,l^{\prime}}$ and only the diagonal entries contribute to the expectation value.

The standard deviation of stochastic trace of an operator scales with the total square magnitude of the off-diagonal entries which enter in the expectation value with random phases Weiße et al. (2006):

[TABLE]

where $\sigma(X)=\sqrt{\mathbbm{E}(|X|^{2})-|\mathbbm{E}(X)|^{2}}$ is the standard deviation and we used that $\sigma(e^{\mathrm{i}\phi_{\bm{x},l}})=1$. We also use that the standard deviation of the sum of independent random variables obeys $\sigma\left(\sum_{i}X_{i}\right)=\sqrt{\sum_{i}\sigma\left(X_{i}\right)^{2}}$.

We are concerned with the stochastic trace of topological markers, such as the Chern and mirror Chern operators. In order to draw conclusions, we need to know the scaling of the off-diagonal matrix elements with respect to the relevant length scales in the problem. There are three length scales, the lattice constant $a$, the localization length $\xi$ and the system 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 other two distances, and, as explained in the main text, we are interested in systems whose size is proportional to the localization length, so we will set $L=c\xi$ in the end. We introduce $a$ as the lattice constant here for clarity, but it is an arbitrary reference length scale we can define in fully disordered (e.g. amorphous) systems as well, for example as the typical spacing of sites. As it cancels from the final result, this argument does not rely on the assumption of an underlying regular lattice.

We start with the Chern operator in 2D

[TABLE]

where the $a^{-2}$ prefactor is included to measure all distances in units of $a$, this way the sum of the diagonal entries of $\hat{C}$ on a site coincides with the Chern number in a clean system. We numerically verify (see Fig. 5) that the off-diagonal matrix elements scale as

[TABLE]

where $f$ is a dimensionless function and we suppressed the dependence on the internal degrees of freedom. $f$ is a quickly decaying function for $(x-x^{\prime})/\xi\gg 1$ in insulating systems, as matrix elements of $\hat{P}$ also decay at the length scale of $\xi$.

The Chern marker averaged over a square region $S$ of size $L$ around the origin is given by

[TABLE]

where we still measure length in units of $a$. Substituting (10) we find for the standard deviation of $C$

[TABLE]

In the second line we took the limit of $\xi\gg a$, so we can replace sums over sites with integrals as $\sum_{x\in[-L/2,L/2]}=1/a\int_{-L/2}^{L/2}dx$. In the third line we changed the integration variables to $\tilde{\bm{r}}=\bm{r}/\xi$. The result is only dependent on the ratio of the localization length and the system size $c=L/\xi$ and the number of random vectors $R$, but not $\xi$. As the integral is proportional to the integration area, $c^{2}$ for $c\gg 1$, the overall scaling of the error with $c$ is $1/\sqrt{Rc^{2}}$.

We find a similar scaling for other topological markers, such as the 3D winding number in chiral classes Song and Prodan (2014); Mondragon-Shem et al. (2014). In general, the standard deviation of the stochastic trace evaluation of the invariant depends only on the ratio of the system size and the localization length:

[TABLE]

For the mirror Chern operator

[TABLE]

the scaling is

[TABLE]

where $\xi_{z}$ is the localization length in the $z$ direction. This form is justified by the fact, that the contributions to the mirror Chern number are centered on the invariant planes, but are spread out on layers in a thickness proportional to $\xi_{z}$, see Fig. 6. The total for all layers is, however, constant, hence the $a/\xi_{z}$ prefactor. This is the key difference compared 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^{\prime}$. The mirror Chern marker averaged over a square region of size $L$ is given by

[TABLE]

and we find using a similar derivation for the standard deviation (setting $\xi=\xi_{z}=L/c$)

[TABLE]

which is only dependent on the ratio of the localization length and the system size $c$ and the number of random vectors $R$.

In the numerical calculations we split the stochastic trace in 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^{\prime}$, resulting in a constant factor reduction in the error. Splitting the stochastic trace into more regions (e.g. separate for each layer parallel to the mirror plane) results in further reduction in the error, at the cost of increased computational effort. The overall scaling of the computational time with $\xi$ for a fixed standard deviation is the same up to a constant factor for all of these schemes, $\xi^{d+1}$.

Appendix C Geometry used in the numerics

As described in Section III.0.0.6 of the main text, we build a tight-binding model with PBC using translation vectors $W[1,1,0]$, $L_{1\overline{1}0}[1,-1,0]$ and $L_{z}[0,0,1]$. This geometry preserves the reflection symmetry with $[1,1,0]$ normal and contains $W\times L_{1\overline{1}0}\times L_{z}$ unit cells with 36 degrees of freedom each. The averaging region of the stochastic trace extends the full width of the system in the $[1,1,0]$ direction and contains half of the linear size in the perpendicular directions, as depicted by the inner box in Fig. 7. Imposing PBC in all directions eliminates gapless surface states, and the (mobility) gap guarantees that the Fermi projector is short-ranged.

However, imposing PBC in the direction normal to the mirror planes results in two mirror invariant planes. As argued in the main text, in a sample with open boundary conditions in the other directions, this results in a doubling of the interface modes, around the edges of the mirror invariant planes. On the other hand, $C_{M}$ counts helical modes, while the 3D mirror Chern number is given by the number of chiral modes, this factor of $1/2$ cancels the previous factor of $2$. We conclude that Eq. (6) is applicable to this 3D $M_{z}$ symmetric geometry with PBC in $z$.

Appendix D 18-orbital tight-binding model

D.0.0.1 specify 18-orbital model

We use the 18-orbital model of \ceSnTe and \cePbTe derived in Lent et al. (1986). The cubic rock-salt structure has two sublattices A and C (referring to the anion and cation nature of the atoms occupying them), the first occupied by \ceTe and the second by \ceSn or \cePb atoms. Each site hosts spinful $s$, $p$ and $d$ orbitals, 18 degrees of freedom in total, with annihilation operators $\bm{c}_{l}$ ($l=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 integrals $H_{ll^{\prime}m\bm{d}}$ in the linear combination of atomic orbitals (LCAO) methodSlater and Koster (1954), where $l$ and $l^{\prime}$ is the total angular momentum of the orbitals connected on the two sites and $m$ is the angular momentum of the bonding along the bonding axis $\bm{d}$ ($m=0,\;1,\;2$ for $\sigma$, $\pi$, $\delta$ bonding respectively). The matrices $H_{ll^{\prime}m\bm{d}}$ are $2(2l+1)\times 2(2l^{\prime}+1)$ and are proportional to the identity in spin space. The onsite terms contain different onsite energies for the various orbitals $E_{l}$ and ${\bm{L}\cdot\bm{S}}$ SOC terms with strength $\lambda_{l}$. The tight-binding Hamiltonian reads:

[TABLE]

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 (\ceSn or \cePb) is chosen randomly with probability $1-x$ and $x$. The value of $E_{l\bm{r}}$ is assigned accordingly to be $E_{lc}^{\text{SnTe}}$ and $E_{lc}^{\text{PbTe}}$ respectively. The superscripts \ceSnTe and \cePbTe refer to the two sets of parameters for the two pure materials. If $\bm{r}\in A$, we use a weighted average $E_{l\bm{r}}=\left[nE_{la}^{\text{SnTe}}+(6-n)E_{la}^{\text{PbTe}}\right]/6$ where $n$ is the number of nearest neighbor sites occupied by \ceSn atoms. The second term is the ${\bm{L}\cdot\bm{S}}$ spin-orbit coupling, $\bm{L}_{l}$ is the vector of angular momentum-$l$ operators ([math] for $l=0$). The values of $\lambda_{l\bm{r}}$ 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 $\bm{r}\in A$ and $\bm{r}^{\prime}\in C$. Depending on the atoms at sites $\bm{r}$ and $\bm{r}^{\prime}$ the value of $V_{l,l^{\prime},m,\bm{r},\bm{r}^{\prime}}$ is set to $V_{l,l^{\prime},m}^{\text{SnTe}}$ if one of the sites is \ceSn or $V_{l,l^{\prime},m}^{\text{PbTe}}$ if one of the sites is \cePb.

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 1.

Because of the identical outer shell electronic structure of \ceSn and \cePb, the alloy composition $x$ does not affect the doping level, therefore, we set the Fermi level $E_{F}$ to ensure half filling for all compositions, see Fig. 1 on the main text.

Appendix E 6-orbital tight-binding model

E.0.0.1 specify 6-orbital model

We adopt the 6-orbital model of \ceSnTe and \cePbTe originally described in Mitchell and Wallis (1966), and used in Refs. Hsieh et al., 2012; Fulga et al., 2016; Schindler et al., 2018; Sessi et al., 2016. Each site hosts spinful $p$-orbitals, 6 degrees of freedom in total with a vector of annihilation operators $\bm{c}$. This Hamiltonian is formally identical to (18) but only includes $p$-orbitals and $pp\sigma$ hopping, while the hopping range is extended to second neighbors. The tight-binding Hamiltonian reads:

[TABLE]

The first term is the onsite energy (also termed “mass term”), and it is the main source of disorder in our simulation. $m_{\bm{r}}$ takes the value of $m_{\ce{Te}}$ on the $A$ sublattice, while for the $B$ sublattice a value is chosen between $m_{\ce{Sn}}$ and $m_{\ce{X}}$ with probability $1-x$ and $x$. The second term is the ${\bm{L}\cdot\bm{S}}$ spin-orbit coupling, its value depends on the sublattice only (identical for \ceSn and \ceX). The third term is a $pp\sigma$ type of hopping Slater and Koster (1954) that only connects $p$-orbitals oriented along the direction of the bond ${\hat{\bm{d}}_{\bm{r},\bm{r}^{\prime}}}$.

We include first neighbor $[001]$ and second neighbor $[110]$ hoppings, with amplitudes that depend on the sublattices. We summarize the parameter values used for numerical results in Table 2.

Because of the identical outer shell electronic structure of \ceSn and \cePb, the alloy composition $x$ does not affect the doping level, therefore, we set the Fermi level $E_{F}$ to ensure half filling for all compositions.

Appendix F Evaluation and uncertainty of the finite-size scaling parameters

We start from the data $C_{M}=f_{i}(x)$ of the the mirror Chern number as a function of the concentration of \cePb_1-xSn_xTe where $i=1,\dots,N$ correspond to different system sizes $L_{i}$. This data is obtained by averaging multiple disorder realizations. We wish to perform a finite-size scaling collapse of the data. The changes of variables $\tilde{x}_{i}=(x-x_{c})L_{i}^{1/\nu}$ define new functions $\tilde{x}\mapsto\tilde{f}_{i}(\tilde{x})$, where both the critical concentration $x_{c}$ and the correlation length exponent $\nu$ are parameters to estimate so that the curves for different $L$ collapse onto a universal master curve. Hence, we seek to minimize the distance between $\tilde{f}_{i}(\tilde{x})$ over their common domain of definition. To do so, we interpolate the discrete data to perform the change of variable, evaluate $\tilde{f}_{i}(\tilde{x})$ over a fixed interval, and compute the variance of the data. This defines a cost function $\Phi(x_{c},\nu)$ that is minimal at optimal parameters $(x_{c}^{*},\nu^{*})\simeq($0.28$,$0.9$)$. The cost function is plotted as a function of the parameters $(x_{c},\nu)$ in figure 9. To evaluate the uncertainty on the optimal parameters, we evaluate the cost function $\Phi(x_{c}^{*},\nu^{*})$ at the optimal parameters $(x_{c}^{*},\nu^{*})$ for different disorder realizations [that give different values of the initial data $C_{\text{M},L}(x)$]. The standard deviation gives an estimation of the uncertainty $u(\Phi)$ on the cost function. The domain $D=\left\{(x_{c},\nu)\mid\Phi(x_{c},\nu)\leq\Phi(x_{c}^{*},\nu^{*})+u(\Phi)\right\}$ where the cost function is closer to its minimum than the uncertainty $u(\Phi)$ is approximately an ellipse (see figure 9) from which the uncertainties $u(x_{c})\simeq$0.03 and $u(\nu)\simeq$0.6 can be directly estimated.