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

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.

Introduction. -Topological materials have attracted continuing interest from both the fundamental physics and the material science communities for the last decade [1,2]. The program to theoretically classify non-interacting crystalline insulators has been completed, tabulating possible topological phases in all space groups [3][4][5]. Recent efforts focus on automated high-throughput methods to discover and classify new topological materials [6][7][8], culminating in the production of comprehensive databases of topological insulators and semimetals [8][9][10][11][12].
However, not all topological insulators are compounds with perfect stoichiometry. The first three-dimensional topological insulator to be predicted [13] and experimentally realized [14] was an alloy-Bi x Sb 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 [15,16]. This approach ignores the intrinsic disorder in alloys, and it is insufficient to explain topological transitions that appear at strong disorder [17,18], or accurately find critical concentrations.
The topological invariant converges to its bulk value in samples larger than the localization length ξ. This is the main limitation in resolving topological phase transitions as ξ diverges. Therefore, the asymptotic scaling of the computational cost with ξ is the main distinction between different numerical approaches.
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 [19][20][21][22][23][24][25][26][27]. However, these methods involve solving at least one eigenvalue equation with size equal to the number of degrees of freedom, resulting in the complexity ξ 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 [28], with a complexity scaling as ξ 3(d−1) , allowing for maximum sizes of 5 × 10 5 degrees of freedom.
We present an algorithm to efficiently identify topological phases of strongly disordered systems using the kernel polynomial method (KPM) [29][30][31], 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 [22] 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 [29] from a small number of random vectors, combined with KPM, is suited for this task. With ξ d+1 scaling of the computational effort, it is the most efficient method in three dimensions.
As a concrete example, we apply our method to lead tin telluride Pb 1 -x Sn x Te 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 Pb 1 -x Sn x Te [32][33][34], we find a critical concentration that matches the one found experimentally [35][36][37][38][39][40][41].
Review of existing algorithms. -We focus on the vicinity of disorder-driven phase transitions, where the localization length ξ diverges. The finite but potentially large value of ξ defines the relevant length scale of our problem. Two regions that are closer than ξ feel each other's presence. Therefore, with open boundary conditions, the bulk has to be further than ξ to the edge; analogously, with periodic boundary conditions when the Figure 1. Top: surface spectra of a 20 × 80 × 80 sample of Pb1-xSnxTe in the trivial (left) and topological (right) phase. The presence of a gapless surface Dirac cone indicates the mirror Chern insulator phase. Bottom: transition between trivial and mirror Chern phase when varying x for Pb1-xSnxTe calculated using our method with various system sizes. Inset: Finite size collapse of the curves with xc 0.28(3) and ν 0.9 (6).
linear size of the system is smaller than ξ 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 ξ. Fluctuations of local quantities resulting from disorder scale with ξ, so averaging over a larger sample provides a good approximation of the thermodynamic limit.
The momentum-space formalism of topological invariants applies to disordered systems by studying a periodic system with a large disordered supercell of volume ξ d . This is equivalent to taking a finite torus and threading fluxes through its cycles [20,42,43]. 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 [21], topological markers [22], pseudospectra [23], and noncommutative index theorems [24][25][26][27]. 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 ξ 3d , restricting them to small system sizes in three dimensions.
The scattering invariant formalism [44] 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 [45] and scales as ξ 3d−3 for d > 1.
General strategy. -Computing the exact spectral by full diagonalization of the Hamiltonian is numerically expensive. Instead, we approximate the projector using the kernel polynomial method with the Jackson kernel [29], 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 finiterange Hamiltonian. In an insulating system, the projector is a local operator with matrix elements x|P |x ∝ exp (−|x − x |/ξ) that have a decay length ξ. Hence, the error of the approximation scales as exp(−M/ξ), and the number of moments necessary for fixed precision scales linearly with the localization length as M ∼ ξ [46,47]. See Appendix A. We use the topological marker formalism introduced by Bianco and Resta [22]. All Z-valued topological markers are a partial trace per unit volume of a local operatorν where the sum runs over the sites x inside the subsystem S with volume |S|, and their internal degrees of freedom λ. The operatorν is a polynomial of the spectral projector, position and symmetry operators, such that ν 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 [22]. An example of topological marker is the real space expression for the Chern number [22]: Here, A is the area of the subsystem,x andŷ are the two components of the position operator, and [·, ·] is the commutator. Topological markers for all strong and weak Z-valued topological invariants have similar algebraic expressions of the projected position operators [2,48,49], making a straightforward application of our method to these cases [50]. We are not aware of similar formulations of Z 2 topological indices suitable for KPM. To estimate the trace per volume, we use the stochastic trace approximation [29] Tr where |r i are random phase vectors localized in the region S. The standard error of this approximation scales as ξ 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). We build a supercell of size L > ξ. To ensure that the invariant ν 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ν is localized with ξ < 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.
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 ∼ ξ and M ∼ ξ), results in an algorithm with a computational cost scaling of ξ d+1 .
Application to mirror Chern number. -Our method provides better scaling than existing approaches in d ≥ 3. We apply it to disordered 3D mirror Chern insulators. These are a widely studied class of topological crystalline materials with a Z topological classification that relies on reflection symmetry [51]. Several experimental realizations are known, including alloys [34,[37][38][39][40][41].
In a reflection-symmetric system of fermions, all wave functions are eigenstates of the mirror operatorM z , with eigenvalues ±i. The Chern numbers C ± for mirror-even and mirror-odd wave functions are wherex ± =M ±PxPM± , andỹ ± =M ±PŷPM± are the projected position operators restricted to the mirror-even or mirror-odd subspaces. HereM ± 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 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 oper-atorM z = i M + −M − commutes with the projector P and position operatorsx andŷ, we express the mirror Chern number as 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 system [52].
Tight binding model of Pb 1−x Sn x Te. -Topological crystalline insulators (TCI) protected by reflection symmetry [36] were theoretically predicted [51] and experimentally observed [34,[37][38][39][40][41] in Pb 1 -x Sn x Te alloys. They host metallic surface states on the surfaces that are symmetric with respect to the mirror plane [2,51,53]. Lead tin telluride was studied using either the virtual crystal approximation (VCA) [32,34] or ab initio methods [33], 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.
We consider the substitutional disorder of the lead tin telluride alloy Pb 1 -x Sn x Te coming from replacing some Pb ions for Sn ions. This disorder is nonmagnetic, and it preserves the reflection symmetry on average, which is sufficient to protect the gapless surface states [28,54]. We disregard other types of symmetry breaking disorder appearing naturally in Pb 1 -x Sn x Te [54], such as ferroelectric structural distortion [51,54,55], and magnetic dopants [51,54,56].
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 tightbinding parameters for both SnTe and PbTe derived from ab initio simulations [32].
To simulate the alloy, we substitute randomly Sn for Te with probability x. We incorporate substitutional disorder by using the hopping amplitudes of SnTe for Sn-Te bonds and PbTe amplitudes for Pb-Te bonds. The onsite parameters of Te atoms are slightly different in SnTe and PbTe; we therefore use a weighted average of these depending on the local environment. We also use the appropriate onsite terms, including L · S spin-orbit coupling (SOC), depending on the type of the Sn or Pb atom. Further details can be found in Appendix D.
When investigating the onsite energy dependent phase diagram of X 1 -x Sn x Te alloys, we use a simplified model that only includes 6 spinful p orbitals, with 12 bands [51,57,58]. We include L · 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 Sn and X sites. For more details, see Appendix E.
Results. -We define the Hamiltonians and perform the KPM expansions using the Kwant software package [59].
The code to reproduce the figures in this article is available in Ref. [60]. First, we study the topological phase transition in the realistic 18-orbital model of Pb 1 -x Sn x Te. We build a tight-binding model with PBC that preserves reflection symmetry and contains W × L 110 × 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 eV and band gap of less than 0.3 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.
We perform finite size collapse of the data [61] (see Appendix F) and find the transition point at x c = 0.28 (3), with critical exponent ν = 0.9(6) accurately describing the transition, see Fig. 2. This result improves significantly from the virtual crystal approximation (VCA) result [32,34,40], and ab initio [33], but is consistent with the most precise experimental data [39,41]  To study a larger parameter space that includes other possible alloys of the X 1 -x Sn x Te 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 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.
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.
We apply our method to study the critical concentration of Pb 1 -x Sn x Te, and find an estimate that agrees with the most precise experimental data [63]. Whereas it is not possible to reconcile any of the other theoretical studies (known to us) with all the experimental data, specially Ref. [41].
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 ξ 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 ξ d .
Beyond Chern numbers, our formalism allows calculation of all Z valued strong and weak topological invariants in all dimensions [50]. 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 [33], tight-binding parameters obtained for all local environments would serve as a more accurate input for disordered models [47].
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 [50]. 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 [64][65][66][67]. We expect our method to perform well in disordered time-reversal breaking Weylsemimetals 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 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 [68]. A similar approach could also be applied to higher order topological insulators to calculate multipole moments of the charge density [69]. 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.
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 Pb 1 -x Sn x Te tight-binding models. All authors took part in the analysis of the method and writing the manuscript.
The band projector is expanded with the kernel polynomial method (KPM) [29] 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/∆, the full bandwidth divided by the mobility gap. This can be a large ratio even if ξ is small, for example in an atomic insulator with a large range of on-site energies and vanishing hoppings. Typically W/∆ and ξ 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 ξ.
The kernel polynomial method provides a stable and efficient method to expand the action of any function of an operatorf that depends on the Hamiltonian H and a set of parameters λ, on a vector |v [29,30]. The expansion up to order M iŝ The coefficients, called moments in the context of KPM expansions, are defined as and the vectors |v m satisfy the recursion relation We approximate the projector operator defined as the step functionP Equipped with the KPM expanded projector, we proceed to evaluate matrix elements of topological markers.
These are finite polynomials ofP 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 |r i that are only nonzero inside the region S, x, l|r i = δ x∈S exp(iφ x,l,i ) with φ x,l,i ∈ [0, 2π] independent random phases for all sites and orbitals. The trace of an operator O equals the expectation value where E denotes the expectation value over random vector realizations and we introduced the notation Tr stÔ for the random variable giving the stochastic trace of operator O. The above equality is proved by using that the random phases are independent, hence E e i(φ x,l −φ x ,l ) = δ x,x δ l,l 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 offdiagonal entries which enter in the expectation value with random phases [29]: where σ(X) = E(|X| 2 ) − |E(X)| 2 is the standard deviation and we used that σ(e iφ x,l ) = 1. We also use that the standard deviation of the sum of independent random 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 ξ 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ξ 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 where the a −2 prefactor is included to measure all distances in units of a, this way the sum of the diagonal entries ofĈ 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 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 )/ξ 1 in insulating systems, as matrix elements ofP also decay at the length scale of ξ.
The Chern marker averaged over a square region S of size L around the origin is given by where we still measure length in units of a. Substituting (B4) we find for the standard deviation of C In the second line we took the limit of ξ a, so we can replace sums over sites with integrals as x∈[−L/2,L/2] = 1/a L/2 −L/2 dx. In the third line we changed the integration variables tor = r/ξ. The result is only dependent on the ratio of the localization length and the system size c = L/ξ and the number of random vectors R, but not ξ. As the integral is proportional to the integration area, c 2 for c 1, the overall scaling of the error with c is 1/ √ Rc 2 . We find a similar scaling for other topological markers, such as the 3D winding number in chiral classes [48,49]. 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: For the mirror Chern operator the scaling is where ξ 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 ξ z , see Fig. 6. The total for all layers is, however, constant, hence the a/ξ 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 . The mirror Chern marker averaged over a square region of size L is given by and we find using a similar derivation for the standard deviation (setting ξ = ξ z = L/c) 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 offdiagonal entries with z = −z , 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 ξ for a fixed standard deviation is the same up to a constant factor for all of these schemes, ξ d+1 . 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
We use the 18-orbital model of SnTe and PbTe derived in [32]. 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 Te and the second by Sn or Pb atoms. Each site hosts spinful s, p and d orbitals, 18 degrees of freedom in total, with annihilation operators 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 md in the linear combination of atomic orbitals (LCAO) method [70], where l and l 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 d (m = 0, 1, 2 for σ, π, δ bonding respectively). The matrices H ll md are 2(2l + 1) × 2(2l + 1) and are proportional to the identity in spin space. The onsite terms contain different onsite energies for the various orbitals E l and L · S SOC terms with strength λ l . The tight-binding Hamiltonian reads: 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 E lr is assigned accordingly to be E SnTe lc and E PbTe lc respectively. The superscripts SnTe and PbTe refer to the two sets of parameters for the two pure materials. If r ∈ A, we use a weighted average E lr = nE SnTe la + (6 − n)E PbTe la /6 where n is the number of nearest neighbor sites occupied by Sn atoms. The second term is the L · S spin-orbit coupling, L l is the vector of angular momentum-l operators (0 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 V l,l ,m,r,r is set to V SnTe l,l ,m if one of the sites is Sn or V PbTe l,l ,m if one of the sites is Pb. All of the onsite energies and hopping terms are spinindependent, 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 Sn and Pb, 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.   Table II. Tight-binding parameters in electronvolts for SnTe used in the 6-orbital model. We use the same Hamiltonian as Ref.
Sessi et al. [72], but the numerical values of the parameters differ due to different normalization and sign conventions. Logarithm of the cost function Φ(xc, ν) as a function of the parameters (xc, ν). The region D allowed by uncertainties around the minimum is delimited by a thick black dashed line.
We start from the data C M = f i (x) of the the mirror Chern number as a function of the concentration of Pb 1 -x Sn x Te where i = 1, . . . , 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 variablesx i = (x − x c )L 1/ν i define new functionsx →f i (x), where both the critical concentration x c and the correlation length exponent ν are parameters to estimate so that the curves for different L collapse onto a universal master curve. Hence, we seek to minimize the distance betweenf i (x) over their common domain of definition. To do so, we interpolate the discrete data to perform the change of variable, evaluatef i (x) over a fixed interval, and compute the variance of the data. This defines a cost function Φ(x c , ν) that is minimal at optimal parameters (x * c , ν * ) (0.28, 0.9). The cost function is plotted as a function of the parameters (x c , ν) in figure 9. To evaluate the uncertainty on the optimal parameters, we evaluate the cost function Φ(x * c , ν * ) at the optimal parameters (x * c , ν * ) for different disorder realizations [that give different values of the initial data C M,L (x)]. The standard deviation gives an estimation of the uncertainty u(Φ) on the cost function. The domain D = {(x c , ν) | Φ(x c , ν) ≤ Φ(x * c , ν * ) + u(Φ)} where the cost function is closer to its minimum than the uncertainty u(Φ) is approximately an ellipse (see figure 9) from which the uncertainties u(x c ) 0.03 and u(ν) 0.6 can be directly estimated.