Many-impurity scattering on the surface of a topological insulator

We theoretically address the impact of a random distribution of non-magnetic impurities on the electron states formed at the surface of a topological insulator. The interaction of electrons with the impurities is accounted for by a separable pseudo-potential method that allows us to obtain closed expressions for the density of states. Spectral properties of surface states are assessed by means of the Green’s function averaged over disorder realisations. For comparison purposes, the configurationally averaged Green’s function is calculated by means of two different self-consistent methods, namely the self-consistent Born approximation (SCBA) and the coherent potential approximation (CPA). The latter is often regarded as the best single-site theory for the study of the spectral properties of disordered systems. However, although a large number of works employ the SCBA for the analysis of many-impurity scattering on the surface of a topological insulator, CPA studies of the same problem are scarce in the literature. In this work, we find that the SCBA overestimates the impact of the random distribution of impurities on the spectral properties of surface states compared to the CPA predictions. The difference is more pronounced when increasing the magnitude of the disorder.

Since the pioneering work of Anderson on the absence of diffusion in random lattices 1 , different models of disorder have played a major role in understanding optical and transport properties of real solids with point defects. The advent of two-dimensional (2D) Dirac materials, such as the surface of topological insulators, graphene and carbon nanotubes, has brought renewed interest in low-dimensional disordered systems. One of the most salient features of Dirac materials is the appearance of a gapless energy spectrum that depends linearly on momentum (Dirac cones). This dispersion makes electrons behave as massless fermions with a Fermi velocity much lower than the speed of light. The single-parameter hypothesis of disordered systems, introduced by Abrahams et al. 2 , led to the general belief that all electron states were exponentially localized in 2D systems. Although this prediction works nicely when the energy spectrum depends quadratically on momentum, it turns out that extended states may arise in 2D systems with linear dispersion where quasi-particles undergo a localisation-delocalisation transition by varying the magnitude of disorder 3 . Therefore, it becomes apparent that electron dynamics in disordered 2D Dirac materials may substantially differ from what is known in conventional solids.
Single-particle spectral properties of disordered systems, such as the density of states (DOS), can be assessed by means of the Green's function averaged over disorder realisations 4,5 . In general, the configurationally averaged Green's function cannot be calculated exactly and various approximations of different degree of sophistication are employed. Among them, self-consistent methods stand out because they correctly explain the main features of the DOS as inferred from photo-emission and soft X-ray experiments 6 .
Impurities and other point defects are common sources of disorder in 2D Dirac materials [7][8][9][10] . Electron scattering by impurities yields spectral features, such as circular s-wave resonances, that can be targeted by scanning tunnelling experiments (see reference 8 and references therein). Theoretical treatments of many impurities are often based on the self-consistent Born approximation (SCBA) [11][12][13][14][15][16][17][18][19] . In the SCBA the imaginary part of the self-energy in the bare Green function is replaced by the self-energy of the full Green function, that renders the problem self-consistent. If the impurity potential is assumed short ranged, the SCBA leads to particularly simple expressions for the average Green's function, from which the DOS is readily determined. The so-called coherent potential approximation (CPA) represents another example of a self-consistent approach routinely used for the theoretical analysis of conventional disordered matter 4,20 . However, CPA studies of spectral and transport properties of disordered 2D Dirac materials are still scarce in the literature [21][22][23] , particularly in the context of surface states of topological insulators.
In this work we study many-impurity electron scattering on the surface of a topological insulator by means of self-consistent methods, namely SCBA and CPA, with the aim of comparing their predictions. The analysis and conclusions can be trivially extended to any 2D material where electron dynamics can be described by the massless Dirac equation. The interaction of the electron with the scatterers is accounted for by a separable pseudo-potential model [24][25][26][27][28][29][30][31] . In spite of its seemingly more complicated form, the separable pseudo-potential model is amenable to analytical solution and allows us to obtain closed expressions for the average Green's function within the SCBA and CPA frameworks. In particular, short-range potentials approaching the δ-function limit, frequently used in previous works 12,18,[32][33][34] , can be viewed as limiting cases of the separable pseudo-potential model. We will show that the SCBA average Green's function can also be obtained from the CPA calculations in the limit of diluted impurities and small magnitude of disorder. However, the main conclusion of this work is that the SCBA overestimates the impact of point-like scatterers on the spectral properties, compared to the CPA predictions. The discrepancy becomes greater when increasing the impurity concentration.

Results
Theoretical model. The Hamiltonian operator of an electron in a pristine surface of a topological insulator will be denoted as H 0 . It is diagonal in the basis of plane waves where 35 Here v is a matrix element having dimensions of velocity (for example, v = 4 × 10 5 m/s for Bi 2 Te 3 36 , and v = 5 × 10 5 m/s for Bi 2 Se 3 37 ), σ x and σ y are Pauli matrices and k = (k x , k y ) is the in-plane momentum. The corresponding bands are simply given as E k = ± v|k| (Dirac cones). Notice that we are restricting the study to a single Dirac cone and neglecting intervalley scattering as well as the contribution of non-linear terms in momentum. Moreover, we are dealing with the case of disorder that modifies the spectrum but does not collapse the bulk gap. Therefore, our results are valid as long as the 2D effective model holds.
Let us address how the electron interacts with impurities located at the surface of the topological insulator. We will assume that they are placed on a regular square lattice of parameter a. Notice that a is not related to the size of the unit cell of the crystal structure of the topological insulator. In fact, electrons do not see the crystal structure since we are using a continuous approximation for the Hamiltonian (1). We will focus on binary disorder hereafter. To this end, two different species of impurities A and B are considered. A given site of the square lattice is occupied by an impurity A with probability c or by an impurity B with probability 1 − c . Therefore, the separable pseudo-potential operator can be cast in the form 25,29 The index n runs over all sites R n of the square lattice and ω(r − R n ) = � r | ω n � will be referred to as shape function. n is the coupling constant that takes on two values A and B at random, with probability c and 1 − c respectively. Hence, the probability distribution in this model of binary disorder is The electron Hamiltonian in the presence of the impurities is the sum of the Hamiltonian H 0 corresponding to the translationally invariant system and the random part V , namely H = H 0 + V . The retarded Green's function operators (resolvents) corresponding to H and H 0 are where z = E + i0 + . Notice that G 0 (z) is diagonal in the basis of plane waves with by virtue of Eq. (1). We will concern ourselves with the ensemble average G(z) av of the Green's function operator in the random medium. The subscript ' av ' indicates the average over the probability distribution (3).
The knowledge of G(z) av allows us to obtain the spectral properties of an electron on the surface of a topological insulator scattered off by a random array of impurities. In general, the average Green's function operator cannot be obtained exactly and some approximations are needed. The conceptually simplest way of finding an approximation to G(z) av is by introducing an effective, translationally invariant medium represented by a Green's function operator G eff (z) such that G eff (z) = � G(z)� av . The first level of approximation is reached in the case of very weak scattering by assuming that the array of impurities is periodic with a coupling constant given as the following average www.nature.com/scientificreports/ This approach is known as the Virtual Crystal Approximation (VCA) (see, e.g., reference 4 ). The VCA is a reasonably good description only if c → 0 (or equivalently c → 1 ) and A ≃ B . However, more elaborated, self-consistent methods have a much wider range of validity. Within these methods, the VCA appears usually as a first constant term in the expansion of the Green's function, as described in the following sections.
Once the effective Green's function is obtained, important physical quantities can be calculated. In particular, the average DOS per unit area is easily computed by the following expression We will take S = 1 and referred to ρ(E) as the DOS hereafter.
Self-consistent Born approximation. In this section we consider the effects of disorder within the SCBA. In the framework of this approximation, the Green's function operator of the effective medium is taken as where the self-energy operator � SCBA (z) is diagonal in the basis of plane waves The self-energy � SCBA (k, z) is to be determined self-consistently from the following equation Here C(k − k ′ ) is the disorder correlator that depends on the transferred momentum. In the case of the separable pseudo-potential model (2) we get where ω(k) = d 2 r e ik·r ω(r) is the Fourier transform of the shape function and = A − B is the magnitude of disorder.
For convenience, we define the self-energy as the product of an effective coupling constant SCBA (z) and the shape function ω(k) as follows Therefore, Eq. (8c) is written as Notice that Eq. (10b) is valid for any shape function and consequently it is suitable for the study of finiterange impurity potentials. However, particularly simple expressions are found for point-like impurities, namely when ω(k) becomes independent of k . Since the resulting integral is divergent at large momenta, we impose a momentum cutoff k c (or, equivalently, we introduce a finite bandwidth) and set where θ is the Heaviside step function and the impurity lattice constant a is introduced for convenience. The assumption of point impurities, albeit crude, is widely used in the context of continuous models. It is worth stressing, however, that the pseudo-potential model introduced in this work does not require this assumption and more structured electron-impurity potential can be handled within the same footing. In particular, lattice distortion around the impurity site can be accounted for by setting a finite-range function ω(r).
Since H 0 (k) in the numerator of (5b) is an odd function of momentum, the corresponding integration vanishes. Moreover, we find it more convenient to express the results in terms of the coupling constant obtained within the VCA (6) by defining � SCBA (z) = SCBA (z) − VCA with z = z − VCA . Hence, the self-consistent equation for SCBA can be expressed in a compact form as follows where we have defined www.nature.com/scientificreports/ Coherent potential approximation. The CPA traces back to the sixties and has proven to be a successful mean field theory for the study of various elementary excitations (electrons, phonons, excitons, magnons) in disordered systems [38][39][40][41] . The CPA combines two basic ideas. On one side, the average Green's function of the disordered system is calculated by introducing a periodic (translationally invariant) effective medium. On the other hand, this effective medium is determined by demanding that the fluctuations of the Green's function average out to zero, thus leading to a self-consistency condition 5 . In the single-site CPA combined with the separable pseudo-potential model, the electron motion in the effective medium is represented by the following Hamiltonian 25,26 where CPA (z) is in general complex and will be determined self-consistently from the condition In contrast to H = H 0 + V with V given by Eq. (2), the effective Hamiltonian H eff has the full symmetry of the impurity lattice since CPA (z) is taken to be independent of the site. The difference between both Hamiltonians can be expressed as To proceed, we consider the t-matrix operator associated with a single site 4 It can be proven that the requirement G eff (z) = � G(z)� av yields the well-known CPA condition [4][5][6] Therefore, from (15) we finally get where, within the one-band approximation (see "Methods" for more details), we have It is worth mentioning that � ω n | G eff (z) | ω n � becomes site independent since the effective medium is translationally invariant. Thus, the ensemble average in the case of binary disorder (3) poses no problem and (17a) leads to The above expression is valid for any shape function. In particular, in the case of point-like impurities (11) one gets where F (z) is defined in (12b). Once more, it is more convenient to express the left-hand side of Eq. (18) in terms of the coupling constant obtained within the VCA (6) by defining whence Notice that, expanding the CPA self-consistent equation given by (20), we can get the SCBA. This can be obtained by solving for � CPA (z) and expanding the result in a Taylor series for small c and up to third order In fact, the SCBA can be obtained as a truncation of the series of the CPA. This is further clarified in the diagrammatic formalism with Feynman rules. The SCBA takes into account the two irreducible diagrams shown in Fig. 1a for the self-energy. The first diagram is the constant VCA term while the second one describes the double scattering off by a single impurity with a dressed internal propagator. On the other hand, CPA sums all the diagrams with any number of scattering events on the same impurity that, upon a proper re-summation 42 , gives the self-consistent equation (20) [see Fig. 1b]. where refers either to SCBA or CPA , α(Ē) is real and Ŵ(Ē) > 0 corresponds to a disorder-induced broadening. Notice that we express the results as a function of Ē that is the shifted energy after taking the limit of z →Ē + i0 + . By analysing the symmetry properties of the CPA self-consistent equation (20), we find that it is invariant under the exchange (α, Ŵ,Ē, c) → (−α, Ŵ, −Ē, 1 − c) and (�, c) → (−�, 1 − c) . Therefore, we can restrict ourselves to � > 0 and 0 ≤ c < 0.5 since all the other scenarios can be obtained from the former range of parameters. For the SCBA self-consistent equation (12a), different symmetries are obtained due to the truncation of the series expansion of the CPA (see Fig. 1). On the one hand, once high-order terms in c are neglected ( c → 0 ), the symmetry (�, c) → (−�, 1 − c) is lost. Notice that in order to investigate the range c ≥ 0.5 , the expansion in Eq. (21) must be performed around 1 − c instead of c. Hence, the expression (12a) can be used only for 0 ≤ c < 0.5 . On the other hand, due to the truncation shown in Eq. (21), an artificial symmetry in energy is generated in the SCBA self-energy, resulting in a symmetric DOS. This symmetry, already reported in the literature 12,43,44 , is due to the presence of a single type of non trivial diagram, as explained in more detail in the diagrammatic approach shown in Fig. 1. In CPA, odd terms in are considered in the expansion of the self-energy, leading to an asymmetric DOS. This asymmetry is consistent with purely numerical findings in Dirac-like systems 45,46 .
In fact, apart from the constant VCA term, the SCBA depends on a single parameter related to disorder, namely c 2 , while the CPA needs both c and separately. We define the SCBA disorder parameter β as follows Hence, expressing the energies in units of v/a , we can write Eq. (12a) as a function of a single dimensionless disorder parameter given by Eq. (23) and the energy cut-off E c ≡ vk c where it is understood that � SCBA = � SCBA (Ē) . This equation is invariant under the exchange (α, Ŵ,Ē) → (−α, Ŵ, −Ē) . Hence, the real part of the effective coupling constant is an odd function of the shifted energy, α(Ē) = −α(−Ē) , while the imaginary part is an even function, Ŵ(Ē) = Ŵ(−Ē) . With these considerations in mind, the SCBA equation (24) can be solved explicitly in the case of Ē = 0 , finding the zero-energy solution This exponential behaviour is opposite to the case of single-node Weyl semimetals studied in reference 18 , where a critical point signals a disorder-induced phase transition. In the Dirac-like Hamiltonian (1b), no critical behaviour is observed as a function of the magnitude of disorder, as discussed later.
The SCBA self-consistent equation can be solved analytically for energies |Ē| ≪ E c and small disorder |� SCBA (Ē)| ≪ E c (see "Methods" for details). In this regime, the coupling constant can be approximated as where W(x) is the Lambert-W function 47 . Figure 2 shows a comparison of the analytic expression for the coupling constant with the numerically solved SCBA equation as a function of energy for weak disorder. Notice that the approximated solution agrees exceedingly well with the numerics, as long as the range of parameters considered www.nature.com/scientificreports/ fulfils all the conditions for the approximation to be valid. For the sake of completeness, the CPA results are plotted in solid lines as well. It is worth mentioning that the CPA and the SCBA results coincide very nicely for small disorder, as predicted by the series-expansion interpretation of the SCBA introduced in the previous section. Most importantly, upon increasing the magnitude of the disorder and the energy, the SCBA tends to overestimate the impact of the impurities on the DOS. In the following, we analyse the limit of Ē → 0 . As already mentioned, the SCBA predicts an exponentiallike broadening given by Eq. (25) that resembles, for small disorder parameter, a purely exponential decay [see "Methods" for further details] The above expression of the broadening allows us to write explicitly the value of the DOS in the zero-energy limit within the SCBA [see the DOS expression in "Methods" for further details] Figure 3 shows a comparison of the analytic limit and the results of the numerically solved Ŵ(Ē) and DOS as a function of the disorder parameter β within the SCBA. The absence of a disorder-induced phase transition is patent and it is a crucial difference between the 2D and 3D cases. In 3D Weyl semimetals, the phase transition is firmly established from analytical and numerical methods [48][49][50][51][52] , while in the case of 2D Dirac-like materials, the phase transition is not observed. In graphene, Dirac-like approaches and tight-biding approximations lead to the exponential behaviour shown in Eq. (27) reported in the literature 43,44,53 . On the other hand, 2D models with mass terms such as the Bernevig-Hughes-Zhang (BHZ) model show more complex phase diagrams [54][55][56] .
The absence of a phase transition is observed in the CPA results as well. In fact, we find a smooth dependence of the DOS at zero energy Ē on the disorder magnitude and the fraction c of A impurities, as seen in Fig. 4. Notice that in the CPA both parameters are needed and they can not be combined into a single disorder   Figure 5 shows in more detail the range of equivalence of both approximations. For small c and ( c 0.2 and 0.5 in the figure) the SCBA and CPA coincide whereas for higher values of disorder the overestimation of the SCBA becomes noticeable. Notice that, in the range of weak disorder, the analytic limit given by Eq. (28) is accurate and the DOS follows the exponential trend −1/ ln ρ(Ē) ∼ c� 2 . In a wider range of magnitude of disorder and concentration c, the disagreement becomes apparent, leading to an excess of the DOS of the order of the value itself, as seen in Fig. 4. Moreover, the SCBA predicts a threshold for non-zero DOS smaller than the one predicted by the CPA, as shown in the Fig. 4d, where the DOS is plotted as a function of .
For non-zero energy, the tendency remains the same and the SCBA results in an overvaluation of the effect of the impurities. Due to the strictly non-zero DOS for |Ē| > 0 , we can compute the relative error defined as Figure 6 shows the relative error at a given energy as a function of and c. We observe that the discordance grows with the energy, as previously shown in (see Fig. 2).
We conclude by stressing the range of validity of the CPA. The CPA has been proven to reliably obtain the self-energy for a wide range of scenarios. It yields the correct result in the weak scattering limit (where it coincides www.nature.com/scientificreports/ with the SCBA), in the strong limit and in the dilute limit 42,57 . In fact, the only approximation assumed in the CPA condition is that if the averaged single-site t-matrix is zero [Eq. (16)], then the averaged T-matrix of the whole system is zero. This approximation is correct whenever the spatial correlation of disorder is negligible. The single-site CPA incorrectly treats multiple scattering terms associated with clusters of fixed number of neighbour sites 5 . Diagrammatically, it corresponds to the fact that the self-energy in CPA does not include wigwam diagrams with crossing lines, whose contribution is negligible as long as the scattering length of the impurity potential is smaller than a 57 . Therefore, if the impurities are diluted and short-range order is absent, the results of the CPA are essentially exact. Finally, let us stress the validity of the results obtained in this work for the understanding of other 2D Dirac materials. After a trivial rotation, the electron Hamiltonian (1a) is basically the same that of a low-energy electron in graphene. Hence, our results are of interest in the description of graphene impurities 46 specially in the non-magnetic impurities case. Starting from the seminal work by Noro et al. 43 , graphene disordered sheets have been studied extensively within the SCBA approach 14,44 , showing a sizeable effect of the disorder present in the samples. The numerical findings also show the behaviour presented here for the DOS 58 . Moreover, proposals have been made in order to obtain the averaged DOS of those systems by measuring the quantum capacitance 59 .

Discussion
We have solved the effective medium approximation for a many-impurity scattering problem on a 2D surface of a topological insulator within the SCBA and CPA. Moreover, we have analysed in detail the differences, weaknesses and strengths of both methods. The simplicity of the SCBA allows us to extend the analytic calculations, bringing almost exact analytic results for small magnitude of disorder and low concentration of impurities without the need for the numerical solution of the self-consistency conditions. On the other hand, the CPA enables us to exactly solve the problem for any number of single-impurity scattering events, yielding reliable results even in the non-pertubative limit. Moreover, as expected by the correspondence of SCBA and CPA for weak disorder, both approximations coincide in the range of dilute and weakly-interacting impurities.
A reliable determination of the effective coupling constant, or equivalently, the self-energy, is of central importance since it allows us to calculate all physically meaningful quantities. Aiming to achieve this, it is crucial to use the appropriate method matching the regime of concentration and disorder strength properly. In conclusion, our finding thus not only calls for a revision of current theories based on the SCBA, but also provides a reliable implementation of the (more accurate) CPA for studying impurity scattering of 2D Dirac matter.

Methods
One-band approximation. Starting from Eq. (13), the Green's function operators associated to H eff and H 0 satisfy 5 We now take into account the closure relation of the plane waves being the identity operator, and Eq. (5a) to obtain where the index K runs over the vectors of the reciprocal lattice of the impurity lattice. In the one-band approximation, the Fourier transform of the shape function is assumed to vanish outside the Brillouin zone 25,26 . In this way, we only retain the term K = 0 in the expansion (32). Therefore where � CPA (k, z) = CPA (z)|ω(k)| 2 /a 2 .
Using the closure relation (31) we get where S is the area of the system. After converting the sum over k into an integration we finally obtain (17b).
Calculation of the coupling constant in the SCBA. As mentioned in the text, the SCBA self-consistent equation can be solved exactly at Ē = 0 and approximately in the case of weak disorder. In the case of Ē = 0 , considering the symmetry properties of Eq. (24) in the main text, we find that the real part of the coupling constant must be zero. Therefore, replacing � → −iŴ , we conclude that the self-consistent condition reduces to whose solution is given by Eq. (25).
In the weak disorder regime and for energies Ē ≪ E c , the coupling constant fulfils |� SCBA |(Ē) ≪ E c . Therefore, we can expand the SCBA equation as where � SCBA = � SCBA (Ē) . Considering solutions with Im (� SCBA ) < 0 , we obtain Eq. (26). This approximate solution resembles the exact case at Ē = 0 for small β . In fact, at zero energy, we obtain Eq. (27), which corresponds to the first term in the series expansion of Eq. (25) for β ≪ 1.
Expression for the DOS. The DOS per unit area is obtained from the Green's function using Eq. (7). In the case of the 2D effective Hamiltonian we are dealing with, this expression is written as After some algebra and expressing the energy in units of v F /a and the coupling constant �(Ē) as given by Eq. (22), when ω(k) = aθ(k c − k) we obtain the following expression. For the sake of simplify, hereafter we omit the dependence on Ē in α ≡ α(Ē) and Ŵ ≡ Ŵ(Ē) .