Local and global patterns in quasiparticle interference: a reduced response function approach

TL;DR

This paper introduces a reduced response function approach to analyze quasiparticle interference patterns, linking local and global features to electronic structure geometry, topology, and topological invariants in physical systems.

Contribution

It develops a new RRF method for interpreting QPI patterns, connecting them to electronic topology and demonstrating robustness of topological indicators.

Findings

Global QPI patterns can indicate topological numbers.

The imaginary part of RRF is the generalized joint density of states.

Robustness of topological indicators against local QPI features.

Abstract

A physical system exposes to us in a real space, while its description often refers to its reciprocal momentum space. A connection between them can be established by exploring patterns of quasiparticles interference (QPI), which is experimentally accessible by Fourier transformation of the scanning tunneling spectroscopy (FT-STS). We here investigate how local and global features of QPI patterns are related to the geometry and topology of electronic structure in the considered physical system. A reduced response function (RRF) approach is developed that can analyze QPI patterns with clear physical pictures. It is justified that the generalized joint density of states, which is the imaginary part of RRF, for studying QPI. Moreover, we reveal that global patterns of QPI may be indicators of topological numbers for gapless systems, and demonstrate that robustness of such indicators against distractive local features of QPI for topological materials with complicated band structures.

Equations

$$n_{\alpha }(\mathbf{r},\omega )=-\frac{1}{\pi }\mathrm{Im}\{\mathrm{Tr}[{\sigma }_{\alpha }G(\mathbf{r},\mathbf{r},\omega )]\},$$

$$G({\mathbf{k}}^{\prime },\mathbf{k},\omega )=G_0(\mathbf{k},\omega ){\delta }_{\mathbf{k}{\mathbf{k}}^{\prime }}+G_0({\mathbf{k}}^{\prime },\omega )T_{{\mathbf{k}}^{\prime }\mathbf{k}}(\omega )G_0(\mathbf{k},\omega ),$$

$$\delta n_{\alpha \beta }(\mathbf{r},\omega )=-\frac{1}{\pi }\mathrm{Im}\{\mathrm{Tr}[\sum _{{\mathbf{k}}^{\prime }\mathbf{k}}{\sigma }_{\alpha }{G}_0({\mathbf{k}}^{\prime },\omega )V_{{\mathbf{k}}^{\prime },\mathbf{k}}^{\beta }{\sigma }_{\beta }{G}_0(\mathbf{k},\omega )]\}.$$

$$\delta n_{\alpha \beta }(\mathbf{p},\omega )$$

$${\Lambda }_{\alpha \beta }(\mathbf{p},\omega )=\sum _{\mathbf{k}}{V}_{\mathbf{k}+\mathbf{p},\mathbf{k}}^{\beta }\mathrm{Tr}\left[{\sigma }_{\alpha }{G}_0(\mathbf{k}+\mathbf{p},\omega ){\sigma }_{\beta }{G}_0(\mathbf{k},\omega )\right].$$

$${\Lambda }_{\alpha \beta }(\mathbf{p},\omega )=\sum _{\mathbf{k}s{s}^{\prime }}\frac{P_{\alpha \beta }{F}_{\alpha \beta }^{s{s}^{\prime }}(\mathbf{k}+\mathbf{p},\mathbf{k},\omega )}{({\omega }_{\mathbf{k}+\mathbf{p}}-s^{\prime }{d}_{\mathbf{k}+\mathbf{p}}+i{0}^{+})({\omega }_{\mathbf{k}}-s{d}_{\mathbf{k}}+i{0}^{+})},$$

$$P_{\alpha \beta }{F}_{\alpha \beta }^{s{s}^{\prime }}=\mathrm{Tr}\left[V_{\mathbf{k}+\mathbf{p},\mathbf{k}}^{\beta }(u_{\mathbf{k}+\mathbf{p},s^{\prime }}^{†}{\sigma }_{\alpha }{u}_{\mathbf{k},s})(u_{\mathbf{k},s}^{†}{\sigma }_{\beta }{u}_{\mathbf{k}+\mathbf{p},s^{\prime }})\right].$$

$$A_s(\mathbf{k},\omega )\equiv -\frac{1}{\pi }\mathrm{Im}[\frac{1}{{\omega }_{\mathbf{k}}-s{d}_{\mathbf{k}}+i{0}^{+}}]=\delta ({\omega }_{\mathbf{k}}-s{d}_{\mathbf{k}})$$

$$B_s(\mathbf{k},\omega )\equiv -\frac{1}{\pi }\mathrm{Re}[\frac{1}{{\omega }_{\mathbf{k}}-s{d}_{\mathbf{k}}+i{0}^{+}}],$$

$$\delta n_{\alpha \beta }(\mathbf{p},\omega )$$

$${\mathcal{R}}_{\alpha \beta }(\mathbf{p},\omega )={\mathcal{S}}_{\alpha \beta }(\mathbf{p},\omega )+i{\mathcal{J}}_{\alpha \beta }(\mathbf{p},\omega ),$$

$${\mathcal{S}}_{\alpha \beta }(\mathbf{p},\omega )$$

$${\mathcal{J}}_{\alpha \beta }(\mathbf{p},\omega )$$

$$n_{\alpha \beta }(\mathbf{p},\omega )\simeq \frac{1}{2}(P_{\alpha \beta }^{\ast }{\mathcal{R}}_{\alpha \beta }(\mathbf{p},\omega )+h.c.).$$

$$P_{\alpha \beta }{F}_{\alpha \beta }^{s{s}^{\prime }}(\mathbf{k}+\mathbf{p},\mathbf{k},\omega )$$

$$\sum _{s{s}^{\prime }}{r}_{AB}^{s{s}^{\prime }}(\mathbf{p},\omega )=$$

$$\sum _{\mathbf{k}s{s}^{\prime }}\left[A_s(\mathbf{k}+\mathbf{p},\omega )B_{s^{\prime }}(\mathbf{k},\omega )+A_s(\mathbf{k}+\mathbf{p},\omega )A_{s^{\prime }}(\mathbf{k},\omega )\right]$$

$$\sum _{s{s}^{\prime }}{r}_{AB}^{s{s}^{\prime }}(\mathbf{p},\omega )$$

$$=\sum _{s{s}^{\prime }}{\int }_{E_{\mathbf{k}}^s=\omega }\frac{1}{\omega -E_{\mathbf{k}+\mathbf{p}}^{s^{\prime }}+i{0}^{+}}\frac{1}{|{\nabla }_{\mathbf{k}}{E}_{\mathbf{k}}^s|}dl.$$

$$\omega -E_{{\mathbf{k}}_{\omega }^s(t_0)+\mathbf{p}}^{s^{\prime }}=0$$

$${\dot{E}}_{{\mathbf{k}}_{\omega }^s(t_0)+\mathbf{p}}^{s^{\prime }}={\frac{d{E}_{{\mathbf{k}}_{\omega }^s(t)+\mathbf{p}}^{s^{\prime }}}{dt}|}_{t=t_0}=0.$$

$$r_{AB}^{s{s}^{\prime }}(\mathbf{p}+\mathbf{\delta },\omega )\sim -\frac{1}{|{{\mathbf{v}}^{\mathbf{s}}}_{{\mathbf{k}}_0}||{{\mathbf{v}}^{{\mathbf{s}}^{\prime }}}_{{\mathbf{k}}_0+\mathbf{p}}|}\frac{1}{|\mathbf{\delta }\cdot {\mathbf{\varrho }}_{{\mathbf{k}}_0+\mathbf{p},{\mathbf{k}}_0}^{s{s}^{\prime }}{|}^{1/2}}\{\begin{array}{rlrl}& \text{sgn}({\mathbf{n}}_{{\mathbf{k}}_0+\mathbf{p}}\cdot {\mathbf{\varrho }}_{{\mathbf{k}}_0+\mathbf{p},{\mathbf{k}}_0}^{s{s}^{\prime }}),& \mathbf{\delta }\cdot {\mathbf{\varrho }}_{{\mathbf{k}}_0+\mathbf{p},{\mathbf{k}}_0}^{s{s}^{\prime }}>0& \\ & i,& \mathbf{\delta }\cdot {\mathbf{\varrho }}_{{\mathbf{k}}_0+\mathbf{p},{\mathbf{k}}_0}^{s{s}^{\prime }}<0& \end{array},$$

$${\mathbf{\varrho }}_{{\mathbf{k}}^{\prime },\mathbf{k}}^{s{s}^{\prime }}={\kappa }_{{\mathbf{k}}^{\prime }}^{s^{\prime }}{\mathbf{n}}_{{\mathbf{k}}^{\prime }}^{s^{\prime }}-{\kappa }_{\mathbf{k}}^s{\mathbf{n}}_{\mathbf{k}}^s.$$

$$J_{\alpha \beta }(\mathbf{p},\omega )$$

$$=\sum _{{\mathbf{k}}_{0}s{s}^{\prime }}{F}_{\alpha \beta }^{s{s}^{\prime }}(\mathbf{k}+\mathbf{p},\mathbf{k}){\left|\frac{\partial (E_{\mathbf{k}+\mathbf{p}}^{s^{\prime }},E_{\mathbf{k}}^s)}{\partial (k_x,k_y)}\right|}_{\mathbf{k}={\mathbf{k}}_0}^{-1},$$

$$=\sum _{{\mathbf{k}}_{0}s{s}^{\prime }}\frac{F_{\alpha \beta }^{s{s}^{\prime }}({\mathbf{k}}_0+\mathbf{p},{\mathbf{k}}_0)}{\left|{\mathbf{v}}_{{\mathbf{k}}_0+\mathbf{p}}^{s^{\prime }}\times {\mathbf{v}}_{{\mathbf{k}}_0}^s\right|},$$

$$E_{{\mathbf{k}}_0+\mathbf{p}}^{s^{\prime }}=\omega ,$$

$$E_{{\mathbf{k}}_0}^s=\omega .$$

$${\mathbf{v}}_{{\mathbf{k}}_0+\mathbf{p}}^{s^{\prime }}\times {\mathbf{v}}_{{\mathbf{k}}_0}^s=0,$$

$${\mathcal{H}}_{\mathbf{k}}=k^n(\cos n\theta {\sigma }_1+\sin n\theta {\sigma }_2)$$

$${\mathcal{H}}_{TI}(\mathbf{k})=v(k_x{\sigma }_y-k_y{\sigma }_x)+\lambda k^3\cos 3{\theta }_{\mathbf{k}}{\sigma }_z,$$

$$H_0(k_x,k_y)$$

$$+\Lambda (k)(3k_x^2-k_y^2)k_y{\sigma }_z,$$

Full text

I Introduction

Interference of quasiparticles in the presence of impurities leads to a modulation of the local density of states (LDOS). The LDOS is accessible experimentally with scanning tunneling spectroscopy (STS)Crommie et al. (1993); Hasegawa and Avouris (1993); Sprunger et al. (1997), and the pattern of modulation is further extracted by its Fourier transformation. There are local and global patterns of QPI that depend on the band structure of the underlying physical system. Concisely, local patterns of QPI depend on the geometry of the dispersion Capriotti et al. (2003), and are further modified (suppressed or enhanced) by the internal structure of quasiparticles Wang and Lee (2003); Hanaguri et al. (2009). Global patterns of QPI, on the other hand, depend on the topology of band structure, e.g., a global structure of spin-momentum locking Guo and Franz (2010). Consequently, we can use QPI to infer the geometry and topology of the band structure, rendering QPI and FT-STS as facilities bridging the $\bf{r}$-space observations with the $\bf{k}$-space description of the underlying system.

There have been intensive studies of QPI for quantum materials Avraham et al. (2018), including metals Crommie et al. (1993); Sprunger et al. (1997), superconductors Hoffman et al. (2002); Capriotti et al. (2003); Wang and Lee (2003); Pereg-Barnea and Franz (2003); Hanaguri et al. (2009); Lee et al. (2009a); Zhang et al. (2013), graphene Rutter et al. (2007); Pereg-Barnea and MacDonald (2008); Brihuega et al. (2008); Dombrowski et al. (2017); Jolie et al. (2018), surface states of topological insulators Zhou et al. (2009); Lee et al. (2009b); Beidenkopf et al. (2011); Kohsaka et al. (2015, 2017), Weyl semimetals Batabyal et al. (2016); Inoue et al. (2016); Mitchell and Fritz (2016); Zheng and Zahid Hasan (2018), and nonsymmorphic materials Topp et al. (2017); Queiroz and Stern (2018); Zhu et al. (2018). The connection between QPI and ${\bf k}$-space spectral information is established through the joint density of states (JDOS) and its generalization (GJDOS) that takes internal structures of quasiparticles into consideration. The JDOS approach has been proposed at an earlier stage to explain hot spots in patterns of QPI for d-wave superconductor Hoffman et al. (2002); Capriotti et al. (2003); Wang and Lee (2003); McElroy et al. (2003), and has been applied for gapless topological systems in terms of spin-selective scattering probability (SSP) recently Beidenkopf et al. (2011); Inoue et al. (2016). A more rigorous treatment, however, should refer to Fourier transform of LDOS (FT-LDOS). Ref. Derry et al. (2015); Kohsaka et al. (2017) point out that GJDOS may give some false features, and is not appliable in general. Nevertheless, as a convenient tool GJDOS allows us to intuitively analyze and understand QPI patterns from band structures directly. It is desirable for a clarification of the applicability of GJDOS. Moreover, a systematic treatment that can relate local and global patterns of QPI to the geometry and topology of the band structure is still awaited.

In this paper, we adopt a reduced response function (RRF) approach for analyzing QPI patterns under different scattering and probing channels. The RRF includes GJDOS as its imaginary component, and is a faithful encoding of the QPI information. We will prove that GJDOS shares the same singularities with FT-LDOS, with few exceptions that can be excluded at first hand. Moreover, we reveal the existence of higher-degree singularities, which can give rise to more significant features (hot spots) in QPI. With the justified GJDOS, we derive its analytical expression that clearly shows how its singularities and coherent factors jointly determine the QPI patterns. Based on this expression, we propose indicators of topological numbers in ideal topological systems. We then testify and demonstrate the robustness of such indicators of topological numbers against distracting local features in QPI arising from complicated geometry of the band structure, by exploring the QPI patterns in several representative topological materials, including topological insulators $\text{Bi}_{2}\text{Te}_{3}$ and BiTeI, as well as the graphene family.

The paper is organized as follows. We first propose the reduced response function and derive its singular behavior in Sec. (II). Then, we give an analytical expression for GJDOS and apply it for analyzing global patterns in Sec. (III). We then numerally study QPI for topological materials in Sec. (IV), including surface states of topological insulators and the graphene family.

II Reduced Response function for quasiparticle interference

II.1 The reduced response function

For simplicity without loss of generality, we restrict our discusses of physical systems described by two-band Hamiltonians $\mathcal{H}({\bf k})=E_{0}({\bf k})+\mathbf{d}({\bf k})\cdot\boldsymbol{\sigma}$. Here $\boldsymbol{\sigma}=(\sigma_{x},\sigma_{y},\sigma_{z})$ denotes a vector of Pauli matrix for a spin (or a pseudospin), which can represent systems with two internal degrees of freedom, including spin-half, sublattice, or particle-hole. Impurity is described as $V(r)=\sum_{\beta}V^{\beta}(r)\sigma_{\beta}$, where $\beta=0,x,y,z$ stands for scattering channels. In the presence of an impurity, interference between scattered-in and out quasiparticles leads to a perturbation to the local density of states. The local density is related to the Green function,

[TABLE]

where $\sigma_{\alpha}$ ($\alpha=0,x,y,z$) represents probe channels Guo and Franz (2010). The Green function in ${\bf k}$-space can be written as

[TABLE]

where $G_{0}({\bf k},\omega)=[\omega-\mathcal{H}({\bf k})+i0^{+}]^{-1}$ is the unperturbed Green function with positive infinitesimal $0^{+}$. For a $\sigma_{\beta}$ impurity and under Born approximation, the T-matrix takes the form $T_{{\bf k}^{\prime}{\bf k}}(\omega)=V^{\beta}_{{\bf k}^{\prime}{\bf k}}\sigma_{\beta}$, where $V^{\beta}_{{\bf k}^{\prime},{\bf k}}$ is a Fourier transform of the potential strength for scattering channel $\sigma_{\beta}$. The scattering potential may be independent of ${\bf k}^{\prime}$ and ${\bf k}$, e.g., for a impurity with a delta potential $\delta({\bf r})$. Or it can depend on ${\bf k}^{\prime}-{\bf k}$ Capriotti et al. (2003). We mention that in general it can depend on both ${{\bf k}+{\bf p}}$ and ${\bf k}$, e.g., for spin-orbit scattering the scattering potential is $V_{{{\bf k}+{\bf p}},{\bf k}}=V_{0}\{1+ic[({{\bf k}+{\bf p}})\times{\bf k}]\cdot\boldsymbol{\sigma}\}$ Lee et al. (2009b); Kohsaka et al. (2017).

The perturbed local density for probing channel $\sigma_{\alpha}$ and scattering channel $\sigma_{\beta}$ turns to be

[TABLE]

QPI patterns can be captured by a Fourier transformation of the perturbed local density,

[TABLE]

where the response function reads

[TABLE]

In the presence of centrosymmetry, $\Lambda_{\alpha\beta}({\bf p},\omega)=\Lambda_{\alpha\beta}(-{\bf p},\omega)$, Eq. (4) reduces to $\delta n_{\alpha\beta}({\bf p},\omega)=-\frac{1}{\pi}{\mathrm{Im}}\,[\Lambda_{\alpha\beta}({\bf p},\omega)]$. We consider the centrosymmetry in this paper since it is applicable for many QPI, while the case $\Lambda_{\alpha\beta}({\bf p},\omega)\neq\Lambda_{\alpha\beta}(-{\bf p},\omega)$ is left for further investigation.

To establish a connection between FT-LDOS and GJDOS, we first decompose the response function $\Lambda_{\alpha\beta}({\bf p},\omega)$. Note H({\bf k})|\psi_{{\bf k}}^{\pm}\big{>}=E_{{\bf k}}^{\pm}|\psi_{{\bf k}}^{\pm}\big{>}, with eigenstates $\psi_{{\bf k}}^{\pm}(\mathbf{r})\sim u_{{\bf k},\pm}\text{exp}(i{\bf k}\cdot\mathbf{r})$ and eigenvalues $E_{{\bf k}}^{\pm}=E_{0}\pm d_{\bf k}$. Here we write $\mathbf{d(\bf{k})}=d_{\bf k}\hat{\bf{d}}_{\bf k}$ with $d_{\bf k}=|\mathbf{d(\bf{k})}|$ . Then $\Lambda_{\alpha\beta}$ can be explicitly written as,

[TABLE]

where

[TABLE]

Here we have briefly written $\omega_{\bf k}=\omega-E_{0}({\bf k})$. $F^{ss^{\prime}}_{\alpha\beta}$ is the spin coherent factor and we set as a real function (see Sec. (II.2) for more details), and the complex component has been put in $P_{\alpha\beta}$, which is momentum independent if we assume $V^{\beta}_{{{\bf k}+{\bf p}},{\bf k}}=e^{i\phi_{\beta}}|V^{\beta}_{{{\bf k}+{\bf p}},{\bf k}}|$, e.g., the phase shift in the scattering is independent of ${\bf p}$. The dominator in $\Lambda_{\alpha\beta}({\bf p},\omega)$ contributes to singularities and the numerator accounts for further enhancements or depressions, and a combination of both can explain hot-spot features in QPI. Introducing

[TABLE]

we can write FT-LDOS (omitting a factor of $-\frac{1}{\pi}$),

[TABLE]

Since an autocorrelation of $AA$-type dominates typically a $BB$-type one in the second line of Eq. (II.1), we chose only the $AA$ term which corresponds to the joint density of states. As for an approximate yet faithful encoding of $\delta n_{\alpha\beta}({\bf p},\omega)$, we propose a so-called reduced response function $\mathcal{R}_{\alpha\beta}({\bf p},\omega)$

[TABLE]

where

[TABLE]

Note that $\mathcal{S}$ and $\mathcal{J}$ are related by the Hilbert transformation (or the Kramers-Kronig relation). The generalized joint density of states (GJDOS) $\mathcal{J}$ incorporates spin coherent factor into JDOS, and has been applied to analyze QPI for gapless topological systems in terms of spin-selective scattering probability Beidenkopf et al. (2011); Batabyal et al. (2016). The reduced response function contributes to FT-LDOS $\delta n_{\alpha\beta}({\bf p},\omega)$ by a projection onto the direction $P_{\alpha\beta}$ in the complex plane, namely

[TABLE]

II.2 Spin coherent factor

Now we derive the expression of $P_{\alpha\beta}F^{ss^{\prime}}_{\alpha\beta}$ and reveal the meaning of the spin coherent factor $F^{ss^{\prime}}_{\alpha\beta}$. We exploit a representation of density matrix $\rho_{{\bf k},s}\equiv u_{{\bf k},s}u_{{\bf k},s}^{+}=\frac{1}{2}(1+s\hat{{\bf d}}_{\bf k}\cdot\boldsymbol{\sigma})$. Since ${\rm Tr}[abc]={\rm Tr}[bca]$, we have

[TABLE]

It should be reminded that $j,l=x,y,z$ while $\alpha,\beta=0,x,y,z$. Some interesting results can be immediately derived for the summation term denoted as $t^{ss^{\prime}}_{\alpha\beta}({{\bf k}+{\bf p}},{\bf k})$. For $\alpha$ and $\beta$, there are three different cases:

$\alpha=\beta$. Then $t^{ss^{\prime}}_{\alpha\alpha}=2(1+ss^{\prime}\hat{{\bf d}}_{{\bf k}+{\bf p}}\cdot\hat{R}_{\alpha}\hat{{\bf d}}_{\bf k})$. Here $\hat{R}_{\alpha}$ is a mirror reflection with the $\alpha=x,y,z$ axis and $R_{0}$ means no operation. 2. 2.

$\alpha=0,\beta\neq 0$. Then $t^{ss^{\prime}}_{0\beta}=2iss^{\prime}L^{\beta}$, where $(L^{x},L^{y},L^{z})=\hat{{\bf d}}_{{\bf k}+{\bf p}}\times\hat{{\bf d}}_{\bf k}$. By symmetry we have $t^{ss^{\prime}}_{\beta 0}=t^{ss^{\prime}}_{0\beta}$. 3. 3.

$\alpha\neq 0,\beta\neq 0,\alpha$. Then $t^{ss^{\prime}}_{\alpha\beta}=2ss^{\prime}(\hat{d}^{\alpha}_{\bf k}\hat{d}^{\beta}_{{\bf k}+{\bf p}}+\hat{d}^{\beta}_{\bf k}\hat{d}^{\alpha}_{{\bf k}+{\bf p}})$.

Here the second case is remarkable. An factor of $i$ appears here. If $V^{\beta}_{{{\bf k}+{\bf p}},{\bf k}}$ is a real number, then we have $P_{0\beta}=i$, meaning that the GJDOS should account for quasiparticle interference and thus is expected to give sharp QPI features. However, those features may be depressed, for instance, in backscattering processes, where $t^{ss^{\prime}}_{0\beta}(-{\bf k},{\bf k})=0$. Nevertheless, GJDOS can be used to explain QPI from scattering between quasiparticles locating at cusps of CCE for time-reversal breaking systems, such as d-wave superconductors Capriotti et al. (2003); Wang and Lee (2003).

II.3 Singular behaviors of the reduced repsonse function

The reduced response function can be written as two parts. The first part corresponds to terms of $AB+iAA$ type and the second part corresponds to $BA+iAA$ type, as written in Eq. (9) and Eq. (II.1). In the following we give a derivation of the first part and the case of the second part can be obtained similarly. Our result shows that both the real and imaginary parts share the same singularities. To focus on the singular behavior, we omit the factor $F_{\alpha\beta}^{ss^{\prime}}$ temporarily. The first part of RRF is

[TABLE]

Using $1/(x+0^{+})=P(1/x)-i\pi\delta(x)$, the first part of RRF can be rewritten as

[TABLE]

The second part $\sum_{ss^{\prime}}r^{ss^{\prime}}_{BA}({\bf p},\omega)$ is obtained by a replacement ${{\bf k}+{\bf p}}\leftrightarrow{\bf k}$.

The integral is along CCEs of $E^{s}_{\mathbf{k}}=\omega$. We may parametrize the $n$-th CCE as $\mathbf{k}^{n,s}_{\omega}(t)=(f^{n,s}_{\omega}(t),g^{n,s}_{\omega}(t))$ with $t\in[0,2\pi)$. For brevity we would omit the index $n$. After a tedious derivation we can show that the integral is singular for ${\bf p}$ under a given energy $\omega$ only if there exists $t_{0}$ that

[TABLE]

The first condition promises that the scattered quasiparticle is on the CCE, and the second one points out that ${\bf p}$ should be specified. An illustration can be seen in Fig. (1) to show that backscattering meets the above conditions. Moreover, Eq. (II.3) servers as a mathematical foundation for exploring singular behavior analytically.

To study the singular behavior, we evaluate $r^{ss^{\prime}}_{AB}(\mathbf{p}+\bm{\delta},\omega)$ with $|\bm{\delta}|\sim 0$. An expression can be found as

[TABLE]

where we have defined the joint curvature,

[TABLE]

Here $\kappa^{s}_{\mathbf{k}}$ is the curvature of the CCE $\omega=E^{s}_{\bf k}$ at $\mathbf{k}$, $\mathbf{n}^{s}_{\mathbf{k}}=\mathbf{v}^{s}_{\mathbf{k}}/|\mathbf{v}^{s}_{\mathbf{k}}|$ the direction vector of the group velocity, $\mathbf{k}_{0}$ the real solution of Eq. (II.3). We remark a similar expression can be found for $r^{ss^{\prime}}_{AB}(\mathbf{p}+\bm{\delta},\omega)$ by an exchange of ${\bf k}_{0}+{\bf p}\leftrightarrow{\bf k}_{0}$.

The expression of singular behavior of Eq. (14) in terms of the joint curvature provides a clear geometrical tool for analysis: (1) The real and imaginary parts of $R(\mathbf{p},\omega)$ share the same singularities and diverge complementarily in a inverse-square-root manner, unless the joint curvature is zero. (2) If the joint curvature of certain singularity is zero, the exponent of the power-law divergence is $-2/3$ and can be even higher. In Fig. (2), we illustrate singularities with nonzero and zero joint curvatures, and they gives divergence of $-\frac{1}{2}$ and at least $-\frac{2}{3}$, separately.

II.4 Stability of singularities for finite life-time quasiparticles

We consider the stability/robustness of singularities of $\mathcal{R}$ when quasiparticles have finite life-time. We argue that singularities are not stable for $\mathcal{S}$ when a length of CCE is approximate flat. We simply take $0^{+}\rightarrow\eta$, where $\eta$ is a positive finite number accounting for the finite life-time for quasiparticles.

A finite life-time of quasiparticle will smooth singularities. Note that $B({\bf k},\omega)$ and $A({\bf k},\omega)$ exhibit distinct behavior around $\omega=E_{\bf k}$(see Fig. (3a)). Remarkably, $B({\bf k},\omega)$ is zero at $\omega=E_{\bf k}$ and changes sign at two sides. The real part of RRF, $\mathcal{S}$, is an autocorrelation of $A({\bf k},\omega)$ and $A({\bf k},\omega)$, has a large weighting at $|\omega-E_{\bf k}|<\eta$. For a large-momentum backscattering in Fig. (3b), scattered-out quasiparticle always has energy $\omega>E_{\bf k}$, and thus gives a positive contribution to $\mathcal{S}$ . However, when scattering occurs on an approximated flat CCE (compared to $|{\bf p}|$), (see Fig. (3c)), then scattered-out quasiparticle of both $\omega>E_{\bf k}$ and $\omega<E_{\bf k}$ contribute to $\mathcal{S}$, and those contributions cancel out. The case for near-nesting can be analyzed similarly. In contrast, the joint density of states $\mathcal{J}$ is an autocorrlation of two $A({\bf k},\omega)$, which receive always positive contributions around and mainly from $\omega=E_{\bf k}$. Thus, when QPI pattern should be given by $\mathcal{S}$, $\mathcal{J}$ would give false features for scattering on an approximately flat band. For instance, a hot spot at ${\bf p}=0$ always presents in $\mathcal{J}$ but is absent in $\mathcal{S}$ for $\eta\neq 0$. This can clarify a discrepancy between FT-LDOS and GJDOS for explaining QPI Derry et al. (2015).

III Global patterns in QPI for ideal topological systems

In the above section, we have revealed that $\mathcal{S}$ and $\mathcal{J}$ share same singularities ideally. While in a real application (where quasiparticles have finite life-time), singularities of $\mathcal{S}$ disappear for approximately flat CCE, but are remained in $\mathcal{J}$, leading to a discrepancy. However, those discrepancies correspond to rather specified QPI patterns, such as hot spot or an asterisk-like pattern around ${\bf p}=0$ Kohsaka et al. (2015). Thus, we can take them as exceptions at first hand when using $\mathcal{J}$. After having justified the applicable of $\mathcal{J}$, we will derive a simple expression of $\mathcal{J}$ that can be used directly for analyzing QPI pattern in a geometrical fashion. We then use this expression to analyze global patterns in QPI for ideal topological Dirac points with varied topological numbers.

III.1 Analytical expressions of GJDOS

Due to the properties of $\delta$-function， $J^{ss^{\prime}}_{\alpha\beta}({\bf p},\omega)$ can be analytically evaluated, :

[TABLE]

where the sum is taken over all $\mathbf{k}_{0}$’s in the real solution set of the following two equations,

[TABLE]

The singularity is obtained by setting the dominator

[TABLE]

which could be seen as a generalization of Van Hove singularity with regards to joint density of states. The condition $\eqref{jacobi}$ is reached when $\mathbf{v}^{s}_{\mathbf{k}_{0}}=0$ or $\mathbf{v}^{s^{\prime}}_{\mathbf{k}_{0}+\bf{p}}=0$, or $\mathbf{v}^{s}_{\mathbf{k}_{0}}\parallel\mathbf{v}^{s^{\prime}}_{\mathbf{k}_{0}+\bf{p}}$. The latter, for example, explains backscattering in metals where quasiparticle of ${\bf k}$ is scattered to $-{\bf k}$. The condition of singularity also could be viewed as an envelope curve for one-parameter $\mathbf{p}$ curve family $E^{s^{\prime}}_{\mathbf{k}+\bf{p}}=\omega$ and $E^{s}_{\mathbf{k}}=\omega$. This allows us to directly draw QPI pattern based on contours of constant energy. It should be pointed out that the condition Eq. (19) can be related to the stationary phase approximation Liu et al. (2012), which has been applied in the study of QPI for surface state of $\text{Bi}_{2}\text{Te}_{3}$.

The expression of Eq. (16), combined with geometrical representation of $F^{ss^{\prime}}_{\alpha\beta}$ as discussed in Sec. (II.2), allows an analysis of QPI pattern from the band structure directly. In the following, we will first give examples for simple models of gapless topological systems.

III.2 Global patterns in QPI

For a simple circle-like CCE, singularities come from backscattering from ${\bf k}_{0}$ to $-{\bf k}_{0}$. The topology of Dirac points manifests in their structures of spin-momentum binding, and will have effects on QPI patterns through the factor $F^{ss^{\prime}}_{\alpha\beta}$ that may suppress the singularities. For illustration, we choose the Hamiltonian,

[TABLE]

where $n$ is an arbitrary integer, as a model for a Dirac point with topological charge/number $Q=n$ Zhao and Wang (2013); Min and MacDonald (2008). Those systems have circle-like CCEs. Quasiparticle with eigenenergy $\pm|{\bf k}|^{n}$ has the wavefunction $u_{{\bf k},\pm}\text{exp}(i{\bf k}\cdot\mathbf{r})$, where $u_{{\bf k},\pm}=(1,\pm e^{-in\theta})^{T}$. We chose positive $\omega$, and indexes $s$ and $s^{\prime}$ in the RRF framework are omitted.

We first consider the charge scattering and charge probe channel ($\alpha=\beta=0$) for systems with different topological numbers $Q$. In Fig. (4) we demonstrate how to use GJDOS of Eq. (16) to analyze QPI patterns. A distinct behavior can be identified between $Q=1$ and $Q=2$ that origins from their distinct structures of spin-momentum locking: QPI patterns are suppressed for $Q=1$ while restored for $Q=2$ for backscattering.

We now consider $\alpha=\beta$ channel. Note that $F_{\alpha\alpha}=2(1+\hat{{\bf d}}_{{\bf k}+{\bf p}}\cdot\hat{R}_{\alpha}\hat{{\bf d}}_{\bf k})$. As the spin lies on $x-y$ plane, $R_{z}$ will not change $\hat{{\bf d}}_{\bf k}$. Interestingly, $F_{zz}$ is the same as that of $F_{00}$. Thus, although charge probe channel can not detect magnetic impurity $\sigma_{z}$ (the well-known prohibition of backscattering for Dirac fermion), a spin $\sigma_{z}$ resolved probe can Guo and Franz (2010). On the other hand, $R_{x/y}$ will change the sign of $\hat{d}^{y/x}_{\bf k}$. As a result, for channel $\alpha=\beta=x/y$ suppressions of singularities will become direct dependent. In Fig. (5) we show the case $\alpha=\beta=x$ for both $Q=1$ and $Q=2$, and QPI patterns have two and four separated hot arcs, separately. One can derive that the number of hot arcs turns to be $2n$ if the topological number $Q=n$, making it a global QPI pattern that can reveal topological number of the underlying system.

Following the above examples we can exploit global QPI patterns from different combinations of scattering and probe channels as indicators for topological numbers of the underlying systems. We give some topological-number indicators as following:

Odd-even indicator. For scattering channel $\alpha=0$ and probe channel $\beta=0$, the spin coherent factor reads as $F_{00}\sim(1+\cos n(\theta_{-{\bf k}_{0}}-\theta_{{\bf k}_{0}}))$. As $\theta_{-{\bf k}_{0}}-\theta_{{\bf k}_{0}}=\pi$, then $\mathcal{J}_{00}\sim\lim\limits_{\theta\rightarrow\pi}\frac{F_{00}}{\sin\theta}\sim 0$ is zero only for odd $n$. This can be understood from the picture of spin-momentum locking, spins would be antiparallel/parallel for Dirac point with odd/even topological number at a pair $(-{\bf k}_{0},{\bf k}_{0})$. Thus, circle-like QPI patterns remain/disappear for even/odd topological numbers. 2. 2.

Integer indicator. For scattering channel $\alpha=x$ and probe channel $\beta=x$, the spin coherent factor reads as $F_{xx}\sim(1+\cos n(2\theta_{{\bf k}_{0}}+\pi))$. There is a $2n$-fold periodicity with angle $\theta$. Consequently, there are $2|n|$ disconnected hot arcs in the QPI pattern. For $\alpha=\beta=y$, the conclusion holds while hot arcs rotate $\frac{\pi}{2n}$ compared to the case $\alpha=\beta=x$. 3. 3.

Positive-negative indicator. Counting hot arcs alone can not tell whether topological charge is positive or negative, thus we need an additional indicator. This can be achieved by rotating the scattering channel and probe channel a little, e.g., rotating clockwise, and observing how QPI pattern would rotate. If it rotates clockwise (anti-clockwise), then topological charge is positive (negative).

Further, we can explore quasiparticle interference of two different Dirac points. For instance, intervalley scattering in graphene where two valleys have opposite topological numbers Mallet et al. (2012). Consider two Dirac points with topological charge $Q_{1}=n$ and $Q_{2}=m$, following the analysis of $F_{\alpha\beta}$, we find that the odd-even indicator($\alpha=\beta=0$) would tell odd or even for $n-m$, and integers indicator would give $n+m$ disconnected hot arcs. Those results provide further clues to identify topological charges of different Dirac points, which would be demonstrated on QPI patterns of the graphene family.

IV Applications to topological matters

We now apply the reduced response function approach for topological materials Hasan and Kane (2010); Qi and Zhang (2011); Zhao and Wang (2013). We chose three representative examples, including surface states of topological insulators $\text{Bi}_{2}\text{Te}_{3}$ Fu (2009); Zhang et al. (2009) and BiTeI Crepaldi et al. (2012); Kohsaka et al. (2015), as well as the graphene family GEIM and NOVOSELOV . For those real materials there may be deviations of both dispersion and spin-momentum locking from ideal topological Dirac points described by ${\mathcal{H}}_{n}({\bf k})$ (Eq. (20)). Thus, they provides playgrounds for studying both local and global patterns of quasiparticle interference.

IV.1 Surface states of topological insulator: $\text{Bi}_{2}\text{Te}_{3}$

In the surface state of $\text{Bi}_{2}\text{Te}_{3}$, wrapping terms appears Fu (2009). At the low energy limit, the CCE can be approximated as a circle, while at larger energy it becomes non-convex. Such distinct geometries of CCEs will lead to different patterns in quasiparticles interference at different energies Zhou et al. (2009); Lee et al. (2009b). On one hand, we expect remarkable new features in QPI; on the other hand, we can testify whether indicators of topological number are robust under complex CCEs when the energy increases.

To study the effect risen up by the complicated geometry of CCE at large energy, we consider the surface state of topological insulator $\text{Bi}_{2}\text{Te}_{3}$ which can be modeled as Fu (2009); Lee et al. (2009b),

[TABLE]

where $\theta_{\bf k}=\arctan{k_{y}/k_{x}}$. The system ${\mathcal{H}}_{TI}({\bf k})$ has time-reversal symmetry and belongs to AII class. It owns a $Z_{2}$ type topological number of $Q=1$ Zhao and Wang (2013). Dispersions of two branches $s=\pm 1$ are opposite and we take only the positive one ($s=1$) for consideration in RRF. At a low energy, the CCE is circle-like, and QPI patterns under different combinations of scattering and probe channels are consistent with those of ideal Dirac point with $Q=1$ (see Fig. (6) for details).

For a sufficient large energy, the wrapping term becomes important and CCE turns to be nonconvex and exhibit six cusps. A scattering between a pair of cusps can satisfy the condition of singularity Eq. (19) but is not backscattering-type. Remarkably, ${\bf p}_{1}$ and ${\bf p}_{4}$ correspond to zero joint curvature. Moreover, ${\bf p}_{3}$ is a backscattering scattering between near-nesting arcs of CCE. Those vectors (${\bf p}_{1},{\bf p}_{3},{\bf p}_{4}$) thus have higher degrees of singularity than $\frac{1}{2}$, evidenced by hotter spots than ${\bf p}_{4}$ in joint density of states $J({\bf p},\omega)$ (See Fig. (7)). QPI under different combinations of scattering and probe channels exhibits more complex patterns as a consequence of both topological and geometrical aspects of the band structure. Numeral evaluation of $S_{00}({\bf p},\omega)$ shows hot spots locating at ${\bf p}_{1},{\bf p}_{4}$ that are absent for small energy when CCE is convex. Due to the spin coherent factor, spot is darker for ${\bf p}_{4}$ than that of ${\bf p}_{1}$. Those results fit very well with the FT-STS experiment Roushan et al. (2009); Beidenkopf et al. (2011). We also present $J_{00}({\bf p},\omega)$ as a comparison. It can be seen that $J_{00}({\bf p},\omega)$ has some sharp features around the center. These features come from joint densities of states with the same arc of the CCE, and are not stable singularities. Thus, a direct application of JDOS should take those as false features, when QPI should be given by the real part of the RRF.

We then investigate QPI in the other situations, e.g, $S_{xx}({\bf p},\omega)$ and $S_{zz}({\bf p},\omega)$. Compared with QPI at low energy, we can see that while local features are different, the global patterns persist. For $S_{zz}({\bf p},\omega)$, the bright closed curve results from backscattering. For $S_{xx}({\bf p},\omega)$, QPI pattern is cut into two parts in the middle, showing a 2-fold pattern. Remarkably, ${\bf p}_{4}$ is enhanced while ${\bf p}_{1}$ is suppressed, due to the effective $R_{x}$ reflection of the spin of scattered-in quasiparticle. These distortions of local features plus persistences of global patterns suggest the robustness of indicators for topological numbers. To reveal these QPI patterns, magnetic impurity with specified direction as well as spin-resolved STS techniques should be required, which are expected to be fulfilled in future FT-STS experiments Wiesendanger (2009); Jeon et al. (2017); Cornils et al. (2017).

IV.2 Surface state of polar semiconductor: BiTeI

We continue to study another topological state: surface state of polar semiconductor BiTeI Crepaldi et al. (2012); Kohsaka et al. (2015, 2017). There are two concentric CCEs, due to a contribution from the bulk state. Those two CCEs correspond to two branches of $s=1,-1$, and thus are opposite in spin-momentum locking along the CCE (see Fig. (8)). Scattering between two concentric CCEs can not be ignored. Moreover, the outer CCE has a larger distortion from the circle. Such a system thus exhibits unconventional properties of the band structure that can lead to novel patterns of quasiparticle interference.

Following Ref. Kohsaka et al. (2015), the surface state can be modeled as

[TABLE]

where $I$ is the identity matrix, and $k=\sqrt{k_{x}^{2}+k_{y}^{2}}$. The last term of Eq. (22) respects the C3v symmetry of BiTeI Ishizaka et al. (2011); Crepaldi et al. (2012). Functions are $E(k)=1+\alpha_{4}k^{2}+\alpha_{6}k^{4}$, $V(k)=v(1+\beta_{3}k^{2}+\beta_{5}k^{4})$ and $\Lambda(k)=\lambda(1+\gamma_{5}k^{2})$. Parameters as set as $m=0.0168$ eV*-1Å-2*, $\alpha_{4}=-2.03$ Å*-2*, $\alpha_{6}=87.5$ Å*-4*, $v=3.13$ eVÅ*-1*, $\beta_{3}=-2.01$ Å*-2*, $\beta_{5}=323$ Å*-4*, $\lambda=-41.7$ eVÅ*-3*, $\gamma_{5}=2.43$ Å*-2*, and $E_{0}=-0.352$ eV.

In the treatment of RRF, we should take all inter and intra CCEs scattering by a summation over $s=\pm 1$ and $s^{\prime}=\pm 1$ in Eq. (12). Numeral simulations are given in Fig. (8) for $\mathcal{S}_{00},\mathcal{S}_{zz},\mathcal{S}_{zx},\mathcal{S}_{xx}$. As we can see in $\mathcal{S}_{00}$, intra-CCE backscattering is prohibited for both the inner and the outer CCEs, due to a reversal of spin directions. However, inter-CCE scattering is allowed and a bright hexagon appears. QPI patterns in the case of $\mathcal{S}_{zz}$ is reversed: bright hexagon appears only for intra-CCE scattering. The most interesting case is spin $\sigma_{x}$ scattering and $\sigma_{x}$ probing channels: there are four disconnected bright arcs that respect a 2-fold symmetry. The outer arcs come from intra-CCE scattering: since $\hat{R}_{x}$ reflection will not change spin in $x$ direction, quasiparticle interference for ${\bf p}_{3}$ is still suppressed due to an approximately reversal of spins. However, hot arcs arises for wavevector ${\bf p}_{2}$, as $\hat{R}_{x}(\uparrow)=\downarrow$, thus the spin coherent factor now becomes nonzero. The inner arcs result from inter-CCE scattering: quasiparticle interference of wave vector ${\bf p}_{1}$ is enhanced while ${\bf p}_{3}$ is suppressed. The outer and inter hot arcs locate at horizon and vertical directions, respectively. Thus, they still respect the 2-fold symmetry. We also show $\mathcal{S}_{zx}$ that also respect the 2-fold symmetry. We remark that the $2$-fold symmetry can be considered as global patterns for indicting topological numbers $Q=1$, instead of merely counting the number of disconnected hot arcs.

IV.3 The graphene family

Graphene and its relatives of N-layer graphene host two topological Dirac points with opposite topological numbers Min and MacDonald (2008). Moreover, the topological number can vary with the layer. For example, there are two gapless points with $Q=\pm 2$ for AB stacked bilayer graphene and $Q=\pm 3$ for ABC stacked trilayer graphene. Those make the graphene family an ideal platform for studying global patterns of quasiparticle interferences.

We consider N-layer graphene with ABC stacking. For the i-th layer two sublattices are denoted as $a_{i}$ and $b_{i}$, then there are isolated atoms $a_{1}$ and $b_{N}$ that are not directly coupled to other layers. Taking $a_{1}$ and $b_{N}$ as a pseudospin, we have an effective Hamiltonian $H_{G,N}({\bf k})=({\mathrm{Re}}\,{f({\bf k})})^{N}\sigma_{x}+({\mathrm{Im}}\,{f({\bf k})})^{N}\sigma_{y}$ to model ABC stacking N-layer graphene Min and MacDonald (2008). Here $f({\bf k})=1+e^{i{\bf k}\cdot\bf{G}_{1}}+e^{i{\bf k}\cdot\bf{G}_{2}}$ with $\bf{G}_{1},\bf{G}_{2}$ being primitive vectors for graphene. Topological Dirac points locate at $\mathbf{K}$ and $\mathbf{K}^{\prime}$ with topological charges $Q=\pm N$ respectively. At low energy limit $H_{G,N}({\bf k})$ can be approximated by $H_{\pm N}({\bf k})$ with $f(k)\sim k^{N}$. It should be noted that the pseudospin here corresponds to sublattice, and the physical meaning of scattering and probing channel should be adjusted accordingly. For point-like impurity locating at site $a/b$ the scattering channel is $\tau_{a/b}=\frac{1}{2}(\sigma_{0}\pm\sigma_{z})$ (here we have shortnoted $a/b$ for $a_{1}/b_{N}$). Similarly, LDOS is measured at site $a/b$ and probing channel is $\tau_{a/b}$. The reduced response functions $\mathcal{R}_{00}({\bf p},\omega)$ and $\mathcal{R}_{zz}({\bf p},\omega)$ can be obtained as follows, $\mathcal{R}_{00}({\bf p},\omega)=\frac{1}{2}(\mathcal{R}_{aa}({\bf p},\omega)+\mathcal{R}_{bb}({\bf p},\omega)+\mathcal{R}_{ab}({\bf p},\omega)+\mathcal{R}_{ba}({\bf p},\omega)+)$,$\mathcal{R}_{zz}({\bf p},\omega)=\frac{1}{2}(\mathcal{R}_{aa}({\bf p},\omega)+\mathcal{R}_{bb}({\bf p},\omega)-\mathcal{R}_{ab}({\bf p},\omega)-\mathcal{R}_{ba}({\bf p},\omega))$, where $\mathcal{R}_{\mu\nu}$($\mu=a,b$ and $\nu=a,b$) correspond to probe channel $\tau_{\mu}$ and scattering channel $\tau_{\nu}$.

As two valleys $\bf{K}$ and $\bf{K^{\prime}}$ have the opposite topological charges, it is expected to see distinct QPI patterns between the center of ${\bf p}$-space due to intravalley scattering and those due to intervalley scattering away from the centre, such as around $\pm 2\mathbf{K}$, $\pm 2\mathbf{K}^{\prime}$ and $\pm 2(\mathbf{K}-\mathbf{K}^{\prime})$. We numeral evaluate $\mathbf{S}_{00}({\bf p},\omega)$ and $\mathbf{S}_{zz}({\bf p},\omega)$ for ABC stacked N-layer graphene for $N=1,2,3$ respectively, as seen in Fig. (9). For both channels, distinct $2N$ pieces of disconnected hot arcs appear for intravalley scattering and are rotated relatively. Those distinct features of QPI both for intervalley scattering indicate opposite topological charges for Dirac point located at the valleys. The simulations fit well with the FT-STS experiments Brihuega et al. (2008); Mallet et al. (2012) for monolayer graphene, while better resolution is required for bilayer graphene. For intravalley scattering, a bright circle appears for channels $\alpha=\beta=z$ and is absent $\alpha=\beta=0$ for monolayer graphene. The persistent presence of bright circle around the center for bilayer and trilayer graphene may result from high density of states from dispersions of $k^{N}$ ($N>1$). We mention that joint density of states give false hot spots at ${\bf p}=0,\pm 2\mathbf{K},\pm 2\mathbf{K}$, and $\pm 2(\mathbf{K}-\mathbf{K}^{\prime})$.

V Conclusions

We have developed a reduced response function approach for analyzing and simulating quasiparticle interference under different scattering and probing channels. The applicability of the generalized joint density of state has been clarified. We have analytically shown how singularities and spin coherent factor jointly determine QPI patterns. Remarkably, those local features and global patterns in QPI can be used to infer geometry details of dispersions and topology of band structures for the underlying system. We have also proposed indicators of topological numbers from global patterns of QPI for topological Dirac points. We have numerally simulated QPI of different topological materials, whose complicated geometrical features yet robust global patterns are evidenced and can be nicely captured in the RRF framework.

We remark that advances in spin-resolved STS technique Wiesendanger (2009); Jeon et al. (2017); Cornils et al. (2017) may experimentally visualize novel QPI patterns explored in this paper. Further investigations along this line include finding QPI that should be explained by the generalized joint density states, and an extension of the reduced response function approach for systems beyond two-band Hamiltonian description.

VI acknowledgement

The work was supported by the GRF (No. HKU173057/17P) and CRF (No. C6005-17G) of Hong Kong, and the NKRDP of China (Grant No. 2016YFA0301800).