Ab-initio Theory of Fourier-transformed Quasiparticle Interference Maps and Application to the Topological Insulator Bi$_2$Te$_3$

The quasiparticle interference (QPI) technique is a powerful tool that allows to uncover the structure and properties of electronic structure of a material combined with scattering properties of defects at surfaces. Recently this technique has been pivotal in proving the unique properties of the surface state of topological insulators which manifests itself in the absence of backscattering. In this work we derive a Green function based formalism for the ab initio computation of Fourier-transformed QPI images. We show the efficiency of our new implementation at the examples of QPI that forms around magnetic and non-magnetic defects at the Bi$_2$Te$_3$ surface. This method allows a deepened understanding of the scattering properties of topologically protected electrons off defects and can be a useful tool in the study of quantum materials in the future.


I. MOTIVATION
The scanning tunneling microscopy (STM) experiments of Crommie, Lutz and Eigler [1] and Hasegawa and Avouris [2], revealing standing density waves of the Cu(111) and Au(111) surfacestate electrons near defects, have pioneered a very powerful, direct method of imaging the surface electron liquid of metals. Together with Scanning Tunneling Spectroscopy (STS), the method gives unique insight on the quasiparticle interference (QPI), scattering phase shifts, and lifetime.
Especially when augmented by the Fourier-transformed quasiparticle interference map, as proposed by Petersen et al. [3], the method unveils scattering properties of quasiparticles off surface defects, giving information on the scattering vectors among points of the band structure. In this context, Fourier-transformed QPI maps provided one of the first experimental proofs of the existence of topological insulators [4], because it revealed the absence of intensity at back-scattering vectors, just as predicted by theory.
From a theoretical point of view, the calculation of QPI maps has been largely based on model methods, e.g. on topological insulator surfaces [5], where the surface band structure can be approximated by simple model Hamiltonians. In general, however, density-functional based methods are necessitated for a realistic description of the surface electronic structure and in particular of the impurity potential, where the charge relaxation around the impurity plays a major role in the correct description of the scattering phase shifts. A difficulty in density-functional calculations is that the density oscillations induced by the defect are very long-ranged, reaching tens or even hundreds of nanometers, so that supercell methods cannot practically reach this limit. These challenges can only be met by ab-initio Green function embedding methods, like the Korringa-Kohn-Rostoker (KKR) method.
As an example of an application, we refer to the calculations by Lounis et al. [6] of QPI on Cu(111) and Cu(001) surfaces due to an isolated impurity buried under the surface. These results show that ab-initio calculations of QPI maps over quite large surface areas are feasible with Green function techniques. However, in the case of Fourier-transformed QPI map, it is practical to express the result directly by a convolution of Green functions [7], avoiding the intermediate step of calculating the real-space map in a large surface area.
In this paper, we approach this problem and give applications in the field of topological insulators. In Sec. II we outline the formalism for real-space and Fourier-transformed QPI maps within the KKR method. Furthermore, we discuss the Fourier-transformed QPI for the practical case of multiple impurities and argue that the many-impurity problem is well approximated by the single-impurity result. We also discuss the extended joint density of states approach (exJDOS). In Sec. III we apply our formalism on the topological insulator Bi 2 Te 3 with surface impurities. This is implemented in the JuKKR code package [8]. Finally, we conclude with a summary in Sec. IV.

II. FORMALISM
Within the Tersoff-Hamann approximation [9], the STM differential tunneling conductivity at a bias voltage U is related to the space-and energy-resolved density of states n(r; E) at energy E and at the position r of the tip: n(r; E F + eU ) ∝ dI dU (U ) (E F is the Fermi level). In a QPI experiment we are interested in the difference of density induced by an impurity with respect to the pristine host surface, ∆n(r; E) = n imp (r; E) − n host (r; E). (1)

A. Green function and T -matrix approach
The Green function and T -matrix approach is a well-established method of calculating the density of systems with impurities. It has been applied to the QPI problem in ab-initio and model calculations, e.g. in Refs. 5, 7, 10. Here we give the formalism in an explicit real-space representation, because it forms the basis of the formalism and discussion in subsequent sections.
The difference in the local density of states is connected to the one-electron Green function of the system with impurity, G imp (r, r ; E), and to the one of the pristine host, G host (r, r ; E), by the well-known identity The trace is implied with respect to spin indices; in the case of a relativistic formalism with Dirac four-vector states, the trace includes also the large and small components. By virtue of the Dyson in terms of the T -matrix, ∆V = V imp − V host is the difference in the potential between the impurity and host systems. The T -matrix has the advantage that it is confined in the small region where the potential difference does not vanish and has to be calculated only once, irrespective of the range of r in the QPI calculation. For the calculation of G host , the translational invariance of the surface allows us to use the Bloch theorem. Decomposing the position vector as r = R + x, where R is a lattice translation vector parallel to the surface plane, while x is a vector in the primitive cell, we write where the Fourier-transform of the Green function obeys the spectral representation Here, Ψ αk is the host wavefunction, α is the band index, Ω rec = (2π) 2 /Ω cryst , with Ω cryst the total crystal surface area, and i0 represents an infinitesimal imaginary energy. The difference in Green functions, Equation (3), takes then the form The sum over R , R and the integration over d 3 x , d 3 x is confined to the sites where the T -matrix (and the impurity perturbation ∆V ) is non-vanishing. The former expression, Equation (7), includes phase factors e −ik·R e ik ·R for the inter-lattice-site propagation of the host Green function, when the impurity spreads over many lattice sites. The latter expression (8) is a more compact form where the matrix elements T αkα k (E) = (Ψ αk , T (E)Ψ α k ) were introduced (note that the summation includes all states, not just the ones at energy E). It leads to the stationary phase approximation [6] pinning the energy to the energy-shell E, if the observation point is far from the impurity (|R| → ∞).
For the Fourier-transformed QPI we need the Fourier transformation of the Green function along a surface parallel to, and at vertical distance z from, the crystal surface: In the step from (9) to (10) we used Equation (7) and we removed lattice sum R by virtue of (Ω BZ is the surface Brillouin zone area). The result represents the convolution of two Green functions, as expected from the Fourier transform of their products. Employing Equation (2), we arrive at the following expression for the Fourier transformed QPI: where it is implied that the complex conjugation operation ∆G * is done after the Fourier transformation. This result can easily be generalized for the spin density with σ∆G in the place of ∆G (σ is the vector of Pauli matrices).
The strongest density change measured by the STM is induced directly "above the impurity," i.e., at a vertical distance z > 0 from the position where ∆V = 0. This region is often excluded in from the Fourier transformation in experiment [11], otherwise the image is dominated by the transform of the impurity shape [12], while one seeks the scattering vectors. Additionally, the region close to the impurity is also excluded sometimes, because it may produce spurious background effects in the experiment [13,14]. In the calculation, the contribution of the excluded region (indicated by Ω excl ) must be subtracted explicitly, because the form (10) already includes a summation over all lattice sites. Thus we define ∆n(z; q; E) = ∆n(z; q; E) − Ω excl d 2 r n(z; q; E) e −iq·r (12) = ∆n(z; q; E) + 1 π ImTr However, since Ω excl is finite-sized, the integration is straightforward in real space.

B. Expression in the KKR formalism
In the KKR method, the Green function is expanded in site-dependent scattering wavefunctions at sites n. The vacuum is also described in a site-centered way by a continuation of the lattice structure beyond the surface, with the corresponding "empty sites" containing no atoms but a finite electron density. We denote the general position by r = X n + x = R i + χ µ + x, where the combined index n = (i; µ) defines a site X n by the lattice-vector R i and the sub-lattice vector χ µ , and where x is a position vector in the atomic site with respect to the site center. We employ the regular, R n L (x; E), and irregular, H n L (x; E), solutions of the scattering problem of the potential in the vicinity of the site R n , where L comprises angular momentum and spin indices of the incoming wave. R n L and H n L are (2 × 1) column-vectors in Pauli-Schrödinger theory and (4×1) column vectors in Dirac theory. Also the corresponding left-hand side solutions are needed, denoted byR n L (x; E) andH n L (x; E), respectively, that are row-vectors. The expansion breaks the Green function down into a single-site term and a multiple scattering term, where κ = (2mE) 1/2 / in Pauli-Schrödinger theory and κ = (2mE + E 2 /c 2 ) 1/2 / in Dirac theory. The single-site term describes the Green function of the potential at site n embedded in free space. It only depends on the local potential and its contribution to ∆G vanishes outside the impurity region. The second term describes the multiple scattering over all sites, expressed by the structural Green function coefficients G nn LL (E). These form a matrix G(E) that obeys an algebraic Dyson equation. This reads for the host system where we have expressed everything in reciprocal space. g(k; E) contains the structural Green function coefficients of free space and t host (E) is a site-diagonal (k-independent) matrix con- The analogon of the T -matrix in the KKR method is the scattering path operator of the impurity with respect to the host, τ (E) . It is expressed in terms of the single-site T -matrices of impurity and host, ∆t n LL (E) = t imp;n LL (E) − t host;n LL (E), and of the host structural Green function, by the Dyson-type equation τ = ∆t + ∆t G host τ . It is not site-diagonal, and has nonvanishing elements τ nn LL (E) only between sites (n, n ) for which ∆t n = 0 and ∆t n = 0. The structural Green function of the system with impurity is then expressed by in analogy to Equation (3).
Since the vacuum region is geometrically described by layers of empty sites parallel to the crystal surface, it is convenient to approximate the Fourier integration over a surface at distance z by an integration over a vacuum layer of volume Ω scan , centered at z: d 2 r → n∈Ωscan n d 3 r.

Expression (10) then becomes
is the difference of the single-site part of the Green function between the impurity and the host system, and is taken only in the impurity region Ω imp (it vanishes outside). In the above expression, the terms x; E) are 2 × 2 or 4 × 4 matrices (depending if the Pauli-Schrödinger or the Dirac theory is used) and must be traced to form the density [see Equation (2)].
Conveniently, they show no k-dependence and thus must be calculated only once at each energy; the same is true for the matrix elements of the scattering path operator, τ iµ;i µ L L (E). The only quantities that need to be calculated for a dense set of k-points (which implies a large numerical effort) are the host structural Green functions (Equation 15). Fortunately, by virtue of the principal layer and decimation techniques [15][16][17], the latter can be computed with a numerical effort that grows linearly with the number of atomic layers in the film, making possible the accurate simulation of the QPI in thick films (of the order of hundreds of atomic layers, if necessary) or semi-infinite geometries.
If we wish to calculate the quantity ∆n(z; q; E) (Equation 13), i.e., exclude the impurity and its immediate surroundings (indicated by Ω excl in Equation 13) from the Fourier transformation, then Equation (17) changes. The single site term ∆G S , vanishes automatically outside Ω imp . However, we must also explicitly subtract the contribution of the multiple-scattering term in Ω excl . The result is given by replacing ∆G S (q; E) by the following correction to the multiple scattering part For the calculation of the integrals we expand the wavefunctions in spherical harmonics, as is normally done in the KKR method [18], and we do the same for the exponential by the identity where Y lm are spherical harmonics. Thus the integral is decomposed in a spherical and an angular part. The integrals containing irregular functions that contribute to ∆G SL (q; E) are handled in an analogous way.
In summary, in the KKR method we calculate the quantities ∆n(q; E) (Equation 11) and ∆n(q; E) (Equation 13) by: C. Multiple scattering among impurities The experimental QPI Fourier transform is usually performed over a large surface area comprising many impurities. A direct simulation of this experiment should account for the contribution of the multiple-scattering events between impurities to the QPI. However, this is numerically expensive, since it involves the calculation of a large T -matrix, corresponding to the collection of all impurities in a large supercell, and perhaps even a statistical average over many impurity configurations. Fortunately, the Fourier-transformed QPI of a single impurity is an excellent approximation to the result of a random impurity distribution, simplifying the calculations. Fang et al. [19] have shown this approximation to hold to lowest order in the potential difference ∆V , i.e., in the Born approximation to the scattering amplitude. Here we argue that the approximation holds in general, allowing for the treatment of strong, e.g., resonant, scattering.
First we discuss the form of the multiple-scattering T -matrix. Let t n (E) be the T -matrix of a single impurity at site n and expressed in a matrix form in a localized basis set. Then, the full T -matrix of a collection of impurities obeys the expansion [20] T nn = t n δ nn + t n G host where the matrixG nm = G host nm (1 − δ nm ) contains the site-off-diagonal part of the Green function (a proof is given in the Appendix). The above expression includes all multiple-scattering events among impurities, while avoiding sequential scattering off the same impurity (since t n contains the sequential site-diagonal scattering to all orders of ∆V ).
We assume that the impurities are non-overlapping (which is a reasonable approximation at low concentration) and identical and are thus characterized by the same matrix t n = t ∀n. We also apply in part the stationary phase approximation for the host Green function [6], which is valid at long distances. This is justified at low impurity concentrations since the largest part of the surface, where the Fourier transform is performed, is covered by host atoms and is far from the impurities. Within this approximation, the host Green function may be approximated by G host (R n + x, R n + x ; E) ≈ K nn (x, x ; E) e ik nn ·R nn , where k nn is a stationary point on the constant energy surface E k nn = E, and is defined by the property that the group velocity v k nn must be parallel to the vector R nn = R n − R n . The quantity K nn (x, x ; E) contains the rest of the Green function, including a power-law decay with distance (K nn ∝ |R nn | −1/2 in two dimensions). The important consequence of this approximation for our purposes is that the phase of the long-distance propagation, e ik nn ·R nn , is governed only by the stationary point (in the case of multiple stationary points, a summation over the corresponding contributions is implied). Then we argue that, in the Fourier transform of Equation (10), the contribution of the first term of the rhs of Equation (22) is dominant and equal to the single-impurity contribution, while the remainder (the contribution of tG T ) is negligible.
We decompose the Green function difference [Equation (3)] in two terms corresponding to the decomposition of the T -matrix (22): All first-order terms, G host nm t G host mn , give identical contributions to the Fourier transform ∆G(q), because both G host mn and G host nm depend on the sites n and m only via the difference R nm . If N is the number of impurities, and setting one impurity at position m = 0, we have , which can be shown by changing the summation over R m to R nm . This is, however, not true for the higherorder terms. Applying the stationary phase approximation to G host n n in Equation (24) (last term), we obtain ∆G (2) nn ≈ m G host nm t n n G mn T n n K n n e ik n n ·R n n . In the Fourier transformation, the translational symmetry of the host allows us again to place m = 0 by a re-indexing, but the contribution of the last phase, e ik n n ·R n n , cannot be lifted. Since the impurities are randomly placed, the total contribution of the random phases over all impurities practically cancels in the A frequently used approach to the Fourier-transformed QPI is the Joint Density of States (JDOS) [4,7,21,22] or extended JDOS (exJDOS) [23] approach, which is applied if the constantenergy contours {k|E k = E} at energy E are known (e.g. from calculations or from angularresolved photoemission experiments), but the full Green function G host or the full T -matrix are not known. Motivated by Equation (10), one defines the quantity which a weighted convolution of the spectral amplitude at the pristine crystalline surface, proximation M k,k−q = |T k,k | 2 ∝ 1 + cos(s k , s k ) has been proposed [4], where s k is the spin polarization vector of the state at k. Additionally, a factor γ STM k,k−q = 1 − cos(v k , v k−q ) was introduced in Ref. [23] in order to promote standing wave formation by back-scattering (opposite group velocities), in the spirit of the stationary phase approximation. At the end, exJDOS(q; E) is expected to approximately reproduce |∆n(q; E)|, since both should peak at the scattering vectors.
The JDOS approach has been introduced as an ad-hoc model. The question is if it also constitutes an approximation to the theory expressed by Eqs. (7), (10), and (11). We find that exJDOS(q; E) ∝ |∆n(q; E)| 2 , under a number of assumptions that are scrutinized in the following. From Eqs. (10) and (11) we have (dropping the variables E and z) |∆n(q)| 2 = ∆n(q) ∆n(−q) We employ the stationary phase approximation to the Green function, according to which = kk a k aku k (r) T k,k u † k (r) e i(k−k)·r /|r| (28) at large distances |r| from the impurity (which is placed at r = 0). We omitted the band index α in order to simplify the notation. The discrete summation for a given r in (27) runs over k-points that are stationary with respect to the host Green function phase, i.e., they are pinned with respect to the energy, E k = Ek = E, and also pinned at such positions on the constant energy contour, that the group velocity v k is in the direction r and the group velocity vk is in the opposite direction, −r. The symbols u k and uk stand for the lattice-periodic part of the wave-function, while k and k run on the constant-energy contour. The last expression (28) is a convenient abbreviation with obvious shorthand notation for a k and ak. The Fourier transformation reads In the last step we converted the integration variable from d 2 r to rdrdθ and subsequently dθ to  (26), products of wavefunctions of the type u † k u k expressing the density, as well as products of T -matrix elements expressing the transition rate, will appear. In order to comply with the JDOS (25), the mixed-k (i.e., k = k ) density terms must be dropped. The terms a k and ak should be also ignored (or set to a constant). Finally, the weighting factor γ STM should be set in the definition of the exJDOS (25) in order to account for the stationary phase approximation.
The above discussion shows that the JDOS and exJDOS approaches are not quantitative approximations but qualitative models. Still, they comprise essential parts of the information that one usually seeks in Fourier-transformed QPI spectra and therefore constitute a useful tool for their analysis.

III. APPLICATIONS
To showcase the use of our newly developed method we apply it to the topological insulator Bi 2 Te 3 , which hosts nontrivial surface states characterized by spin-momentum locking that are protected by time-reversal symmetry against backscattering.
In our density functional calculations within the relativistic full-potential KKR Green function framework [27][28][29][30] we considered a 6 quintuple layer thick film of Bi 2 Te 3 using the experimental lattice constant [31], which was chosen such that the "top" and "bottom" surface states of the thin film decouple. We used an angular momentum cutoff of max = 3 including corrections for the exact shape of the cells [32,33] and the local spin density approximation [34] for the exchangecorrelation functional. The Fermi level was set such that it resides inside the bulk band gap, which ensures that the Fermi surface consists of the topological surface state alone without projections of bulk bands. Such a situation can be achieved experimentally by doping or gating. Figure 1 summarizes the setup of the calculation, where the bandstructure (a) and Fermi surface of the system at hand (b) are shown together with the spin-polarization of the topological surface state.
Single substitutional impurities were embedded into the Bi 2 Te 3 host at the outermost Bi site as indicated in Figure 1(c). We considered non-magnetic Te Bi (the subscript indicates the substituted position) defect, which occurs naturally as an anti-site defect, as well as magnetic Mn Bi and Fe Bi impurity atoms. It has been previously shown [24][25][26] that the substitutional Bi site is a thermodynamically stable position of transition metals in Bi 2 Te 3 , indicating the relevance of our findings. The impurities were embedded into the host system self-consistently making use of the Dyson equation [28] while neglecting structural relaxations around the impurities. The first shell of nearest neighbours was included in the calculation for a correct screening of the charge of the impurities. The resulting density of states for the three defects are shown in Figure 1(d). While the density of states at the Fermi level is small for the Te Bi defect, the incompletely filled d-shell of the transition metal impurities results in a higher density of states at the Fermi level. In accordance to Hund's rule, we find that the Mn Bi defect is close to half-filling whereas the Fermi level bisects the d-state resonance of the Fe Bi impurity. From analyzing the impurity density of states and in the spirit of the Friedel sum rule, which was recently demonstrated to hold in this class of materials [35], one expects considerable differences in the scattering properties of the different impurities.
In addition, the magnetic nature of the Mn and Fe defects is expected to reopen the forbidden backscattering channel due to breaking of time-reversal symmetry.
The strong hexagonal warping in Bi 2 Te 3 leads to a snowflake-like shape of the Fermi surface that enables two major scattering channels [12,36], q 1 and q 2 , which are depicted in Figure 1(e). The backscattering channel (q 1 ) is expected to be suppressed by time-reversal symmetry for scattering off non-magnetic defects and can be re-opened by scattering off magnetic defects.
The trivial scattering channel q 2 is however always possible. Next we compare the different signatures (intensities at q 1 , q 2 )in the FT-QPI image among the different impurities. Scattering off the non-magnetic Te Bi defect is characterized by a strong focus in forward (i.e., small-angle or q ≈ 0) scattering direction. This is a consequence of the forbidden backscattering that disallows the scattering vector q 1 as seen by looking at the scattering rate [27,29,30], which is illustrated for one particular incoming wave, k 0 , in Figure 2(j). The strong focus to small angle scattering leads to the absence of a signal at q 1 and only a moderate intensity at q 2 compared to the dominant star-like feature around q = 0. In contrast, the intensity at q 2 is stronger for the Mn Bi and Fe Bi impurities [see Figure 2(h,i)]. This is a consequence of the lesser focused scattering in forward direction as exemplified in Figure 2 In summary, our simulations of the impurity specific FT-QPI reveals that not only information on the host's electronic structure can be extracted but also the scattering properties of different impurities are accessible.

B. Comparison to joint density of states approaches
Usually the interpretation of experimental FT-QPI images is done by comparison to calculations based on the joint density of states. The comparison between the (S)JDOS in Figure 1(f,g) on the one hand and the impurity specific results of the FT-QPI presented in Figure 2 on the other hand reveals that a proper description of the impurity scattering is crucial for quantitative understanding of the scattering process off defects. In particular the intensity at the backscattering signal (q 1 ) is overestimated in the simple (S)JDOS approaches and the strong focus in small angle scattering is underestimated. A significant improvement was the exJDOS [see Equation (25)], that does include the correct scattering information of the different impurities and which is shown for the three defects (Te Bi , Mn Bi and Fe Bi ) in Figure 3. Qualitatively the correct FT-QPI can be reproduced and a strong suppression of the q 1 signal is found in accordance to the results of Figure 2.
This gives an a posteriori justification of the use of the exJDOS model in the interpretation of experimental QPI images and allows for the extraction of scattering information as the interplay between the electronic structure of the host system and the embedded impurities. It can, however, not be excluded that in other systems, the terms that are dropped in the JDOS-approaches (e.g., mixed k contributions or anisotropic scattering rate) become important.

IV. SUMMARY
We have derived a Green function and T -matrix based formalism for the calculation of Fourier-   tions because of a cancellation of the multiple-scattering wavefunction phase.
We have applied our KKR-based implementation to non-magnetic and magnetic defects embedded in the surface of the topological insulator Bi 2 Te 3 , providing microscopic insights into the scattering properties of topological surface state electrons and their response to time-reversal conserving or breaking defects.
We rewrite the first form as Let us denote the impurity sites with the index n. It is convenient to use a notation where we collect the single-site t-matrices in a site-diagonal matrix T d with (T d ) nn = t n δ nn . In analogy, we collect the diagonal part of the Green function in the matrix (G host d ) nn = G host nn δ nn . ∆V is trivially site-diagonal. Then, for the site-diagonal parts we use the second form of Equation (A.1) that yields T d = ∆V + T d G host d ∆V , i.e., Eliminating ∆V in (A.2) and (A.3) we obtain which is equivalent to (22).