Spin-orbit scattering visualized in quasiparticle interference

TL;DR

This study demonstrates that spin-orbit scattering significantly influences quasiparticle interference patterns, and introduces a method to detect it using spectroscopic-imaging STM, supported by numerical simulations on BiTeI.

Contribution

The paper presents a novel approach to visualize spin-orbit scattering in QPI patterns and extends simulation techniques to include spin-orbit effects.

Findings

Spin-orbit scattering affects QPI patterns in BiTeI.

Numerical simulations accurately reproduce experimental QPI data.

A selection rule makes spin-orbit scattering detectable in QPI patterns.

Abstract

In the presence of spin-orbit coupling, electron scattering off impurities depends on both spin and orbital angular momentum of electrons -- spin-orbit scattering. Although some transport properties are subject to spin-orbit scattering, experimental techniques directly accessible to this effect are limited. Here we show that a signature of spin-orbit scattering manifests itself in quasiparticle interference (QPI) imaged by spectroscopic-imaging scanning tunneling microscopy. The experimental data of a polar semiconductor BiTeI are well reproduced by numerical simulations with the $T$-matrix formalism that include not only scalar scattering normally adopted but also spin-orbit scattering stronger than scalar scattering. To accelerate the simulations, we extend the standard efficient method of QPI calculation for momentum-independent scattering to be applicable even for spin-orbit scattering. We further identify a selection rule that makes spin-orbit scattering visible in the QPI pattern. These results demonstrate that spin-orbit scattering can exert predominant influence on QPI patterns and thus suggest that QPI measurement is available to detect spin-orbit scattering.

Equations

$$H_0(k_x,k_y)$$

$$\rho (\mathbf{q},\omega )=-\frac{1}{2\pi i}\sum _{\mathbf{k}}\mathrm{Tr}\left\{\hat{G}(\mathbf{k},\mathbf{k}-\mathbf{q},\omega )-{\hat{G}}^{\ast }(\mathbf{k},\mathbf{k}+\mathbf{q},\omega )\right\},$$

$$\hat{G}(\mathbf{k},{\mathbf{k}}^{\prime },\omega )={\hat{G}}_0(\mathbf{k},\omega ){\delta }_{\mathbf{k},{\mathbf{k}}^{\prime }}+{\hat{G}}_0(\mathbf{k},\omega ){\hat{T}}_{\mathbf{k},{\mathbf{k}}^{\prime }}(\omega ){\hat{G}}_0({\mathbf{k}}^{\prime },\omega ),$$

$${\hat{T}}_{\mathbf{k},{\mathbf{k}}^{\prime }}(\omega )={\hat{V}}_{\mathbf{k},{\mathbf{k}}^{\prime }}+\sum _{\mathbf{p}}{\hat{V}}_{\mathbf{k},\mathbf{p}}{\hat{G}}_0(\mathbf{p},\omega ){\hat{T}}_{\mathbf{p},{\mathbf{k}}^{\prime }}(\omega ).$$

$${\hat{V}}_{\mathbf{k},{\mathbf{k}}^{\prime }}=V_0\{I+ic(\mathbf{k}\times {\mathbf{k}}^{\prime })\cdot \mathbf{\sigma }\},$$

$$\frac{dI}{dV}(\mathbf{r},V,V_{\mathrm{set}})\propto \frac{\rho (\mathbf{r},eV)}{{\int }_0^{e{V}_{\mathrm{set}}}\rho (\mathbf{r},E_{\mathrm{F}}+ϵ)dϵ},$$

$$\delta {\rho }_{\mathbf{k},{\mathbf{k}}^{\prime }}(\omega )$$

$$\delta {\rho }_{-{\mathbf{k}}^{\prime },-\mathbf{k}}(\omega )$$

$$\rho (\mathbf{q},\omega )\sim -\frac{1}{2\pi i}\sum _{\mathbf{k}}\mathrm{Tr}\left\{{\hat{G}}_0(\mathbf{k},\omega )\hat{T}(\omega ){\hat{G}}_0(\mathbf{k}-\mathbf{q},\omega )-{\hat{G}}_0^{\ast }(\mathbf{k},\omega ){\hat{T}}^{\ast }(\omega ){\hat{G}}_0^{\ast }(\mathbf{k}+\mathbf{q},\omega )\right\},$$

$$\sum _{\mathbf{k}}\mathrm{Tr}\left\{{\hat{G}}_0(\mathbf{k},\omega )\hat{T}(\omega ){\hat{G}}_0(\mathbf{k}+\mathbf{q},\omega )\right\}=\sum _{m,n}{t}^{(mn)}(\omega )\sum _{\mathbf{k},j}{g}_0^{(jm)}(\mathbf{k},\omega )g_0^{(nj)}(\mathbf{k}+\mathbf{q},\omega ),$$

$$\sum _{\mathbf{k}}{g}_0^{(jm)}(\mathbf{k},\omega )g_0^{(nj)}(\mathbf{k}+\mathbf{q},\omega )=\sum _{\mathbf{r}}{g}_0^{(jm)}(-\mathbf{r},\omega )g_0^{(nj)}(\mathbf{r},\omega )e^{i\mathbf{q}\cdot \mathbf{r}},$$

$${\hat{V}}_{\mathbf{k},{\mathbf{k}}^{\prime }}=\sum _j{\hat{u}}_j(\mathbf{k}){\hat{v}}_j({\mathbf{k}}^{\prime })={}^t\mathbf{u}(\mathbf{k})\mathbf{v}({\mathbf{k}}^{\prime }),$$

$${\hat{T}}_{\mathbf{k},{\mathbf{k}}^{\prime }}^{(n)}(\omega )=\{\begin{array}{ll}{\hat{V}}_{\mathbf{k},{\mathbf{k}}^{\prime }}& (n=1),\\ {\hat{V}}_{\mathbf{k},{\mathbf{k}}^{\prime }}+{\sum }_{\mathbf{p}}{\hat{V}}_{\mathbf{k},\mathbf{p}}{\hat{G}}_0(\mathbf{p},\omega ){\hat{T}}_{\mathbf{p},{\mathbf{k}}^{\prime }}^{(n-1)}(\omega )& (n\ge 2).\end{array}$$

$$K_n({\mathbf{k}}^{\prime },\omega )=\{\begin{array}{ll}M(\omega )\mathbf{v}({\mathbf{k}}^{\prime })& (n=1),\\ M(\omega )\mathbf{v}({\mathbf{k}}^{\prime })+M(\omega )K_{n-1}({\mathbf{k}}^{\prime },\omega )& (n\ge 2),\end{array}$$

$${\hat{T}}_{\mathbf{k},{\mathbf{k}}^{\prime }}^{(n)}(\omega )={}^t\mathbf{u}(\mathbf{k})\sum _{j=0}^{n-1}\{M(\omega ){\}}^j\mathbf{v}({\mathbf{k}}^{\prime }),$$

$$U=\left(\begin{array}{ccc}{V}_{0}I& 0& 0\\ 0& 0& ic{V}_0{\sigma }_z\\ 0& -ic{V}_0{\sigma }_z& 0\end{array}\right).$$

$$M(\omega )$$

$${\rho }^{\mathrm{JDOS}}(\mathbf{q},\omega )\propto \sum _{\mathbf{k}}{\rho }_0(\mathbf{k},\omega ){\rho }_0(\mathbf{k}-\mathbf{q},\omega ),$$

$${\rho }^{\mathrm{SJDOS}}(\mathbf{q},\omega )\propto \sum _{i=0,x,y,z}\sum _{\mathbf{k}}{\rho }_i(\mathbf{k},\omega ){\rho }_i(\mathbf{k}-\mathbf{q},\omega )$$

$$\displaystyle\rho^{\mathrm{SJDOS}}(\bm{q},\omega)\propto\sum_{\bm{k}}\sum_{m,n}\mathrm{Im}\big{\{}g_{n}(\bm{k},\omega)\big{\}}\mathrm{Im}\big{\{}g_{m}(\bm{k}-\bm{q},\omega)\big{\}}\Big{|}\big{\langle}\psi_{m}(\bm{k}-\bm{q})\big{|}\psi_{n}(\bm{k})\big{\rangle}\Big{|}^{2}.$$

$$\displaystyle\rho(\bm{q},\omega)=\sum_{\bm{k}}\delta\rho_{\bm{k},\bm{k}-\bm{q}}(\omega)\propto\sum_{\bm{k}}\sum_{m,n}\mathrm{Im}\big{\{}g_{n}(\bm{k},\omega)g_{m}(\bm{k}-\bm{q},\omega)\big{\}}\Big{|}\big{\langle}\psi_{m}(\bm{k}-\bm{q})\big{|}\psi_{n}(\bm{k})\big{\rangle}\Big{|}^{2}.$$

$$\delta {\rho }_{\mathbf{k},{\mathbf{k}}^{\prime }}(\omega )\equiv -\frac{1}{2\pi i}\mathrm{Tr}\left\{\hat{G}(\mathbf{k},{\mathbf{k}}^{\prime },\omega )-{\hat{G}}^{\ast }({\mathbf{k}}^{\prime },\mathbf{k},\omega )\right\}.$$

$$\delta {\rho }_{\mathbf{k},{\mathbf{k}}^{\prime }}(\omega )=-\frac{1}{2\pi i}\mathrm{Tr}\left[{\hat{G}}_0(\mathbf{k},\omega ){\hat{V}}_{\mathbf{k},{\mathbf{k}}^{\prime }}{\hat{G}}_0({\mathbf{k}}^{\prime },\omega )-{\left\{{\hat{G}}_0({\mathbf{k}}^{\prime },\omega ){\hat{V}}_{{\mathbf{k}}^{\prime },\mathbf{k}}{\hat{G}}_0(\mathbf{k},\omega )\right\}}^{\ast }\right].$$

$$\mathrm{Tr}\left\{{\hat{G}}_0(\mathbf{k},\omega ){\hat{V}}_{\mathbf{k},{\mathbf{k}}^{\prime }}{\hat{G}}_0({\mathbf{k}}^{\prime },\omega )\right\}$$

$$\displaystyle=\sum_{m,n}g_{n}(\bm{k},\omega)g_{m}(\bm{k}^{\prime},\omega)\mathrm{Tr}\left\{\big{|}\psi_{n}(\bm{k})\big{\rangle}\big{\langle}\psi_{n}(\bm{k})\big{|}\hat{V}_{\bm{k},\bm{k}^{\prime}}\big{|}\psi_{m}(\bm{k}^{\prime})\big{\rangle}\big{\langle}\psi_{m}(\bm{k}^{\prime})\big{|}\right\}$$

$$\displaystyle=\sum_{m,n}g_{n}(\bm{k},\omega)g_{m}(\bm{k}^{\prime},\omega)\big{\langle}\psi_{n}(\bm{k})\big{|}\hat{V}_{\bm{k},\bm{k}^{\prime}}\big{|}\psi_{m}(\bm{k}^{\prime})\big{\rangle}\big{\langle}\psi_{m}(\bm{k}^{\prime})\big{|}\psi_{n}(\bm{k})\big{\rangle}.$$

$$\displaystyle\mathrm{Tr}\left\{\hat{G}_{0}(\bm{k}^{\prime},\omega)\hat{V}_{\bm{k}^{\prime},\bm{k}}\hat{G}_{0}(\bm{k},\omega)\right\}^{*}=\sum_{m,n}\big{\{}g_{n}(\bm{k},\omega)g_{m}(\bm{k}^{\prime},\omega)\big{\}}^{*}\big{\langle}\psi_{n}(\bm{k})\big{|}\hat{V}_{\bm{k}^{\prime},\bm{k}}{}^{\dagger}\big{|}\psi_{m}(\bm{k}^{\prime})\big{\rangle}\big{\langle}\psi_{m}(\bm{k}^{\prime})\big{|}\psi_{n}(\bm{k})\big{\rangle}.$$

$$\displaystyle\delta\rho_{\bm{k},\bm{k}^{\prime}}(\omega)=-\frac{1}{\pi}\sum_{m,n}\mathrm{Im}\big{\{}g_{n}(\bm{k},\omega)g_{m}(\bm{k}^{\prime},\omega)\big{\}}\big{\langle}\psi_{n}(\bm{k})\big{|}\hat{V}_{\bm{k},\bm{k}^{\prime}}\big{|}\psi_{m}(\bm{k}^{\prime})\big{\rangle}\big{\langle}\psi_{m}(\bm{k}^{\prime})\big{|}\psi_{n}(\bm{k})\big{\rangle}.$$

$$\delta {\rho }_{-{\mathbf{k}}^{\prime },-\mathbf{k}}(\omega )$$

$$\displaystyle=-\frac{1}{\pi}\sum_{m,n}\mathrm{Im}\big{\{}g_{n}(-\bm{k},\omega)g_{m}(-\bm{k}^{\prime},\omega)\big{\}}\big{\langle}\psi_{m}(-\bm{k}^{\prime})\big{|}\hat{V}_{-\bm{k}^{\prime},-\bm{k}}\big{|}\psi_{n}(-\bm{k})\big{\rangle}\big{\langle}\psi_{n}(-\bm{k})\big{|}\psi_{m}(-\bm{k}^{\prime})\big{\rangle}$$

$$\displaystyle=-\frac{1}{\pi}\sum_{m,n}\mathrm{Im}\big{\{}g_{n}(\bm{k},\omega)g_{m}(\bm{k}^{\prime},\omega)\big{\}}\big{\langle}\psi_{n}(\bm{k})\big{|}\Theta\hat{V}_{-\bm{k}^{\prime},-\bm{k}}{}^{\dagger}\Theta^{-1}\big{|}\psi_{m}(\bm{k}^{\prime})\big{\rangle}\big{\langle}\psi_{m}(\bm{k}^{\prime})\big{|}\psi_{n}(\bm{k})\big{\rangle}.$$

$$\delta {\rho }_{-{\mathbf{k}}^{\prime },-\mathbf{k}}(\omega )$$

$$\displaystyle\bm{S}_{n}(\bm{k})\cdot\bm{S}_{m}(\bm{k}^{\prime})=2\Big{|}\big{\langle}\psi_{n}(\bm{k})\big{|}\psi_{m}(\bm{k}^{\prime})\big{\rangle}\Big{|}^{2}-1.$$

Full text

I Introduction

Spin-dependent scattering has played important roles in many fields of physics for a long time. Spin-dependent asymmetric scattering of electron beams in vacuum provided a foundation of relativistic quantum mechanics Mott (1929); Shull et al. (1943). In condensed matter physics, spin-orbit scattering of electrons propagating in solids contributes to some transport phenomena Hikami et al. (1980); Nagaosa et al. (2010); Sinova et al. (2015). For example, spin-dependent impurity scattering caused by spin-orbit scattering is among the origins of anomalous Hall effect and (extrinsic) spin Hall effect Nagaosa et al. (2010); Sinova et al. (2015). Another direct consequence of spin-orbit scattering is rotation of electron spin, changing interference between wave functions of electrons around an impurity. This effect on quantum interference is known as the origin of weak anti-localization Hikami et al. (1980).

Interference of wave functions results in a periodic modulation of the local density of states (LDOS). This modulation, known as quasiparticle interference (QPI), has been imaged by spectroscopic imaging scanning tunneling microscopy (SI-STM) in a wide variety of materials Hasegawa and Avouris (1993); Crommie et al. (1993); Hoffman et al. (2002); Zhang et al. (2009); Beidenkopf et al. (2011); Zeljkovic et al. (2014); Arguello et al. (2015); Zheng et al. (2016); Inoue et al. (2016). QPI has been studied mostly to acquire momentum-resolved information of electronic states from its characteristic periodicity based on an assumption that the scattering center is a scalar one. This assumption is widely used even for strong spin-orbit coupling systems Zhang et al. (2009); Beidenkopf et al. (2011); Zeljkovic et al. (2014); Zheng et al. (2016); Inoue et al. (2016) where spin-orbit scattering is also likely to be strong. However, the role of spin-orbit scattering in QPI is obscure.

A primal difference between spin-orbit and scalar scattering is that the former depends on both spin and momentum of electrons whereas the latter does not. Consequently, spin-orbit scattering can cause additional enhancement or suppression of QPI that is unanticipated for scalar scattering. Lee et al. theoretically indicate that spin-orbit scattering enhances new scattering channels for the surface states of topological insulator Bi2Te3 Lee et al. (2009). In the experiments, however, such enhancement and resultant multiple branches of QPI have never been observed Zhang et al. (2009); Beidenkopf et al. (2011).

In this paper, exploiting atomic-resolution SI-STM and numerical simulations, we reveal that spin-orbit scattering is predominant for QPI of the quasi-two-dimensional states of a polar semiconductor BiTeI. We performed a detailed analysis of QPI using the standard $T$-matrix formalism with an extended technique to accelerate calculations including not only momentum-independent scalar scattering but also momentum-dependent spin-orbit scattering. All the components of QPI observed by the experiments are successfully reproduced when spin-orbit scattering is considered. The key ingredient is the selection rule due to spin-orbit scattering, which selectively enhances one of the main scattering channels contributing QPI of BiTeI. Our finding demonstrates that spin-orbit scattering is actually observed in QPI and provides a foothold to get better insight into the electronic states through QPI.

II Overview of expermental results of BiTeI

BiTeI has a polar crystal structure with layered stacking of triple layers composed of Te, Bi, and I layers, and hosts giant Rashba-type spin splitting in bulk bands as well as at the surface as observed by angle-resolved photoemission spectroscopy (ARPES) Ishizaka et al. (2011); Crepaldi et al. (2012); Landolt et al. (2012); Sakano et al. (2012, 2013). A domain structure composed of opposite stacking orders is found in this material, resulting in two kinds of termination at a single surface Tournier-Colletta et al. (2014); Butler et al. (2014); Fiedler et al. (2014); Kohsaka et al. (2015). Especially in the Te-terminated areas, differential conductance ($dI/dV$) images show clear QPI of quasi-two-dimensional states split off from the bulk valence band by the spontaneous electric polarization Kohsaka et al. (2015).

To define points to be focused, we briefly summarize QPI of BiTeI. (Experimental details are described in Ref. Kohsaka et al., 2015.) Figure 1(a) shows a typical image of QPI in the Te-terminated area. The Fourier-transformed image shown in Fig. 1(b) reveals that the QPI patterns consist of three major components: the hexagonal ring surrounding the $\overline{\Gamma}$ point, the strong peaks in the $\overline{\Gamma}$-$\overline{\mathrm{M}}$ direction, and the outermost humps in the $\overline{\Gamma}$-$\overline{\mathrm{M}}$ direction. These features are commonly observed in all samples studied. By comparing QPI dispersions with the ARPES results, Ishizaka et al. (2011); Crepaldi et al. (2012); Landolt et al. (2012); Sakano et al. (2012) the hexagonal ring and the strong peaks are assigned to interband scattering between the spin-split bands and intraband scattering between the corners of the hexagonally warped outer branch, respectively, as depicted in Fig. 1(c). The outermost humps change their locations depending on bias voltages for stabilizing the scanning tip, and therefore are attributed to the so-called setpoint effect discussed later. The near-$\overline{\Gamma}$ feature, which varies from one sample to another, originates from nanometer-scale inhomogeneity due to random distribution of defects.

Although the positions of the hexagonal ring and the strong $\overline{\Gamma}$-$\overline{\mathrm{M}}$ peaks are understood as described above, there remains a puzzle in their intensities. The intraband scattering is mostly forbidden because the backscattering from $\bm{k}$ to $-\bm{k}$ is suppressed due to the anti-parallel spin orientations. The $\overline{\Gamma}$-$\overline{\mathrm{M}}$ peaks are nevertheless allowed because of deviation from the backscattering. Meanwhile, the interband scattering giving the hexagonal ring is always allowed because the spin orientations are almost parallel. That is, the spin texture of the band structure is more beneficial for the hexagonal ring than for the $\overline{\Gamma}$-$\overline{\mathrm{M}}$ peaks, although the former is actually weaker than the latter. This inverted intensity is a robust signature of spin-orbit scattering, as revealed below.

III The model and $T$-matrix formalism

To solve the puzzle of intensity, we numerically simulate QPI patterns. To model the quasi-two-dimensional state originating from the bulk valence band predominated by Bi 6$p_{z}$ orbitals, Bahramy et al. (2011) we employ an extended Rashba Hamiltonian,

[TABLE]

where $I$ and $\sigma_{i}$ $(i=x,y,z)$ are the identity matrix and the Pauli matrices, respectively, with $k=\sqrt{k_{x}^{2}+k_{y}^{2}}$. If $E(k)$ and $V(k)$ are constant and $\Lambda(k)=0$, Eq. (1) gives the Bychkov-Rashba Hamiltonian Bychkov and Rashba (1984). The last term of Eq. (1) reflects C3v symmetry of BiTeI Ishizaka et al. (2011); Bahramy et al. (2012). We extend $H_{0}$ up to $k^{6}$ terms, $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})$ so that it is invariant under a three-fold rotation along the $z$ direction, mirror operation about the $xz$ plane ($x$ is along the $\overline{\Gamma}$-$\overline{\mathrm{M}}$ direction), and the time-reversal operation. The higher terms up to $k^{6}$ are required for QPI calculations performed in the whole surface Brillouin zone whereas $k^{3}$ terms are enough to reproduce the ARPES results near the $\overline{\Gamma}$ point Sakano et al. (2013). We choose parameters 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 by fitting experimental data (Fig. 2).

QPI patterns have been calculated with the standard $T$-matrix formalism for a single local impurity. In fact, there are many defects in the field of view of Fig. 1(a). However, the three major features of QPI are independent of details of defect distribution as evidenced by the experimental fact that they are observed in all samples. Therefore, we postulate that multiple impurities work on overall intensity in a statistical manner Arguello et al. (2015) and a single impurity is a good starting point to discuss QPI patterns. The LDOS is written by the retarded Green’s function $\hat{G}$ in momentum space,

[TABLE]

where $\rho(\bm{q},\omega)$ is the Fourier transform of the LDOS, $\rho(\bm{q},\omega)=\int\rho(\bm{r},\omega)e^{-i\bm{q}\cdot\bm{r}}d\bm{r}$. Here we consider the Green’s function in matrix form to include spin. In the presence of a scattering center the potential of which in momentum space is $\hat{V}_{\bm{k},\bm{k}^{\prime}}$,

[TABLE]

where the $T$ matrix satisfies

[TABLE]

Here $\hat{G}_{0}$ is the bare Green’s function, \hat{G}_{0}(\bm{k},\omega)=\sum_{n}g_{n}(\bm{k},\omega)\big{|}\psi_{n}(\bm{k})\big{\rangle}\big{\langle}\psi_{n}(\bm{k})\big{|}, $g_{n}(\bm{k},\omega)=\left(\omega+i\eta-\epsilon_{n}(\bm{k})\right)^{-1}$, where $\epsilon_{n}(\bm{k})$ and $\psi_{n}(\bm{k})$ are the $n$th eigenvalue and the $n$th eigenstate of the bare Hamiltonian, respectively, with $\eta$ being a small broadening factor (10 meV for our simulations).

QPI patterns can be computed in principle with these equations. For momentum-independent scattering (e.g., scalar scattering), a direct calculation of the $\bm{k}$ summation requires $\mathcal{O}(N^{4})$ operations for a single QPI image, where $N\times N$ is the number of grid points. The amount of this calculation can be reduced to $\mathcal{O}(N^{2}\log_{2}N)$ by using fast Fourier transform (FFT) Arguello et al. (2015); Derry et al. (2015). For momentum-dependent scattering (e.g., spin-orbit scattering), however, the FFT-based technique has not been applied and consequently the $\bm{k}$ summation requiring $\mathcal{O}(N^{6})$ operations has been directly calculated. The enormous amount of calculation has hindered precise and comprehensive analysis of QPI; calculations have often been done only in a narrow range and at a low resolution of energy and momentum Lee et al. (2009). We find that the FFT-based technique is still available for momentum-dependent scattering satisfying a certain condition. Because of this method, the amount of calculation for spin-orbit scattering can be reduced to $\mathcal{O}(N^{2}\log_{2}N)$ that greatly accelerates our simulations. Details of the method are described in Appendix A.

IV Numerical results

We begin our simulations with a scalar impurity. In this case, the scattering potential is independent of momentum, $\hat{V}=V_{0}I$, where $V_{0}$ is strength of scattering ($V_{0}=0.1$ eV for all simulations). The $T$ matrix is also simplified to a momentum-independent form. This simplification makes the calculation greatly easy and is why a scalar impurity is widely assumed as scattering center. The simulation result successfully reproduces the hexagonal ring and the $\overline{\Gamma}$-$\overline{\mathrm{M}}$ peaks as shown in Fig. 3(a). The hexagonal ring appears at the location of interband scattering between the spin-split bands and the $\overline{\Gamma}$-$\overline{\mathrm{M}}$ peaks lie slightly outside of intraband scattering of the outer branch, corroborating the peak assignment described above. No prominent feature appears near the $\overline{\Gamma}$ point, being consistent with the experiments. (See Appendix B for details.) Although these basic features are reproduced, the ring is stronger than the peaks, replicating the puzzle. This discrepancy in intensities is robust as long as the band parameters are in a reasonable range, suggesting that the scattering is not a simple scalar one.

We then consider two kinds of scattering, magnetic scattering and spin-orbit scattering. Since both rotate electron spin, they may change the situation that is beneficial for the hexagonal ring. Figure 3(b) shows the calculation result with magnetic scattering $\hat{V}=V_{0}\sigma_{z}$ corresponding to a classical magnetic moment pointing in the $z$ direction. The $\overline{\Gamma}$-$\overline{\mathrm{M}}$ peaks are still weaker than the hexagonal ring and the QPI pattern is strongly suppressed overall. In contrast, noteworthy results are found for spin-orbit scattering,

[TABLE]

where $c$ is the effective spin-orbit coupling parameter Sinova et al. (2015) and denotes strength of spin-orbit scattering relative to that of scalar scattering. The $\overline{\Gamma}$-$\overline{\mathrm{M}}$ peaks are selectively enhanced as $c$ increases and become stronger than the hexagonal ring as shown in Fig. 3(c).

The remaining feature, the humps in the $\overline{\Gamma}$-$\overline{\mathrm{M}}$ direction, can be calculated by taking data acquisition procedures of SI-STM (the setpoint effect) into account. Even if a $dI/dV$ spectrum is proportional to LDOS at each location as generally assumed, the proportional constant is not generally uniform but has a spatial structure reflecting variation of the tip height. The height of a scanning tip is adjusted at each location such that the tunneling current is a set value. The current is determined by the LDOS integrated up to a given bias voltage $V_{\mathrm{set}}$. Consequently, a $dI/dV$ image observed by SI-STM depends on $V_{\mathrm{set}}$ as well as the LDOS,

[TABLE]

where $E_{\mathrm{F}}$ is the Fermi energy Kohsaka et al. (2007). The denominator of Eq. (6) represents the setpoint effect. This effect has been known in the experiments but neglected in the calculations of QPI. Full simulation including spin-orbit scattering and the setpoint effect is shown in Fig. 4. All of the hexagonal ring, the strong $\overline{\Gamma}$-$\overline{\mathrm{M}}$ peaks, and the $\overline{\Gamma}$-$\overline{\mathrm{M}}$ humps are well reproduced. The peak intensities agree with the experiment as shown in Fig. 4(e) with $c=80$ Å2 or a dimensionless parameter $c\pi^{2}/{a_{0}}^{2}=40$ ($a_{0}$ = 4.34 Å, $a$-axis length), indicating predominance of spin-orbit scattering over scalar scattering.

V Discussion

The above results of numerical simulations clearly show that spin-orbit scattering is the crucial ingredient to explain the QPI intensities of BiTeI. The contrasting results of magnetic and spin-orbit scattering can be understood as follows. Electrons with spin-up and spin-down feel potentials of opposite signs for magnetic scattering Zhou et al. (2009). In addition, electrons scattered to the right and the left do as well for spin-orbit scattering. The opposite signs result in suppression of QPI for magnetic scattering whereas the same sign, as a result of a combination of the two effects, is cooperative for QPI in the case of spin-orbit scattering. An essential point of this mechanism is the scattering amplitude. To the first order of $\hat{V}$, a contribution to QPI from a scattering process ($\bm{k}^{\prime}\rightarrow\bm{k}$) and its time-reversal counterpart ($-\bm{k}\rightarrow-\bm{k}^{\prime}$) are written as

[TABLE]

where $\Theta$ is the time-reversal operator. (Derivation of these formulas is written in Appendix C.) The difference between the two processes is found to be the potential in the scattering amplitude. For scalar scattering $\hat{V}_{\bm{k},\bm{k}^{\prime}}=V_{0}I$, Eqs. (7) and (8) are the same, $\delta\rho_{\bm{k},\bm{k}^{\prime}}=\delta\rho_{-\bm{k}^{\prime},-\bm{k}}$. For magnetic scattering $\hat{V}_{\bm{k},\bm{k}^{\prime}}=V_{0}\sigma_{i}$, the scattering amplitude changes sign under time reversal because $\Theta\sigma_{i}\Theta^{-1}=-\sigma_{i}$. Time-reversal processes thus always cancel with each other, $\delta\rho_{\bm{k},\bm{k}^{\prime}}=-\delta\rho_{-\bm{k}^{\prime},-\bm{k}}$, leading to the strong suppression of QPI for magnetic scattering. (This explains why QPI patterns are unchanged even in the presence of magnetic impurities Beidenkopf et al. (2011); Stróżecka et al. (2011).) For spin-orbit scattering $\hat{V}_{\bm{k},\bm{k}^{\prime}}=icV_{0}(\bm{k}\times\bm{k}^{\prime})\cdot\bm{\sigma}$, the scattering amplitude of time-reversal processes has the same sign because $\Theta\hat{V}_{-\bm{k},-\bm{k}^{\prime}}\Theta^{-1}=(-i)cV_{0}(\bm{k}\times\bm{k}^{\prime})\cdot(-\bm{\sigma})=\hat{V}_{\bm{k},\bm{k}^{\prime}}$. Therefore, as for spin-orbit scattering, the two scattering processes cooperatively contribute to QPI without being canceled as in the case of scalar scattering.

The preferential enhancement of the $\overline{\Gamma}$-$\overline{\mathrm{M}}$ peaks (Fig. 3(c)) is attributed to the directional and spin-dependent nature of spin-orbit scattering. For electrons in two dimensions, spin-orbit scattering is written as $\hat{V}_{\bm{k},\bm{k}^{\prime}}=icV_{0}kk^{\prime}\sin\theta_{\bm{k},\bm{k}^{\prime}}\sigma_{z}$, where $\theta_{\bm{k},\bm{k}^{\prime}}$ is the angle from $\bm{k}$ to $\bm{k}^{\prime}$. Since the hexagonal ring mainly consists of $\delta\rho_{\bm{k},\bm{k}^{\prime}}$ with $\theta_{\bm{k},\bm{k}^{\prime}}\sim\pi$ and $\langle\psi_{n}(\bm{k})|\sigma_{z}|\psi_{m}(\bm{k}^{\prime})\rangle\sim 0$, spin-orbit scattering does not contribute to the hexagonal ring whereas it does contribute to the $\overline{\Gamma}$-$\overline{\mathrm{M}}$ peaks because $\theta_{\bm{k},\bm{k}^{\prime}}\sim\pm 2\pi/3$ and $\langle\psi_{n}(\bm{k})|\sigma_{z}|\psi_{n}(\bm{k}^{\prime})\rangle\neq 0$. (Whether the scattering amplitude is zero or not is easily estimated as written in Appendix D.) This selective suppression is a selection rule originating from the scattering amplitude $\langle\psi_{n}(\bm{k})|\hat{V}_{\bm{k},\bm{k}^{\prime}}|\psi_{m}(\bm{k}^{\prime})\rangle=0$, which is distinct from a selection rule stemming from orthogonal wave functions $\langle\psi_{n}(\bm{k})|\psi_{m}(\bm{k}^{\prime})\rangle=0$ Petersen and Hedegård (2000); Pascual et al. (2004); Brihuega et al. (2008); Stróżecka et al. (2011); Beidenkopf et al. (2011).

Spin-orbit scattering exists in principle in any materials and grows with spin-orbit coupling. However, its appearance in QPI depends on details of relevant electronic states. QPI of surface states of Au(111) Hasegawa and Avouris (1993) is insensitive to spin-orbit scattering due to the selection rule, $\theta_{\bm{k},\bm{k}^{\prime}}\sim\pi$ and $\langle\psi_{n}(\bm{k})|\sigma_{z}|\psi_{m}(\bm{k}^{\prime})\rangle\sim 0$, being the same as the hexagonal ring of BiTeI. As for topological surface states of Bi2Te3, Zhang et al. (2009); Beidenkopf et al. (2011) we presume that spin-orbit scattering enhances the $\overline{\Gamma}$-$\overline{\mathrm{M}}$ peaks so it does for BiTeI, but its influence remains to be clarified. In this sense, the quasi-two-dimensional states of BiTeI with two scattering channels, one of which is sensitive to spin-orbit scattering and the other insensitive, are suited to investigate effects of spin-orbit scattering. QPI arising from electronic states near the Brillouin-zone boundary may be subject to spin-orbit scattering because of large momenta and $\theta_{\bm{k},\bm{k}^{\prime}}\sim\pi/2$ Lee et al. (2009). Such candidates are found in topological crystalline insulators and Weyl semimetals Zeljkovic et al. (2014); Zheng et al. (2016); Inoue et al. (2016).

Including other factors affecting QPI intensities would be conducive to better quantifying strength of spin-orbit scattering. A delta-function scattering potential is used for simplicity in our simulations. In a more realistic case, $V_{0}$ in Eq. (5) is changed from a constant to $V(|\bm{q}|)$ for a spherical scattering potential. Since $V(|\bm{q}|)$ usually decreases monotonically with increasing $|\bm{q}|$, QPI intensities are prone to be suppressed at large $|\bm{q}|$, where the amplitude of spin-orbit scattering is large. Finite sharpness of a scanning tip also causes a similar effect. Strength of spin-orbit scattering therefore may be underestimated due to these factors. Nevertheless, the strength of spin-orbit scattering obtained for BiTeI (80 Å2) is much larger than a theoretical value for $n$-GaAs (5.3 Å2), Sinova et al. (2015) being consistent with strong spin-orbit coupling in BiTeI. We note that the obtained value is averaged over many defects of multiple kinds. If defects are separated enough, spin-orbit scattering can be probed at individual defects and may be available for designing and optimizing materials of spin Hall effect.

VI conclusions

In conclusion, using atomic-resolution SI-STM and numerical simulations with the $T$-matrix formalism, we identify a signature of spin-orbit scattering in the QPI patterns of BiTeI. Spin-orbit scattering manifests itself in QPI through a selection rule originating from the scattering amplitude. Our results highlight the importance of the scattering process beyond featureless scalar scattering and, more importantly, suggest a potential capability of QPI measurement as a local, direct (unaffected by scattering time), and quantitative probe of spin-orbit scattering detected heretofore by transport measurements. We believe that including spin-orbit scattering into QPI analysis, which is now readily possible as demonstrated in our simulations, leads to a deeper understanding of electronic states of and more functionality from strong spin-orbit coupling systems.

Acknowledgements.

We thank M. S. Bahramy and Wei-Cheng Lee for fruitful discussions. This work is supported by the Murata Science Foundation.

Appendix A Fast calculation of QPI patterns

The $\bm{k}$ summation in Eq. (2) with the integral equation about the $T$-matrix Eq. (4) requires $\mathcal{O}(N^{4})$ operations for single $\bm{q}$, resulting in $\mathcal{O}(N^{6})$ operations in total to calculate a QPI image. In the case of momentum-independent scattering, the total amount of calculation is reduced to $\mathcal{O}(N^{4})$ and further reduced to $\mathcal{O}(N^{2}\log_{2}N)$ by using FFT. Here we show that, even in the case of momentum-dependent scattering, the amount of calculation can be reduced to $\mathcal{O}(N^{2}\log_{2}N)$ from $\mathcal{O}(N^{6})$ of the $\bm{k}$ summation when the scattering satisfies a certain condition.

To introduce our approach to reduce the amount of calculation, we begin with momentum-independent scattering. Since the $T$ matrix is also independent of $\bm{k}$, Eqs.(2)-(4) are simplified to

[TABLE]

where $\hat{T}(\omega)=\left(I-\hat{V}\sum_{\bm{k}}\hat{G}_{0}(\bm{k},\omega)\right)^{-1}\hat{V}$. Only the inhomogeneous part of LDOS is shown in Eq. (9) for brevity because it gives spatial modulations of QPI patterns. Each matrix element of the $\bm{k}$ summation in Eq. (9) is

[TABLE]

where $g_{0}^{(ij)}$ and $t^{(ij)}$ are matrix elements of $\hat{G}_{0}$ and $\hat{T}$, respectively. The right-hand side of Eq. (10) is the cross-correlation between $g_{0}^{(jm)}$ and $g_{0}^{(nj)}$, and thus can be expressed with Fourier transform,

[TABLE]

where $g_{0}^{(ij)}(\bm{r},\omega)=\sum_{\bm{k}}g_{0}^{(ij)}(\bm{k},\omega)e^{-i\bm{k}\cdot\bm{r}}$. Since the right-hand side of Eq. (11) is the inverse Fourier transform of $g_{0}^{(jm)}(-\bm{r},\omega)g_{0}^{(nj)}(\bm{r},\omega)$, the left-hand side of Eq. (11) can be calculated via FFT without taking the $\bm{k}$ summation Arguello et al. (2015); Derry et al. (2015). The amount of calculation is thus reduced from $\mathcal{O}(N^{4})$ of the direct calculation to $\mathcal{O}(N^{2}\log_{2}N)$ of FFT.

FFT, the essential point to reduce the amount of calculation, is available as in Eq. (11) because Eq. (9) is expressed virtually as a product of two matrices; one is a function of $\bm{k}$ and the other is a function of $\bm{k}+\bm{q}$. At a glance, this condition is not satisfied for general scatterers because the $T$ matrix depends on momentum. However, when $\hat{V}_{\bm{k},\bm{k}^{\prime}}$ is expressed as a sum of products between $\bm{k}$- and $\bm{k}^{\prime}$-dependent matrices

[TABLE]

where $\bm{u}(\bm{k})={}^{t}\!\left(\hat{u}_{1}(\bm{k})\ \hat{u}_{2}(\bm{k})\ \cdots\right)$ and $\bm{v}(\bm{k})={}^{t}\!\left(\hat{v}_{1}(\bm{k})\ \hat{v}_{2}(\bm{k})\ \cdots\right)$, FFT is available to calculate QPI patterns as in the case of momentum-independent scattering.

We rewrite Eq. (4) in a form of a recurrence formula

[TABLE]

By multiplying $\bm{v}(\bm{k})\hat{G}_{0}(\bm{k},\omega)$ from the left and taking a sum with respect to $\bm{k}$, we obtain

[TABLE]

where $K_{n}(\bm{k}^{\prime},\omega)=\sum_{\bm{p}}\bm{v}(\bm{p})\hat{G}_{0}(\bm{p},\omega)\hat{T}^{(n)}_{\bm{p},\bm{k}^{\prime}}(\omega)$ and $M(\omega)=\sum_{\bm{p}}\bm{v}(\bm{p})\hat{G}_{0}(\bm{p},\omega){}^{t}\!\bm{u}(\bm{p})$. Since these equations are summarized as $K_{n}(\bm{k}^{\prime},\omega)=\sum_{j=1}^{n}\{M(\omega)\}^{j}\bm{v}(\bm{k}^{\prime})$, Eq. (13) is also summarized as

[TABLE]

where $\{M(\omega)\}^{0}$ is the identity matrix, ${}^{t}\!\bm{u}(\bm{k})\{M(\omega)\}^{0}\bm{v}(\bm{k}^{\prime})={}^{t}\!\bm{u}(\bm{k})\bm{v}(\bm{k}^{\prime})$. Now the $T$ matrix and thus the second term of Eq. (3) are expressed as a product of $\bm{k}$ and $\bm{k}^{\prime}$ terms, and FFT is available to calculate QPI patterns as described above. We stress that the integral equation of the $T$ matrix [Eq. (4)] is reduced virtually to a sum of $\{M(\omega)\}^{j}$, which is independent of $\bm{k}$. Therefore, the calculation size is almost the same as that of momentum-independent scatterers. Namely, the total amount of QPI calculation is drastically reduced from $\mathcal{O}(N^{6})$ to $\mathcal{O}(N^{2}\log_{2}N)$.

Spin-orbit scattering, $\hat{V}_{\bm{k},\bm{k}^{\prime}}=V_{0}\left\{I+ic(\bm{k}\times\bm{k}^{\prime})\cdot\bm{\sigma}\right\}$, is written in the form of Eq. (12). In two dimension, it is written as $\hat{V}_{\bm{k},\bm{k}^{\prime}}=V_{0}\left\{I+ic(k_{x}k_{y}^{\prime}-k_{y}k_{x}^{\prime})\sigma_{z}\right\}={}^{t}\!\bm{u}(\bm{k})\bm{v}(\bm{k}^{\prime})$, where ${}^{t}\!\bm{u}(\bm{k})=(I\ k_{x}I\ k_{y}I)$, $\bm{v}(\bm{k})=U\bm{u}(\bm{k})$, and

[TABLE]

The $M$ matrix is

[TABLE]

Note that each element of $M$ is a $2\times 2$ matrix independent of $\bm{k}$ and thus $\sum M^{j}$ is readily calculable. Calculation of Fig. 3(c) with this method took only several seconds with a desktop computer. This means direct $\bm{k}$ summation at the same resolution is estimated to take about several years. Even if the direct calculations were done at lower resolutions, it may be difficult to find a reasonable compromise between resolution and time as shown in Fig. 5, highlighting drastic reduction of calculation costs with our method.

Appendix B Limitations of the joint density of states approach

All QPI patterns observed in the experiments show large intensities near $\bm{q}=0$ as shown in Fig. 1(b). Since these near-$\overline{\Gamma}$ features vary from one field of view to another, and one sample to another, they originate from the nanoscale inhomogeneity due to random distribution of defects. The varying near-$\overline{\Gamma}$ feature means that a QPI pattern near $\bm{q}=0$, if any, is small and masked by the nanoscale inhomogeneity. Therefore, no prominent feature near $\bm{q}=0$ in Fig. 3 is consistent with the experiments.

Meanwhile, from the viewpoint of the so-called joint density of states (JDOS), one may expect a large QPI intensity near $\bm{q}=0$; if DOS is large at a given $\bm{k}$, large DOS is also found near $\bm{k}$, resulting in a large QPI intensity near $\bm{q}=0$. However, as revealed by JDOS calculations shown below, such a JDOS-derived pattern near $\bm{q}=0$ has never been observed in BiTeI.

In the JDOS approach, QPI patterns are interpreted to be proportional to JDOS,

[TABLE]

where $\rho_{0}(\bm{k},\omega)$ is the DOS at $\bm{k}$, \rho_{0}(\bm{k},\omega)=-\mathrm{Im}\Big{[}\mathrm{Tr}\big{\{}G_{0}(\bm{k},\omega)\big{\}}\Big{]}. To include a spin effect, spin-dependent JDOS is considered,

[TABLE]

where \rho_{i}(\bm{k},\omega)=-\mathrm{Im}\Big{[}\mathrm{Tr}\big{\{}\sigma_{i}G_{0}(\bm{k},\omega)\big{\}}\Big{]}. The JDOS approach always predicts a large QPI intensity near $\bm{q}=0$ because the JDOS is the auto-correlation of $\rho_{i}(\bm{k},\omega)$. Actually, as shown in Fig. 6, calculations with the JDOS and spin-dependent JDOS approaches show an asterisk-like pattern centered at $\bm{q}=0$, sticking out in the $\overline{\Gamma}$-$\overline{\mathrm{K}}$ direction, and extending close to the hexagonal ring of intraband scattering. However, such a salient pattern has never been observed in the experiments. This discrepancy between the JDOS calculations and the observed QPI patterns highlights limitations of the JDOS approach.

The limitations of the JDOS approach derive from a difference between QPI and the JDOS. The spin-dependent JDOS of Eq. (19) can be written as

[TABLE]

A QPI pattern calculated by the $T$-matrix formalism for scalar scattering is written in a similar form. Equation (7) with $\hat{V}_{\bm{k},\bm{k}^{\prime}}=V_{0}I$ gives

[TABLE]

The product of $\mathrm{Im}\,g$ in Eq. (20) is not included in Eq. (21); in other words, that QPI is not directly related to JDOS as discussed in Ref. Derry et al., 2015.

Appendix C Contribution to QPI from $\bm{k}\rightarrow\bm{k}^{\prime}$ scattering process

We define $\delta\rho_{\bm{k},\bm{k}^{\prime}}(\omega)$ such that its summation with respect to $\bm{k}$ gives $\rho(\bm{k}-\bm{k}^{\prime},\omega)$,

[TABLE]

Here we consider an approximation to the first order of $\hat{V}_{\bm{k},\bm{k}^{\prime}}$ for simplicity. Since the $T$ matrix is $\hat{T}_{\bm{k},\bm{k}^{\prime}}\sim\hat{V}_{\bm{k},\bm{k}^{\prime}}$,

[TABLE]

The first term is

[TABLE]

Similarly, the second term is

[TABLE]

Given $\hat{V}_{\bm{k}^{\prime},\bm{k}}{}^{\dagger}=\hat{V}_{\bm{k},\bm{k}^{\prime}}$ as $\hat{V}$ is Hermitian, Eqs. (25) and (26) are summarized to

[TABLE]

Replacing $\bm{k}$ ($\bm{k}^{\prime}$) with $-\bm{k}^{\prime}$ ($-\bm{k}$) gives the time-reversal counterpart,

[TABLE]

The second follows from three identities; $g_{n}(-\bm{k},\omega)=g_{n}(\bm{k},\omega)$ as $H_{0}$ is time invariant, and \big{\langle}\beta\big{|}{\cal L}\big{|}\alpha\big{\rangle}=\big{\langle}\alpha^{\prime}\big{|}\Theta{\cal L}^{\dagger}\Theta^{-1}\big{|}\beta^{\prime}\big{\rangle} and \big{\langle}\beta\big{|}\alpha\big{\rangle}=\big{\langle}\alpha^{\prime}\big{|}\beta^{\prime}\big{\rangle}, where \big{|}\alpha^{\prime}\big{\rangle}=\Theta\big{|}\alpha\big{\rangle}, \big{|}\beta^{\prime}\big{\rangle}=\Theta\big{|}\beta\big{\rangle}, and $\cal L$ is a linear operator Sakurai (1993). The last two hold because $\Theta$ is antiunitary. Given $\hat{V}_{-\bm{k}^{\prime},-\bm{k}}{}^{\dagger}=\hat{V}_{-\bm{k},-\bm{k}^{\prime}}$ as $\hat{V}$ is Hermitian, we obtain

[TABLE]

Appendix D Relation between the scattering amplitude and spin orientation

Let $\bm{S}_{n}(\bm{k})$ be the expectation value of the Pauli matrices, \bm{S}_{n}(\bm{k})=\big{\langle}\psi_{n}(\bm{k})\big{|}\bm{\sigma}\big{|}\psi_{n}(\bm{k})\big{\rangle}, then the following relation holds

[TABLE]

The $\pi$-rotation operator along the $z$ direction for a spin-1/2 state is given as $e^{-i\pi\sigma_{z}/2}=-i\sigma_{z}$. Therefore, by rotating $\bm{S}_{m}(\bm{k}^{\prime})$ and \big{|}\psi_{m}(\bm{k}^{\prime})\big{\rangle} by $\pi$ along the $z$ direction, we obtain

[TABLE]

where $R_{z\pi}$ denotes $\pi$-rotation of spin orientation along the $z$ direction. Equation (30) means \big{\langle}\psi_{n}(\bm{k})\big{|}\sigma_{z}\big{|}\psi_{m}(\bm{k}^{\prime})\big{\rangle}=0 for two states with the spin orientations parallel and lying in the $xy$ plane.