Visualizing Pairing Symmetry in Nematic Superconductors using spin-polarized spectroscopy of magnetic impurities

TL;DR

This paper investigates how spin-polarized spectroscopy of magnetic impurities can reveal pairing symmetry in nematic superconductors, providing a method to distinguish different pairing states through spectral analysis.

Contribution

It introduces a detailed analysis of spin-polarized local density of states in nematic superconductors, linking spectral features to pairing symmetry and impurity orientation.

Findings

SP LDOS can distinguish pairing symmetries in nematic superconductors.

Fourier-transformed SP LDOS reveals quasiparticle interference patterns.

Method applies to both 2D and 3D materials with different impurity configurations.

Abstract

We study the spin-polarized spectral properties of Yu-Shiba-Rusinov resonance states induced by magnetic impurities in 2- and 3-dimensional nematic superconductors: few layer Bi$_2$Te$_3$ grown on FeTe$_{0.55}$Se$_{0.45}$ (2-dimensional) and Cu$_x$Bi$_2$Se$_3$ (3-dimensional). We focus on the relationship between pairing symmetry and the topograph of spin-polarized spectroscopy. We calculate the spin-polarized local density of states (SP LDOS) and the corresponding Fourier transformation using the $T$-matrix method for both the 2- and 3-dimensional materials. Various situations with different impurity orientations and different SP LDOSs have been investigated. We find that, like the quasiparticle interference spectrum, the SP LDOS can be applied to distinguish other pairings which preserve the threefold rotation symmetry of $D_3$ point group and nematic pairings in these materials.

Tables

Table 1. Table 1: Parameters used for numerical calculation of the tight-binding model Hao  et al. ( 2014 ); Hao and Ting ( 2017 ) , in units of electron volts (eV).

	$C_0$	$C_1$	$C_2$	$M_0$	$M_1$
CuxBi2Se3	$-$0.008	0.06	1.0	$-$0.26	0.3
Bi2Te3	$-$0.18	0.0634	2.59	$-$0.3	0.102

Table 2. Table 2: Matrix representation of the elementary operations for D 3 ​ d subscript 𝐷 3 𝑑 D_{3d} point group. 𝒞 z subscript 𝒞 𝑧 \mathscr{C}_{z} refers to the threefold rotation operation along the z 𝑧 z -axis. 𝒞 y ′ subscript superscript 𝒞 ′ 𝑦 \mathscr{C}^{\prime}_{y} refers to the twofold rotation operation along the y 𝑦 y -axis. ℐ ℐ \mathcal{I} is the inversion operation. Here ζ 0 subscript 𝜁 0 \zeta_{0} and ζ z subscript 𝜁 𝑧 \zeta_{z} are the 2 × 2 2 2 2\times 2 identity matrix and the third Pauli matrix in the Nambu spinor space.

operations	${𝒞}_z$	${𝒞}_y^{\prime }$	$ℐ$
matrix rep. (normal state)	${\sigma }_0\otimes e^{2\pi i{s}_z/3}$	$i{\sigma }_0\otimes s_x$	${\sigma }_x\otimes s_0$
matrix rep. (BdG)	${\sigma }_0\otimes e^{2\pi i{s}_z\otimes {\zeta }_z/3}$	$i{\sigma }_0\otimes s_x\otimes {\zeta }_z$	${\sigma }_x\otimes s_0\otimes {\zeta }_0$

Equations

$$\hat{H}=\sum _{\mathbf{k}}{\psi }_{\mathbf{k}}^{†}h(\mathbf{k}){\psi }_{\mathbf{k}},$$

$$h_{3d}(\mathbf{k})=\epsilon (\mathbf{k}){\Gamma }_0+M(\mathbf{k}){\Gamma }_5+B_0{c}_z{\Gamma }_4+A_0[c_y(\mathbf{k}){\Gamma }_1$$

$$-c_x(\mathbf{k}){\Gamma }_2]+R_1{d}_1(\mathbf{k}){\Gamma }_3+R_2{d}_2(\mathbf{k}){\Gamma }_4.$$

$$h_{2d}(\mathbf{k})=\epsilon (\mathbf{k}){\Gamma }_0+M(\mathbf{k}){\Gamma }_5+A_0[c_y(\mathbf{k}){\Gamma }_1-c_x(\mathbf{k}){\Gamma }_2]$$

$$+R_1{d}_1(\mathbf{k}){\Gamma }_3+R_2{d}_2(\mathbf{k}){\Gamma }_4,$$

$${\hat{\mathcal{H}}}_{\text{BdG}}=\frac{1}{2}\sum _{\mathbf{k}}{\Psi }_{\mathbf{k}}^{†}\left(\begin{array}{cc}h(\mathbf{k})& \Delta \\ {\Delta }^{†}& -h^{\mathsf{T}}(-\mathbf{k})\end{array}\right){\Psi }_{\mathbf{k}},$$

$${\mathcal{G}}_0(E,\mathbf{k})={\left[E+i{0}^{+}-\left(\begin{array}{cc}h(\mathbf{k})& \Delta \\ {\Delta }^{†}& -h^{\mathsf{T}}(-\mathbf{k})\end{array}\right)\right]}^{-1}.$$

$${\mathcal{V}}^{\mathrm{NM}}$$

$${\mathcal{V}}^{\mathrm{M}}$$

$$\mathcal{T}(E)=[{\mathcal{V}}^{-1}-{\mathcal{G}}_0(E,0){]}^{-1},$$

$${\mathcal{G}}_0(E,\mathbf{R})=\frac{1}{N}\sum _{\mathbf{k}}{e}^{i\mathbf{k}\cdot \mathbf{R}}{\mathcal{G}}_0(E,\mathbf{k}),$$

$$\mathcal{G}(E,\mathbf{R})={\mathcal{G}}_0(E,\mathbf{R})+\Delta \mathcal{G}(E,\mathbf{R}),$$

$$\Delta \mathcal{G}(E,\mathbf{R})={\mathcal{G}}_0(E,\mathbf{R})\mathcal{T}(E){\mathcal{G}}_0(E,-\mathbf{R}).$$

$$S_z(E,\mathbf{R})=-\frac{1}{\pi }\text{Im}[\Delta {\mathcal{G}}_{11}-\Delta {\mathcal{G}}_{22}+\Delta {\mathcal{G}}_{33}-\Delta {\mathcal{G}}_{44}],$$

$$S_x(E,\mathbf{R})=-\frac{1}{\pi }\text{Im}[\Delta {\mathcal{G}}_{12}+\Delta {\mathcal{G}}_{21}+\Delta {\mathcal{G}}_{34}+\Delta {\mathcal{G}}_{43}],$$

$$S_y(E,\mathbf{R})=-\frac{1}{\pi }\text{Re}[\Delta {\mathcal{G}}_{12}-\Delta {\mathcal{G}}_{21}+\Delta {\mathcal{G}}_{34}-\Delta {\mathcal{G}}_{43}],$$

$$\delta \rho (E,\mathbf{q})=\sum _{\mathbf{R}}{e}^{i\mathbf{q}\cdot \mathbf{R}}\delta \rho (E,\mathbf{R}),$$

$$\delta \rho (E,\mathbf{R})=-\frac{1}{\pi }\text{Im}[\Delta {\mathcal{G}}_{11}+\Delta {\mathcal{G}}_{22}+\Delta {\mathcal{G}}_{33}+\Delta {\mathcal{G}}_{44}].$$

$${\mathcal{C}}_y^{\prime -1}\left(\begin{array}{cc}0& \Delta \\ {\Delta }^{†}& 0\end{array}\right){\mathcal{C}}_y^{\prime }=\left(\begin{array}{cc}0& \Delta \\ {\Delta }^{†}& 0\end{array}\right).$$

$${\mathcal{C}}_y^{\prime -1}{h}_{\mathrm{BdG}}(k_x,k_y){\mathcal{C}}_y^{\prime }=h_{\mathrm{BdG}}(-k_x,k_y),$$

$$h_{\mathrm{BdG}}(k_x,k_y)=\left(\begin{array}{cc}h(k_x,k_y)& \Delta \\ {\Delta }^{†}& -h^{\mathsf{T}}(-k_x,-k_y)\end{array}\right).$$

$${\mathcal{C}}_y^{\prime -1}{\mathcal{G}}_0(E,k_x,k_y){\mathcal{C}}_y^{\prime }={\mathcal{G}}_0(E,-k_x,k_y).$$

$${\mathcal{C}}_y^{\prime -1}\mathcal{T}(E){\mathcal{C}}_y^{\prime }=\mathcal{T}(E).$$

$${\mathcal{C}}_y^{\prime -1}\Delta \mathcal{G}(E,x,y){\mathcal{C}}_y^{\prime }=\Delta \mathcal{G}(E,-x,y).$$

$$S_z(E,x,y)=-S_z(E,-x,y),$$

$$S_x(E,x,y)=+S_x(E,-x,y),$$

$$S_y(E,x,y)=-S_y(E,-x,y).$$

$$\mathcal{I}\Delta {\mathcal{G}}_0(E,x,y)\mathcal{I}=\Delta \mathcal{G}(E,-x,-y),$$

$$S_z(E,x,y)=+S_z(E,-x,-y),$$

$$S_x(E,x,y)=+S_x(E,-x,-y),$$

$$S_y(E,x,y)=+S_y(E,-x,-y).$$

$${\mathcal{C}}_z^{-1}\Delta \mathcal{G}(E,\mathbf{R}){\mathcal{C}}_z=\Delta \mathcal{G}(E,\tilde{\mathbf{R}}).$$

$$S_z(E,\mathbf{R})=S_z(E,\tilde{\mathbf{R}}),$$

$$\left(-\frac{1}{2}{S}_x+\frac{\sqrt{3}}{2}{S}_y\right)(E,\mathbf{R})=S_x(E,\tilde{\mathbf{R}}),$$

$$\left(-\frac{\sqrt{3}}{2}{S}_x-\frac{1}{2}{S}_y\right)(E,\mathbf{R})=S_y(E,\tilde{\mathbf{R}}).$$

$$E(\mathbf{k})=\pm \sqrt{\alpha \pm 2\sqrt{\beta }},$$

$$\text{where }\alpha$$

$$\beta$$

$$\gamma$$

Full text

I Introduction

The search for topological superconductors (TSCs) is an important topic in recent years Qi and Zhang (2011); Sato and Ando (2017). Nontrivial topology of the superconducting gap in TSC provides a buildable platform for the study of Majorana bound states, exotic quasiparticles which are their own antiparticles; obey Non-Abelian statistics; and have been considered as fundamental building blocks for future topological quantum computation Nayak et al. (2008). Unconventional superconductor Sr2RuO4 has been investigated for many years and been believed to be an intrinsic TSC with spin-triplet pairing Mackenzie and Maeno (2003). Recent study on iron-based superconductor shows that the surface states of FeTe0.55Se0.45 are intrinsic topological superconducting Zhang et al. (2018). In addition to these prominent candidates, another series of materials are investigated in recent years as potential TSCs, CuxBi2Se3 and its variants, SrxBi2Se3, NbxBi2Se3, TlxBi2Se3, etc.

The parent material of these candidates, Bi2Se3, is a strong topological insulator with bulk band gap about $0.3$ eV and gapless surface states protected by time-reversal symmetry Zhang et al. (2009); Liu et al. (2010). When the Fermi energy is tuned into the conduction or valence band by chemical doping, strong spin-orbit coupling may drive these materials into a topological nontrivial superconducting phase at low temperature. Experimental investigations show that Hor et al. (2010), if the copper atoms are intercalated in between the quintuple layers of Bi2Se3, the Fermi energy will be tuned into the conduction band, and there will be a superconducting transition at maximum $T_{C}=3.8$ K for copper concentration $0.12<x<0.15$. Theoretical investigations of Fu and Berg Fu and Berg (2010) show that, the order parameter of this superconductor should be parity-odd; belong to the 2-dimensional $E_{u}$ representation of the $D_{3}$ point group; and support topological nontrivial Andreev bound states. This kind of pairing potential, which preserves the twofold rotation symmetry but spontaneously breaks the threefold rotation symmetry, is generally named nematic pairing potential Fu (2014). Experimental investigations Sasaki et al. (2011); Kirzhner et al. (2012) show that there is an in gap zero-bias conductance peak (ZBCP) in the point contact spectra measurement on the surface of CuxBi2Se3. The existence of the ZBCP confirms that CuxBi2Se3 should be a parity-odd topological superconductor. However, detailed measurements based on scanning tunneling microscope (STM) Levy et al. (2013); Du et al. (2017a) and Andreev reflection spectroscopy Peng et al. (2013) give negative evidences. It is reported that the ZBCP disappears for clean materials. This pronounced peak is more likely induced by Andreev reflection between the point contact and the superconductor. It is further suggested that the pairing potential is $s$-wave symmetric Levy et al. (2013). On the other hand, various experiments support that threefold rotation symmetry is broken and that these materials are nematic pairing Fu (2016), i.e., nuclear magnetic resonance Matano et al. (2016); Nagai and Ota (2016), thermodynamics Yonezawa et al. (2016), upper critical field Yonezawa et al. (2016); Venderbos et al. (2016); Pan et al. (2016); Shen et al. (2017), penetration depth Smylie et al. (2016, 2017), torque magnetometry Asaba et al. (2017), superconducting gap versus in-plane magnetic field measurement Tao et al. (2018), magnetoresistance Yu et al. (2018), and transport experiment Du et al. (2017b), etc. In addition to these 3-dimensional materials, recently few layer Bi2Te3 grown on FeTe0.55Se0.45 substrate is demonstrated to be nematic pairing by using the quasiparticle interference (QPI) and magnetic flux vortex measurements Chen et al. (2018), etc. Despite the experimental evidences of nematic pairing listed above, different pairing potentials are proposed from theoretical considerations, e.g., momentum-dependent pairing potentials proposed in early stage Chen and Wan (2013), chiral pairing Wu and Martin (2017) belongs to the $E_{u}$ representation of $D_{3}$ point group, etc.

The identification of pairing symmetry in CuxBi2Se3 and its variants plays an critical role to distinguish the topological nontrivial nematic pairing and the trivial one. Experimentally, the phase-sensitive tetracrystal measurements Tsuei et al. (1994, 1997), the spacial distribution of Yu-Shiba-Rusinov resonance states induced by magnetic impurities Yu (1965); Shiba (1968); Rusinov (1969), and QPI spectrum Hoffman et al. (2002); Wang and Lee (2003); Hanaguri et al. (2007, 2009) are typical methods applied to identify the pairing symmetry of superconductors. In this work, we propose another method, the spin-polarized spectroscopy (SPS), to discern the nematic pairing. SPS of magnetic impurities in superconductors has been studied in previous works Kaladzhyan et al. (2016a, b) to measure the spin-orbital couplings in superconductors. Here we systematically study the SPS of magnetic impurities located on the surface of CuxBi2Se3 and Bi2Se${}_{3}/$FeTe0.55Se0.45 superconductors. We find that, for these materials with strong spin-orbit coupling, the SPS can give a clear manifestation of the resonance state and help us to distinguish the pairing potentials which belong to different representations of the $D_{3}$ point group.

The rest of this article is organized as follows: theoretical investigations are presented in Sec. II. Model Hamiltonian, the derivation of spin-polarized (SP) LDOS and QPI, and symmetry analysis of $D_{3}$ point group are given in this section. In Sec. III, we show the numerical results of SP LDOS for both the 2- and 3-dimensional systems. A conclusion is given in Sec. IV.

II Theoretical Investigation

II.1 Model Hamiltonian

We use the following two-orbital tight binding models Zhang et al. (2009); Fu (2009); Liu et al. (2010); Hao et al. (2014); Hao and Ting (2017) to describe the normal state of CuxBi2Se3 ($h_{3d}$) and Bi2Te3 grown on FeTe0.55Se0.45 ($h_{2d}$),

[TABLE]

$\psi^{{\dagger}}_{\bm{k}}$ and $\psi_{\bm{k}}$ are the creation and annihilation operators of electron states with wave-vector $\bm{k}$. The explicit expression of $\psi_{\bm{k}}$ is $(a_{\bm{k},\uparrow},a_{\bm{k},\downarrow},b_{\bm{k},\uparrow},b_{\bm{k},\downarrow})^{\mathsf{T}}$, here $a$ and $b$ refer to, respectively, the two low-energy $p_{z}$ orbitals of Se atoms on the top and bottom layers of each Bi2Se3 quintuple unit. $\Gamma_{0}$ to $\Gamma_{5}$ are $4\times 4$ matrices acting in the orbital and spin space. $\Gamma_{0}=\sigma_{0}\otimes{s_{0}}$ is the identity matrix, $\Gamma_{1}=\sigma_{z}\otimes{s_{x}}$, $\Gamma_{2}=\sigma_{z}\otimes{s_{y}}$, $\Gamma_{3}=\sigma_{z}\otimes{s_{z}}$, $\Gamma_{4}=-\sigma_{y}\otimes{s_{0}}$, and $\Gamma_{5}=\sigma_{x}\otimes{s_{0}}$. $\sigma_{\#}$ and $s_{\#}$ (${\#}=x,y,z$) are Pauli matrices in the orbital and spin spaces, respectively. $\varepsilon(\bm{k})$ represents the normal dispersion, the off-diagonal terms proportional to $M(\bm{k})$ represent the band inversion of topological insulator, the terms with $A_{0}$ and $B_{0}$ contained represent the in-plane and out of-plane spin-orbit couplings, the last two terms in Eqs. (2) and (3) induce the hexagonal warping effect of the Fermi surface. For completeness, here we list the details of these terms represented in tight-binding model. $\varepsilon(\bm{k})=C_{0}+2C_{1}[1-\cos(\bm{k}\cdot\bm{\delta}_{4})]+\frac{4}{3}C_{2}[3-\cos(\bm{k}\cdot\bm{\delta}_{1})-\cos(\bm{k}\cdot\bm{\delta}_{2})-\cos(\bm{k}\cdot\bm{\delta}_{3})]$, $M(\bm{k})=M_{0}+2M_{1}[1-\cos(\bm{k}\cdot\bm{\delta}_{4})]+\frac{4}{3}M_{2}[3-\cos(\bm{k}\cdot\bm{\delta}_{1})-\cos(\bm{k}\cdot\bm{\delta}_{2})-\cos(\bm{k}\cdot\bm{\delta}_{3})]$, $c_{x}(\bm{k})=\frac{1}{\sqrt{3}}[\sin(\bm{k}\cdot\bm{\delta}_{1})-\sin(\bm{k}\cdot\bm{\delta}_{2})]$, $c_{y}(\bm{k})=\frac{1}{3}[\sin(\bm{k}\cdot\bm{\delta}_{1})+\sin(\bm{k}\cdot\bm{\delta}_{2})-2\sin(\bm{k}\cdot\bm{\delta}_{3})]$, $c_{z}(\bm{k})=\sin(\bm{k}\cdot\bm{\delta}_{4})$, $d_{1}(\bm{k})=-\frac{8}{3\sqrt{3}}[\sin(\bm{k}\cdot\bm{a}_{1})+\sin(\bm{k}\cdot\bm{a}_{2})+\sin(\bm{k}\cdot\bm{a}_{3})]$, $d_{2}(\bm{k})=-8[\sin(\bm{k}\cdot\bm{\delta}_{1})+\sin(\bm{k}\cdot\bm{\delta}_{2})+\sin(\bm{k}\cdot\bm{\delta}_{3})]$. Here $\bm{\delta}_{j}$ and $\bm{a}_{j}$ represent the nearest-neighboring and in-plane next nearest-neighboring bond vectors, respectively. Their detailed expressions are, $\bm{\delta}_{1}=(\frac{\sqrt{3}}{2}a,\frac{1}{2}a,0)$, $\bm{\delta}_{2}=(-\frac{\sqrt{3}}{2}a,\frac{1}{2}a,0)$, $\bm{\delta}_{3}=(0,-a,0)$, $\bm{\delta}_{4}=(0,0,c)$, $\bm{a}_{1}=\bm{\delta}_{1}-\bm{\delta}_{2}$, $\bm{a}_{2}=\bm{\delta}_{2}-\bm{\delta}_{3}$, $\bm{a}_{3}=\bm{\delta}_{3}-\bm{\delta}_{1}$. $a=4.14\textup{\AA}$ and $3c=28.64\textup{\AA}$ are the in-plane and out-of-plane lattice spacings of CuxBi2Se3; for Bi2Te3 grown on FeTe0.55Se0.45 substrate, $a=4.38\textup{\AA}$. The parameters used for numerical calculation are listed in Table 1. These Hamiltonians, Eqs. (2) and (3), preserve the desired rotation symmetry of $D_{3d}$ point group. $D_{3d}$ point group has 12 elements which are generated by 3 elementary operations: threefold rotation along the $z$-axis, twofold rotation along the $y$-axis, and inversion. The matrix representations of these elementary operations are given in Table 2 for both the normal state and superconducting cases.

For both the 3- and 2-dimensional situations, the superconducting Hamiltonian can be written in the Bogliubov-de Gennes (BdG) formalism as follows by using the Nambu spinor basis $\Psi^{{\dagger}}_{\bm{k}}=(\psi^{{\dagger}}_{\bm{k}},\psi_{-\bm{k}}^{\mathsf{T}})$,

[TABLE]

where $\Delta$ is the superconducting pairing potential, $h(\bm{k})$ equals to $h_{3d}(\bm{k})$ and $h_{2d}(\bm{k})$ for the 3- and 2-dimensional systems respectively. In this work, we focus on the pairing potentials belong to the 2-dimensional $E_{u}$ representation of the $D_{3d}$ point group (the nematic pairing). As shown in previous studies Fu and Berg (2010); Hao et al. (2014); Hao and Ting (2017), these pairing potentials in the $E_{u}$ representation can be written as, $\Delta_{4a}=i\Delta_{a}\sigma_{y}\otimes{s_{0}}$ and $\Delta_{4b}=\Delta_{b}\sigma_{y}\otimes{s_{z}}$, with $\Delta_{a}$ and $\Delta_{b}$ referring to the strength of pairing potentials. Generally, these pairing potentials break the threefold rotation symmetries of $D_{3d}$ point group.

II.2 SP LDOS and QPI

The retarded Green’s function of the BdG Hamiltonian can be evaluated using the standard method, which gives,

[TABLE]

We assume that there is an impurity located at the original point, i.e., $\hat{V}=\delta(\bm{R})\psi^{{\dagger}}(\bm{R})\mathcal{V}\psi(\bm{R})$. For nonmagnetic and magnetic impurities, the impurity potential take the following forms,

[TABLE]

the superscripts, NM and M, refer to the nonmagnetic and magnetic impurity potentials. $J_{0}$ and $\bm{J}=(J_{x},J_{y},J_{z})$ represent the strength of scattering potential for nonmagnetic and magnetic impurities, respectively. Here, we have assumed that the impurity potential is short-ranged.

In the following studies, we consider two different configurations: (1) for 3-dimensional material CuxBi2Se3, described by the Hamiltonian (2), the impurity is located on the surface of the material, i.e., the $xy$-plane at $z=0$; (2) for few layer Bi2Te3 grown on FeTe0.55Se0.45 substrate, described by the 2-dimensional Hamiltonian (3), the impurity is located at the 2-dimensional plane. For both these two cases, $\bm{R}=(x,y)$.

When the Kondo temperature is much smaller than the superconducting gap, the magnetic impurity is not screened by the itinerant electrons, and it can be considered as classical impurity Balatsky et al. (2006). In this scenario, the full Green’s function of the electrons affected by the impurity can be solved by using the standard $T$-matrix method. For the short-ranged impurity potential given in Eqs. (6) and (7), the $T$-matrix is expressed as,

[TABLE]

where $\mathcal{G}_{0}(E,\bm{0})$ is the free Green’s function at $\bm{R}=(0,0)$. For the 2-dimensional Hamiltonian, the non-perturbed retarded Green’s function, $\mathcal{G}_{0}(E,\bm{R})$ is evaluated as the Fourier transformation of Eq. (5),

[TABLE]

where $N$ is the number of wave-vectors in summation. For the 3-dimensional Hamiltonian (the impurity is located at the original point on the surface $z=0$), $\mathcal{G}_{0}(E,\bm{R})$ is the non-perturbed retarded Green’s function of surface states, which can be derived numerically by using the fast iterative method Lopez Sancho et al. (1984). The full Green’s function in real space is given by,

[TABLE]

The SP LDOS is defined as,

[TABLE]

where $\Delta\mathcal{G}_{ij}$ takes the $i$th row and $j$th column component of the matrix $\Delta\mathcal{G}(E,\bm{R})$. Generally, the orientation of magnetic impurity and SP LDOS can take various combinations, e.g., $\bm{J}=(J,0,0)$ (the magnetic impurity is oriented in the $x$-direction) and $S_{z}(E,\bm{R})$ (the SP LDOS of $z$-component is detected).

In addition to the SP LDOS, the QPI is another important method to identify the pairing symmetry of superconductors, which is defined as follows for nonmagnetic impurity,

[TABLE]

Here we need to notify that, generally, $\delta\rho(E,\bm{q})$ and the Fourier transformation of SP LDOS, $S_{\#}(E,\bm{q})=\int\mathrm{d}\bm{R}e^{-i\bm{q}\cdot\bm{R}}S_{\#}(E,\bm{q})$ ($\#=x,y,z$) should take complex numbers, though $\delta\rho(E,\bm{R})$ and $S_{\#}(E,\bm{R})$ are real functions.

II.3 Symmetry analysis

In this subsection, we analyze the symmetries of different SP LDOS patterns if the pairing potential is invariant under some operations of the $D_{3}$ point group. Generally, the normal state Hamiltonian preserves all the symmetries listed in Table 2. Due to the strong spin-orbit coupling in these materials, if the pairing potential breaks some symmetries listed in Table 2, the SP LDOS patterns should reveals these symmetry-breakings.

For the 3-dimensional case, the situation is a bit complicated. We consider only the STM measurable resonance states, i.e., the impurity should be located on the surface of CuxBi2Se3. In this case, the inversion and twofold rotation symmetry related to $z\rightarrow-z$ is broken. Only the threefold rotation symmetry may be survival. In the following symmetry analysis, we focus on the 2-dimensional situation. For both 2- and 3-dimensional situations, numerical investigations and analysis of SP LDOS will be presented in the next section.

Firstly, we consider the twofold rotation operation along the $y$-axis, $\mathscr{C}_{y}^{\prime}$ given in Table 2. Suppose a pairing potential $\Delta$ preserves this symmetry, i.e.,

[TABLE]

In this case, the BdG Hamiltonian should be invariant under this operation,

[TABLE]

here we have used the explicit expression of wave vector $\bm{k}=(k_{x},k_{y})$, and

[TABLE]

Using the definition of unperturbed Green’s function, Eq. (5), we get,

[TABLE]

If the orientation of magnetic impurity takes the form $\bm{J}=(J,0,0)$, the impurity potential $\mathcal{V}$ should be invariant under the $\mathscr{C}^{\prime}_{y}$ operation, so we get $\mathscr{C}_{y}^{\prime-1}\mathcal{V}\mathscr{C}_{y}^{\prime}=\mathcal{V}$, and,

[TABLE]

Substituting Eqs. (20) and (21) into (9) and (11), we find,

[TABLE]

where $(x,y)$ is the explicit expression of location $\bm{R}$. Using the definition of SP LDOSs given in Eqs. (12)-(14), we can find that,

[TABLE]

If the pairing potential is invariant under the inversion operation $\mathscr{I}$ given in Table 2, after a similar derivation, we find,

[TABLE]

and

[TABLE]

for arbitrary impurity spin orientation $\bm{J}=(J_{x},J_{y},J_{z})$.

Finally, if the pairing potential is invariant under the threefold rotation operation $\mathscr{C}_{z}$ given in Table 2 and the spin of magnetic impurity is oriented at $z$-direction, $\bm{J}=(0,0,J)$, straightforward calculation gives,

[TABLE]

where $\tilde{\bm{R}}=(-\frac{1}{2}x-\frac{\sqrt{3}}{2}y,\frac{\sqrt{3}}{2}x-\frac{1}{2}y)$ for $\bm{R}=(x,y)$. Using the expressions of SP LDOS, we get,

[TABLE]

The eigenvalue of twofold rotation ($\mathscr{C}_{y}^{\prime}$) and inversion ($\mathscr{I}$) operations are both $\pm 1$. For a given pairing potential which is parity-odd under these operations and the eigenvalue is $-1$, i.e., $\Delta{\xLongrightarrow{\mathscr{C}_{y}^{\prime}{\rm~{}or~{}}\mathscr{I}}}-\Delta$, we find that, the results (23)-(25) and (27)-(29) are still correct. For the two species of pairing potentials belong to the $E_{u}$ representation of $D_{3}$ point group, we find the following relationships, $\Delta_{4a}\xLongrightarrow{\mathscr{C}_{y}^{\prime}}-\Delta_{4a}$, $\Delta_{4b}\xLongrightarrow{\mathscr{C}_{y}^{\prime}}+\Delta_{4b}$, $\Delta_{4a}\xLongrightarrow{\mathscr{I}}-\Delta_{4a}$, $\Delta_{4b}\xLongrightarrow{\mathscr{I}}-\Delta_{4b}$. Therefor, the results (23)-(25) and (27)-(29) will be manifested if the magnetic impurity is properly oriented and the pairing potential is nematic. However, if the pairing potential is a finite combination of $\Delta_{4a}$ and $\Delta_{4b}$, e.g., the chiral pairing proposed in Wu and Martin (2017), the twofold rotation symmetry and the inversion symmetry are generally broken and the pattern preserves Eqs. (23)-(25) and (27)-(29) will not be observed. On the other hand, the nematic pairing $\Delta_{4a}$ or $\Delta_{4b}$ breaks the threefold rotation symmetry straightforwardly, so the symmetry of SP LDOS given in Eqs. (31)-(33) will be violated for $\Delta_{4a}$ and $\Delta_{4b}$.

III Numerical results

It is difficult to get the analytic expressions of the SP LDOSs for the Hamiltonians given in Eqs. (1)-(3). In order to get a more intuitive understanding and show the quantitative results of the SP LDOS influenced by the pairing symmetry, in this section, we present the numerical results of SP LDOS for the parameters given in Table 1. Before immersing into the detailed discussion of the numerical results, here we give some more details of the numerical calculations. The imaginary part of retarded Green’s function has been chosen to be $0.001$eV. A mesh size with $1024\times 1024$ points has been chosen for the SP LDOS and QPI calculations. In the evaluation of SP LDOS versus energy shown in Figs. 2(g) and 3(g), a mesh size up to $4096\times 4096$ points has been chosen to convergence the results. The fast Fourier transformation is applied to speed up the numerical calculation.

Firstly, we consider the 2-dimensional situation, few layer Bi2Te3 grown on FeTe0.55Se0.45 substrate. The four (twofold degenerated) eigenenergies of the BdG Hamiltonian, Eq. (4), are,

[TABLE]

$E_{F}$ is the Fermi energy, $\Delta_{\#}$ ($c_{\#}$) equals to $\Delta_{a}$ ($c_{y}(\bm{k})$) and $\Delta_{b}$ ($c_{x}(\bm{k})$) for $\Delta_{4a}$ and $\Delta_{4b}$ pairings, respectively. For each pairing potential ($\Delta_{4a}$ and $\Delta_{4b}$), the eigenenergies close to the Fermi surface have two nodes. Numerical calculation finds that these nodes are located at $\bm{k}a\approx(\pm 0.26,0)$ and $(0,\pm 0.26)$ for $\Delta_{4a}$ and $\Delta_{4b}$ pairings, respectively. These nodes make the density of states (DOS) to be V-shape like (see Fig. 1(a) for details), which is consistent with the experimental results shown in Ref. Chen et al. (2018) for few quintuple layers of Bi2Te3 grown on FeTe0.55Se0.45 substrate.

Now we analysis the QPI shown in Fig. 1. In order to get a clear representation, the QPI is plotted in a small region near the $\Gamma$ point. When the energy is located on the Fermi surface, i.e. $E=0$ shown in Fig. 1(d), the QPI manifests strong anisotropy along the $x$-direction, which indicates that the pairing potential $\Delta_{4a}$ has two nodes on the $k_{x}$-axis. Typical momentum transfer $\bm{q}a\approx(\pm 0.52,0)$ shown in Fig. 1(d) is consistent with the position of nodes $\bm{k}a\approx(\pm 0.26,0)$. When the energy is tuned away from the superconducting gap, i.e., $E=-0.03$eV and $0.03$eV shown in Figs. 1(c) and (f), QPI reveals sixfold rotation symmetry, which reflects the fact that the lattice structure is threefold rotation symmetric along the $z$-axis.

Next we study the SP LDOS for a magnetic impurity with spin oriented at $x$-direction, $\bm{J}=(J,0,0)$. Here, the nematic pairing symmetry $\Delta_{4a}$ are chosen for study. Fig. 2(g) shows the spatially averaged SP LDOS which is defined as $S_{x}(E)=S_{x}(E,\bm{q}=0)$. From this plot, one can find two in-gap Yu-Shiba-Rusinov (YSR) resonance states at $E\approx\pm{1.05}$meV. For the negative (positive) resonance energy $E=-1.05$meV ($E=1.05$meV), the resonance state and the magnetic impurity is ferromagnetic (antiferromagnetic) coupled. Figs. 2(a)-(c) show the spatial distributions of the YSR resonance states for different spin orientations respectively. From these plots, one can find that the SP LDOS exist Friedel oscillations with oscillation period on the scale of coherence length. The constraints from twofold rotation symmetry and inversion symmetry are manifested in these plots: (1) $S_{x}$ is invariant under the twofold rotation operation $\mathscr{C}_{y}^{\prime}$, $(x,y)\rightarrow(-x,y)$, and the inversion operation $\mathscr{I}$, $(x,y)\rightarrow(-x,-y)$; (2) both $S_{y}$ and $S_{z}$ change signs after the twofold rotation operation $\mathscr{C}_{y}^{\prime}$; and (3) $S_{y}$ and $S_{z}$ are invariant under the inversion operation $\mathscr{I}$. The combination of twofold rotation and inversion operations give another strong constraint that $S_{y}$ and $S_{z}$ should be vanishing on the $x$- and $y$-axes. These feathers are more significant after the Fourier transformation, the numerical results are presented in Figs. 2(d)-(f), respectively.

Next we study the SP LDOS for another orientation of the magnetic impurity, $\bm{J}=(0,0,J)$. The spatially averaged SP LDOS $S_{z}(E)=S_{z}(E,\bm{q}=0)$ as a function of energy shown in Fig. 3(g) is qualitatively consistent with that given in Fig. 2(g): there are two in-gap YSR resonance peaks at $E\approx\pm{12}$meV respectively; for the negative (positive) resonance energy, the YSR state is ferromagnetic (antiferromagnetic) coupled to the impurity. Now we study the spacial distribution of the YSR resonance state with negative resonance energy, $E\approx-12$meV. Figs. 3(a)-(c) show the contour plots of SP LDOS of $S_{z}(E,\bm{R})$, $S_{x}(E,\bm{R})$ and $\left(-\frac{1}{2}S_{x}+\frac{\sqrt{3}}{2}S_{y}\right)(E,\bm{R})$, respectively. Generally, as discussed in the previous section, if the pairing potential preserves the threefold rotation symmetry, $S_{z}$ should be invariant under an 120∘ rotation operation. However, both the real space SP LDOS 3(a) and its Fourier transformation 3(d) show a clear 180∘ rotation symmetry. The threefold rotation symmetry broken is evident. Furthermore, as shown in Eq. (32), if the pairing potential is invariant under the threefold rotation operation, plot 3(b) will coincide to 3(c) after an 120∘ rotation operation. It is significant that these two plots can not coincide to each other. These observations demonstrate that SP LDOS can be applied to identify the threefold rotation symmetry broken if the pairing potential is nematic.

Finally, we study the SP LDOS for a magnetic impurity located on the surface of CuxBi2Se3. Pairing potential has been chosen to be $\Delta_{4a}$ with pairing strength $\Delta_{a}=0.02$eV. The numerical results are presented in Fig. 4. As investigated in previous studies, due to the nontrivial topology of the pairing potential, there exists robust surface Andreev boundary states (SABSs) in the superconducting gap. Figs. 4(m) and (n) show the spacial averaged SP LDOS, $S_{z}(E)$ and $S_{x}(E)$ versus energy, for magnetic $z$-impurity and $x$-impurity respectively. In these plots, we find that the zero-energy SABSs strongly coupled to the magnetic impurity antiferromagnetically. As shown in Fig. 4(n), for the $x$-impurity, the spacial averaged SP LDOS is complicated. In addition to the antiferromagnetic (AFM) YSR resonance peak at $E=0$, there are two other resonance peaks near the band edge $E=15$meV. Here we focus on the zero-energy YSR resonance states. The first and third rows show the spacial distributions of the zero-energy YSR resonance states induced by $z$-impurity and $x$-impurity, respectively. The second and forth rows present the Fourier transformation of these SP LDOSs correspondingly. Inversion symmetry, $(x,y,z)\rightarrow(-x,-y,-z)$, is spontaneously broken for the surface impurity problem. The spatial distributions of the resonance states are dramatically different from the 2-dimensional case. One can find that Figs. 4(a), (b) and (i) are almost antisymmetric under the surface inversion operation $(x,y)\rightarrow(-x,-y)$. Numerical investigations show that the Fourier transformation of these SP LDOSs have both real and imaginary parts; the real parts are about two magnitudes smaller than the imaginary parts, respectively for each sub-figure. Here we give only the imaginary parts of SP LDOSs, which are plotted in Figs. 4(d), (e) and (l). Due to the intrinsic inversion symmetry broken of the surface states, it is difficult to identify the inversion symmetry of the bulk pairing potentials. Here we focus on the threefold symmetry broken of pairing potential. If the pairing potential and the impurity orientation is threefold rotation symmetric along the $z$-axis, i.e., $\bm{J}=(0,0,J)$, $S_{z}(E,\bm{R})$ and its Fourier transformation should be threefold rotation symmetric. One can find from the Fourier transformation Fig. 4(f) that $S_{z}(E,\bm{q})$ reveals twofold rotation symmetry. This symmetry broken is not as obvious as the 2-dimensional case. This may be induced by the topological surface states (TSSs). For the 3-dimensional material with a small Fermi surface, the TSSs of the parent material Bi2Se3 coexist with the SABSs Wray et al. (2010); Hao and Lee (2011); Hao and Wang (2015). Generally, there are three kinds of scatterings: the scattering between TSSs, the scattering between TSSs and SABSs, and the scattering between SABSs. The first scattering preserves the threefold rotation symmetry. Only the scattering between SABSs can significantly reveal the symmetry broken if the pairing potential is nematic. This signature may be weakened by the inconsequential scatterings.

IV Conclusion

We have analyzed the SPS of magnetic impurity resonance states for 2- and 3-dimensional nematic superconductors, in particular, few layer Bi2Se3 grown on FeTe0.55Se0.45 (2D) and CuxBi2Se3 (3D) respectively. We focus on the relationship between nematic pairing symmetry and the topograph of SP LDOS. Different constraints from different symmetries are investigated, i.e., the SP LDOSs and their Fourier transformations should satisfy some constraints if the pairing potential preserves inversion symmetry, twofold rotation symmetry along the $y$-axis, or threefold rotation symmetry along the $z$-axis. The spacial distributions and the Fourier transformations of the SP LDOSs of YSR resonance states are investigated. For 2-dimensional case (few layer Bi2Se3 grown on FeTe0.55Se0.45), the SP LDOSs and their Fourier transformations reveal significant twofold rotation symmetry, threefold rotation symmetry is evident broken. For 3-dimensional case (CuxBi2Se3), threefold rotation symmetry breaking can be identified by studying the Fourier transformation of the SP LDOS $S_{z}(E,\bm{R})$ for $z$-impurity orientation. In our studies, we use the tight-binding model with realistic parameters fitting the first principle calculations, our results can be tested experimentally.

In this work, we present only the numerical results for one species of the $E_{u}$ representation of the nematic pairing potential, $\Delta=\Delta_{4a}$. For the 2-dimensional case, spatial distribution and corresponding Fourier transformation of resonance states with only negative resonance energy are presented. The numerical calculations for other situations have been accomplished, i.e., pairing potential is chosen to be $\Delta_{4b}$ and the SP LDOS of resonance states with positive resonance energy. We did not find qualitative difference, so the numerical results are not presented here.

Acknowledgment

We appreciate support from the NSFC under Grants No. 11504106 and No. 11447167 and the Fundamental Research Funds for the Central Universities under Grant No. 2018MS049.