Nematic superconductivity in Cu$_{x}$Bi$_{2}$Se$_{3}$: The surface Andreev bound states
Lei Hao C. S. Ting








 ^1Department of Physics and Texas Center for Superconductivity, University of Houston, Houston, Texas 77204, USA 

 ^2Department of Physics, Southeast University, Nanjing 210096, China 





We study theoretically the topological surface states (TSSs) and the possible surface Andreev bound states (SABSs) of Cu_xBi_2Se_3 which is known to be a topological insulator at x=0. The superconductivity (SC) pairing of this compound is assumed to have the broken spin-rotation symmetry, similar to that of the A-phase of ^3He as suggested by recent nuclear-magnetic resonance experiments. For both spheroidal and corrugated cylindrical Fermi surfaces with the hexagonal warping terms, we show that the bulk SC gap is rather anisotropic; the minimum of the gap is negligibly small as comparing to the maximum of the gap. This would make the fully-gapped pairing effectively nodal. For a clean system, our results indicate the bulk of this compound to be a topological superconductor with the SABSs appearing inside the bulk SC gap. The zero-energy SABSs which are Majorana fermions, together with the TSSs not gapped by the pairing, produce a zero-energy peak in the surface density of states (SDOS). The SABSs are expected to be stable against short-range nonmagnetic impurities, and the local SDOS is calculated around a nonmagnetic impurity. The relevance of our results to experiments is discussed.





Nematic superconductivity in Cu_xBi_2Se_3: The surface Andreev bound states
    Lei Hao^1,2 and C. S. Ting^1
    December 30, 2023
===========================================================================






§ INTRODUCTION




Cu_xBi_2Se_3, the first superconductor (SC) realized in a doped topological insulator, has attracted great interest since its discovery <cit.>. After early gathering of experimental and theoretical evidences trying to imply the superconductivity (SC) pairing in this compound to be topologically nontrivial <cit.>, several later experiments seem to conclude that the actual pairing symmetry is the conventional s-wave and topologically trivial <cit.>. However, initiated by a recent nuclear magnetic resonance (NMR) experiment <cit.>, a new surge of researches on this material and its several variants (e.g., Sr_xBi_2Se_3 <cit.>, Nb_xBi_2Se_3 <cit.>, and Tl_xBi_2Se_3 <cit.>) revived the possibility that the SC  phase realized in this group of materials is topologically nontrivial.
In the NMR experiment, a prominent in-plane uniaxial anisotropy in the SC order parameter is observed in the Knight shift measurement <cit.>. The above twofold in-plane rotational symmetry is confirmed further by field-angle dependent specific heat measurement <cit.>. These new experiments indicate that the SC pairing in Cu_xBi_2Se_3 might be the first example of a pairing breaking spontaneously the spin rotation symmetry of the parent material, from the threefold rotational symmetry of the normal phase to the uniaxial twofold in-plane symmetry <cit.>.




Amazingly, the new NMR experiment is claimed to be explainable by an odd-parity pairing proposed earlier <cit.>. The alluded pairing, with a salient in-plane anisotropy, has also been called a nematic SC <cit.>. For a simplified model with spheroidal Fermi surface, this SC pairing is known to be equivalent to the A phase of ^3He  and has two bulk nodes <cit.>. However, by including in the model the terms responsible for hexagonal warping of the Fermi surface, this SC was argued to have a full pairing gap <cit.>. If confirmed, this could be the first three-dimensional topological superconductor with a fully-gapped bulk ever discovered. However, a hallmark of the topological SC is the presence of surface Andreev bound states (SABSs) within the bulk SC gap, and its robustness against nonmagnetic impurities and interactions. In early studies in terms of simplified models based on spheroidal Fermi surface without the hexagonal warping terms, the SC is known to support novel SABSs which are flat along one direction [i.e., (k_x,0)] and linearly dispersive in the perpendicular direction [i.e., (0,k_y)] <cit.>. In addition, the topological surface states (TSSs) are well defined at the chemical potential and well separated from the bulk states<cit.>. The SC pairing would not open a gap in the TSSs and thus the TSSs coexist with the SABSs <cit.>. If the hexagonal warping term is incorporated, then on one hand the two bulk nodes would be gapped out <cit.>, whereas on the other hand the TSSs will remain gapless <cit.>. In addition, the odd-parity topological SC pairing should still supports SABSs, independent of the model parameters <cit.>. As a result, it seems natural to expect the existence of the nontrivial spectral features related to the two types of surface states.



Previous experiments, on the other hand, have made conflicting statements on the surface states in the SC phase. In several point contact spectroscopy (PCS) studies, a pronounced zero-bias peak appeared and was ascribed to the existence of SABSs <cit.>. Later, a scanning tunneling spectroscopy (STS) study reported a standard BCS-like spectrum <cit.>. Moreover, a detailed PCS experiment arrived at the same conclusion that no SABSs existed <cit.>. However, recent experiments clearly indicated that the Fermi surface of the Cu_xBi_2Se_3 compound changes from spheroidal to corrugated-cylindrical surfaces as the doping x increases <cit.>. As a result and without detailed investigations, there exist works <cit.> arguing that the absence of SABSs in the odd-parity (e.g., nematic) SC state of Cu_xBi_2Se_3 is consistent with a (corrugated) cylindrical Fermi surface.




Inspired by the above experimental observations and theoretical arguments, in this work we explore whether the corrugated cylindrical Fermi surface would or would not support the SABSs in the Cu_xBi_2Se_3 compound with the nematic SC pairing proposed for explaining the recent experiments <cit.>. Using the band parameters which fit qualitatively the experimental Fermi surfaces and the TSSs, we find that the SC is fully gapped, and the bulk gap is rather anisotropic. The minima of the bulk SC gap is smaller than the maximum of the bulk gap by two to three orders of magnitude. The bulk quasiparticle spectrum, while in principle is fully gapped, appears to be nodal-like from the point of view of measurements. The SABSs are shown to exist in the clean system for both spheroidal and corrugated cylindrical Fermi surfaces, and the zero-energy Majorana bound state is a characteristic of the SABSs. In addition, we verify the stability of the SABSs against short-range nonmagnetic impurities, both for uniformly distributed bulk impurities and for dilute impurities doped only to the surface layer. As to whether the SABSs could be observed experimentally may depend on the condition of the sample surfaces. For instance, the excessive magnetic Cu^2+ (3d^9) ions or Cu (3d^104S^1) atoms on the surfaces could very much suppress the SABSs. For clean and perfect surface, the SABSs should be detectable. Recently, the possible existence of surface states in superconducting Sr_xBi_2Se_3 is inferred from the Shubnikov-de Hass oscillation measurement <cit.>. In combination with the results of the present work, it is therefore highly desirable to reexamine the existence of the SABSs in Cu_xBi_2Se_3 as a crucial test for the relevancy of the proposed SC pairing <cit.>.







§ MODEL


We consider a two-orbital tight-binging model for the low-energy degrees of freedom of the material

    H_0(𝐤)   =   ϵ(𝐤)I_4+M(𝐤)Γ_5+B_0c_z(𝐤)Γ_4
    +A_0[c_y(𝐤)Γ_1
          -c_x(𝐤)Γ_2]+R_1d_1(𝐤)Γ_3+R_2d_2(𝐤)Γ_4.

The basis vector is taken as ϕ^†_𝐤=[a^†_𝐤↑, a^†_𝐤↓, b^†_𝐤↑, b^†_𝐤↓], where the two orbitals (a and b) are mainly from the two p_z orbitals on the top and bottom Se layers of each Bi_2Se_3 quintuple unit <cit.>. I_4 is the 4×4 unit matrix. Γ_1=σ_3⊗ s_1, Γ_2=σ_3⊗ s_2, Γ_3=σ_3⊗ s_3,
Γ_4=-σ_2⊗ s_0, and Γ_5=σ_1⊗ s_0 <cit.>.
s_i and σ_i (i=1,2,3) are Pauli matrices for the spin and orbital degrees of freedom. The parity operator is defined as P=σ_1⊗ s_0 <cit.>. The above model was obtained previously based on symmetry and comparison with an existing 𝐤·𝐩 model <cit.>. The lattice of Bi_2Se_3 and Cu_xBi_2Se_3, which belong to the D_3d^5 space group, is mapped to a hexagonal lattice in the tight-binding model. The in-plane (labeled as the xy plane) and out-of-plane (labeled as the z direction) lattice parameters, a and c, are taken as a=4.14 Å and 3c=28.64 Å <cit.>. ϵ(𝐤)=C_0+2C_1[1-cos(𝐤·δ_4)]
+4/3C_2[3-cos(𝐤·δ_1)-cos(𝐤·δ_2)
-cos(𝐤·δ_3)]. M(𝐤) is obtained from ϵ(𝐤)
by making the substitutions C_i→ M_i (i=0,1,2). c_x(𝐤)=1/√(3)
[sin(𝐤·δ_1)-sin(𝐤·δ_2)],
c_y(𝐤)=1/3[sin(𝐤·δ_1)+sin(𝐤·δ_2)
-2sin(𝐤·δ_3)], and c_z(𝐤)=sin(𝐤·δ_4).
Finally, d_1(𝐤)=-8/3√(3)[sin(𝐤·𝐚_1)+sin(𝐤·𝐚_2)
+sin(𝐤·𝐚_3)] and d_2(𝐤)=-8[sin(𝐤·δ_1)
+sin(𝐤·δ_2)+sin(𝐤·δ_3)]. Here, the four nearest-neighboring bond vectors of the hexagonal lattice are δ_1=(√(3)/2a, 1/2a, 0),
δ_2=(-√(3)/2a, 1/2a, 0), δ_3=(0, -a, 0), and δ_4=(0, 0, c).
The three in-plane second-nearest-neighboring bond vectors in d_1(𝐤) are 𝐚_1=δ_1-δ_2,
𝐚_2=δ_2-δ_3, and 𝐚_3=δ_3-δ_1. The last two terms in H_0(𝐤) induce hexagonal warping of the Fermi surface and the topological surface states (TSSs)<cit.>.




Before doping with copper, the Fermi surface of Bi_2Se_3 is spheroidal. After intercalating copper to inter-quintuple-layer positions, the material becomes more two-dimensional. According to the experiments <cit.>, the Fermi surface for certain Cu_xBi_2Se_3 becomes (corrugated) cylindrical, although the details of the evolution are still unclear. On the other hand, a common feature of the normal phase of superconducting Cu_xBi_2Se_3 is that the TSSs are well defined and coexist with the Fermi surface <cit.>. In this work, we consider three sets of parameters shown in Table I. The Fermi surface contours (on the k_y=0 plane) and the surface spectral functions for the three sets of parameters are shown in Figure 1. The surface spectral functions are calculated in terms of the iterative Green's function method <cit.>, for the upper xy surface of a sample that can be regarded as consisting of an infinite number of layers. The parameters are chosen here to fit qualitatively three different shapes of the Fermi surfaces and the coexisting TSSs. The second corrugated cylindrical Fermi surface is less corrugated compared to the first corrugated cylindrical Fermi surface. Therefore, the material described by the `Cylindrical 2' is more two-dimensional than the material described by the `Cylindrical 1'. With the three typical sets of parameters, we can study the qualitative evolution of the property of a pairing as the Fermi surface turn from spheroidal to corrugated cylindrical and then becomes even more two-dimensional.



Note that, these parameter sets are chosen to reflect the evolution of the Fermi surface and the coexistence with the TSSs, which are most crucial for the low-energy physics in the superconducting phase. A completely two-dimensional model with zero hopping along the z direction is unsuitable because it cannot give the TSSs observed in experiments <cit.>. In addition, the relative magnitudes of the various parameters in Table I are in agreement with the set of parameters obtained previously by fitting the first-principle band structures for Bi_2Se_3 <cit.>. By increasing the value of R_1 artificially (e.g., to 2 eV) and keeping other parameters unchanged, the topology of the Fermi surface and the coexisting TSSs can still be retained qualitatively. We will discuss the effects of increasing R_1 at the end of Section V-B.
















§ PSEUDOSPIN BASIS




The full model contains the complete information of the phase, but it is hard to work with analytically. On the other hand, it is the states close to the Fermi level that are most important to the superconducting phase. By introducing the pseudospin basis, the full model containing both of the two bands of the model in the normal phase can be projected to a simplified model containing only the band contributing to the Fermi surface  <cit.>. By making this reduction, the low-energy properties of the superconducting phase, in particular the gap structure of the bulk quasiparticle spectrum and the SABSs, can be understood more easily. Here, we follow the approach of Yip, which was originally applied to a simplified version of the model, to construct the pseudospin basis for our tight-binding model <cit.>. This method makes use of the time-reversal symmetry (T) and inversion symmetry (P) of the model, which lead to the Kramers degeneracy of each state. The two pseudospin bases for each Kramers doublet are thus required to be related by the combined action of PT operation. Throughout this work, we assume the chemical potential to lie within the conduction band. The eigenbasis can be constructed by first diagonalizing the model in the spin subspace and then in the orbital subspace. One basis can be taken as

    |𝐤,α'⟩=1/D̃_𝐤N_𝐤[     Ẽ_𝐤; M̃_-(𝐤) ][ A_0c_+(𝐤);    D_-(𝐤) ],

where the first and second vectors are separately spinors in the subspaces of the original orbital and spin degree of freedom. For notational simplicity, here and later in this work we will use the following abbreviations
c_±(𝐤)=c_y(𝐤)± i c_x(𝐤), M̃_±(𝐤)=M(𝐤)± i[B_0c_z(𝐤)+R_2d_2(𝐤)], D_𝐤=√(A_0^2[c_x^2(𝐤)+c_y^2(𝐤)]+R_1^2d_1^2(𝐤)), E_𝐤=√(|M̃_±(𝐤)|^2+D_𝐤^2), Ẽ_𝐤=E_𝐤+D_𝐤, N_𝐤=√(2E_𝐤Ẽ_𝐤), D_±(𝐤)=D_𝐤± R_1d_1(𝐤), D̃_𝐤=√(2D_𝐤D_-(𝐤)). The other pseudospin basis is related to the one listed above by symmetry

    |𝐤,β'⟩=PT|𝐤,α'⟩=1/D̃_𝐤N_𝐤[ M̃_+(𝐤);     Ẽ_𝐤 ][   -D_-(𝐤); A_0c_-(𝐤) ],

In order for the pseudospin basis to have the correct P and T symmetries in the whole BZ, we fix the wave vectors in |𝐤,α'⟩ and |𝐤,β'⟩ to lie on the northern hemisphere (k_z>0). States on the southern hemisphere are obtained by symmetry, namely |-𝐤,α'⟩=P|𝐤,α'⟩ and |-𝐤,β'⟩=P|𝐤,β'⟩=T|𝐤,α'⟩.





The naive choice of the pseudospin basis defined above are not guaranteed to have the correct rotational property of the original model. As a result, they may not be the suitable basis set for studying the symmetry of a specific pairing channel. As has been shown in Ref.<cit.>, a good set of the pseudospin basis can be constructed as a linear combination of |𝐤,α'⟩ and |𝐤,β'⟩ that make the magnetic moment operator expressed under this basis to have the proper transformation property under rotation <cit.>. For the model defined by Eq.(1), the magnetic moment is a linear combination of 𝐬 and σ_1𝐬. Here, following the method in Ref.<cit.>, we choose to focus on the component m_1α=g_1ασ_0+σ_1/2s_α of the magnetic moment <cit.>. α=x,y,z, and g_1x=g_1y=g_1p are assumed. In the space of {|𝐤,α'⟩, |𝐤,β'⟩}, m_1z has the following matrix form

    m_1z(𝐤)/g_1z|W_𝐤|^2=[                 cosθ_𝐤   ie^i(φ_𝐤+2ϕ_𝐤)sinθ_𝐤; -ie^-i(φ_𝐤+2ϕ_𝐤)sinθ_𝐤                -cosθ_𝐤 ],

where the three phase factors, φ_𝐤, ϕ_𝐤, and θ_𝐤, are defined by

    c_+(𝐤)=i√(c_x^2(𝐤)+c_y^2(𝐤))e^-iφ_𝐤=ic(𝐤)e^-iφ_𝐤,


    W_𝐤=Ẽ_𝐤+M̃_+(𝐤)/√(2)N_𝐤=|W_𝐤|e^iϕ_𝐤,


    R_1d_1(𝐤)+iA_0c(𝐤)=D_𝐤e^iθ_𝐤.

m_1z(𝐤) in the basis of  {|𝐤,α'⟩, |𝐤,β'⟩} clearly does not have the desired form of the z-component of an axial vector. The new basis |𝐤,α⟩ and |𝐤,β⟩ are constructed such that m_1z(𝐤) is proportional to the z-component of the Pauli matrix, namely they are the eigenbasis of m_1z <cit.>. We take

    |𝐤,α⟩=h(𝐤)[(1+cosθ_𝐤)|𝐤,α'⟩ -ie^-i(φ_𝐤+2ϕ_𝐤)sinθ_𝐤|𝐤,β'⟩].

From |𝐤,β⟩=PT|𝐤,α⟩, we get the other basis

    |𝐤,β⟩=h^∗(𝐤)[(1+cosθ_𝐤)|𝐤,β'⟩ -ie^i(φ_𝐤+2ϕ_𝐤)sinθ_𝐤|𝐤,α'⟩].

Normalization of the eigenbasis requires

    |h(𝐤)|^2=1/2(1+cosθ_𝐤)=1/4cos^2θ_𝐤/2.

In this basis, we have m_1z(𝐤)=g_1z|W_𝐤|^2ρ_z, where ρ_z is the conventional z-component of the Pauli matrices. The x-component of 𝐦_1 in the new basis has a purely off-diagonal form with the two off-diagonal elements

    [m_1x(𝐤)]_αβ/g_1p|W_𝐤|^2=[m_1x(𝐤)]^∗_βα/g_1p|W_𝐤|^2=
    [2ih^∗(𝐤)cosθ_𝐤/2e^i(φ_𝐤+ϕ_𝐤)]^2.

Taking

    h(𝐤)=i/2cosθ_𝐤/2e^i(φ_𝐤+ϕ_𝐤),

we have

    m_1x(𝐤)=g_1p|W_𝐤|^2ρ_x,


    m_1y(𝐤)=g_1p|W_𝐤|^2ρ_y,

where ρ_x and ρ_y are the conventional x-component and y-component of the Pauli matrices. Therefore, we have shown that the new basis {|𝐤,α⟩, |𝐤,β⟩} defined by Eqs.(8), (9), and (12) can ensure the correct transformation property of the magnetic moment operator, and are thus proper choices in discussing symmetry properties of the system. This basis, employed in the present work, is  shown <cit.> to coincide with the so-called manifestly covariant Bloch basis introduced by Fu <cit.>.







§ PAIRING AND GAP STRUCTURE OF THE BULK QUASIPARTICLE SPECTRUM



In the Nambu basis, ψ^†_𝐤=[ϕ^†_𝐤,(ϕ_-𝐤)^T], and denoting the pairing term generically as Δ(𝐤), the model for a bulk superconducting topological insulator is written as

    Ĥ   =   1/2∑_𝐤ψ^†_𝐤[    H_0(𝐤)-μ I_4            Δ(𝐤);        -Δ^∗(-𝐤) μ I_4-H^∗_0(-𝐤) ]ψ_𝐤
       =   1/2∑_𝐤ψ^†_𝐤H(𝐤)ψ_𝐤,

where μ is the chemical potential. The 1/2 factor accounts for the particle-hole redundancy introduced by the Nambu representation.




The two-fold in-plane rotation symmetry in the Knight shift and field-angle dependent specific heat experiments indicate that the pairing must belong to a multi-dimensional representation of the symmetry group. Because the pairing order parameter for a one-dimensional representation should necessarily respect the three-fold rotational symmetry of the D_3d^5 space group. Presently, most of attention has been paid to the two-dimensional E_u representation of the D_3d^5 space group. One set of the two bases for the E_u representation is Δ_4a(𝐤)=iΔ_aσ_2⊗ s_0 and Δ_4b(𝐤)=Δ_bσ_2⊗ s_3, with Δ_a and Δ_b the pairing amplitudes. For a simplified model without the hexagonal warping terms (i.e., R_1=R_2=0), both of the two components are known to lead to bulk spectrum with point nodes <cit.>. However, it was shown by Fu that the bulk nodes for Δ_4a(𝐤) are gapped out by including the hexagonal warping term proportional to R_1 of Eq.(1) <cit.>. It seems that the fully gapped Δ_4a(𝐤) provides a most natural explanation to the Knight shift and field-angle dependent specific heat experiments.




Here, we study more carefully the excitation gap of the bulk quasiparticle spectrum. Since only the conduction band contribute to the Fermi surface, the gap structure of the quasiparticle spectrum is understood more easily from the low-energy effective model obtained by projecting the full model defined by Eq.(15) to the conduction band <cit.>. The dispersion of the conduction band is ϵ(𝐤)+E_𝐤. The projection is thus achieved by replacing H_0(𝐤)-μ I_4 with (ϵ(𝐤)+E_𝐤-μ)I_2 and transforming the pairing term expressed in the spin-orbital basis to the pseudospin basis derived in the last section. For an arbitrary pairing denoted as Δ(𝐤) in the original spin-orbital basis, its expression in the new pseudospin basis of the conduction band is

    Δ̃(𝐤)=U^†_𝐤ΔU^∗_-𝐤,

where the transformation matrix is U_𝐤=[|𝐤,α⟩,|𝐤,β⟩].



For Δ_4a(𝐤)=iΔ_aσ_2⊗ s_0, we have <cit.>

    Δ̃_4a(𝐤)   =   Δ_0/E_𝐤[-R_1d_1(𝐤)ρ_1 -(B_0c_z(𝐤)+R_2d_2(𝐤))ρ_2
          +A_0c_y(𝐤)ρ_3]iρ_2.

On the Fermi surface, ϵ(𝐤)+E_𝐤-μ=0, the quasiparticle spectrum is determined only by the pairing term

    E(𝐤)=±|det[Δ̃_4a(𝐤)]|
    =±|Δ_a|√(1-M^2(𝐤)+A_0^2c_x^2(𝐤)/E_𝐤^2).

Up to slight hexagonal warping induced by terms proportional to R_1 and R_2, both M(𝐤) and E_𝐤 are approximately symmetrical in the k_xk_y plane. As a result of the c_x(𝐤) term in Eq.(18), the size of superconducting gap is smaller along the k_y=0 contour of the Fermi surface than that along the k_x=0 contour of the Fermi surface. Therefore, the bulk energy spectrum for Δ_4a(𝐤) has a strong anisotropy between the k_x direction and the k_y direction.



Let us focus on the contour of the Fermi surface on the k_y=0 plane, where the minimum of the superconducting gap is attained. Eq.(18) is written as

    E(𝐤)=±|Δ_a|/μ-ϵ_𝐤√(R_1^2d_1^2(𝐤)+B_0^2c^2_z(𝐤)),

where ϵ(𝐤)+E_𝐤-μ=0 has been used. For clarity, we have d_1(𝐤)≃ (k_xa)^3 for k_y=0 and k_x small, and c_z(𝐤)=sin(k_zc). Notice that for both spheroidal and corrugated cylindrical Fermi surfaces, including those shown in Fig.1, d_1(𝐤) and c_z(𝐤) do not attain zero simultaneously. While for R_1=0 there are point nodes on the k_y=0 Fermi surface contour determined by c_z(𝐤)=0, a finite R_1 removes all these nodes <cit.>. One exception is a spheroidal Fermi surface with a point 𝐤=(0,0,π) on it, which marks the transition between a spheroidal Fermi surface and a corrugated cylindrical Fermi surface. For practical purpose, however, we will ignore this special case and so the bulk spectrum of Δ_4a(𝐤) is always fully gapped for R_10.




On the other hand, k_xa on the Fermi surface is small for actual materials. The size of the gap for c_z(𝐤)=0 and k_y=0, which grows like (k_xa)^3 for small k_x, is actually much smaller than Δ_a. In all cases studied, μ-ϵ(𝐤)=E_𝐤 has only a small variation on the Fermi surface. Therefore, the numerator of Eq.(19) determines the qualitative behavior of the superconducting gap. For simplicity and without losing generality, we focus on the k_x≥0 and k_z≥0 portion of the k_y=0 Fermi surface contours shown in Figs.1(a, c, and e). Each point on the chosen portion of the Fermi surface contour can then by labeled by a unique k_x. For the spheroidal Fermi surface shown in Fig.1(a), as we go along the Fermi surface contour from (0,0,k_z1) (k_z1c≃0.29π) to (k_x1,0,0) (k_x1a≃0.077π), c_z(𝐤) decreases monotonously and we have |B_0c_z(0,0,k_z1)|≫|R_1d_1(k_x1,0,0)|. The size of the gap is thus expected to decrease monotonously as we go along the contour from (0,0,k_z1) to (k_x1,0,0). For the two cases with corrugated cylindrical Fermi surfaces, the Fermi surface contour is bounded by two points (k_xi,0,π) and (k_xf,0,0). k_xia≃0.039π (k_xia≃0.074π) and k_xfa≃0.103π (k_xfa≃0.103π) for Fig.1(c) [Fig.1(e)]. As we increase k_x from k_xi to k_xf, k_z changes from π to 0. As a result, we have a nonmonotonous variation of c_z(𝐤), which first increases towards 1 as k_z approaches π/2 and then decreases to 0 afterwards. Because we have |R_1d_1(k_x,0,k_z)/B_0|≪1 along the Fermi surface contour, we expect to get a nonmonotonous variation of the superconducting gap for corrugated cylindrical Fermi surfaces, which first increases and then decreases, with two minima at (k_xi,0,π) and (k_xf,0,0). In comparison to the case for Fig.1(a), the number of gap minima is doubled when the Fermi surface evolves from spheroidal to corrugated cylindrical.



To have a more quantitative understanding on the evolution of the superconducting gap explained above, we plot in Fig.2 simultaneously three functions f_1=4√(R_1^2d_1^2(𝐤)+B_0^2c^2_z(𝐤)), f_2=4(μ-ϵ_𝐤)=4E_𝐤, and f_3=f_1/f_2=|E(𝐤)/Δ_a|, on the k_z≥0 and k_x≥0 portion of the k_y=0 Fermi surface contours. The states are labeled uniquely in terms of the value of k_x. The value of f_3 gives the magnitude of the bulk gap, normalized by the pairing amplitude. For all three sets of parameters considered, the minimal values of the bulk gap are much smaller than the corresponding maximum values. For the experimental transition temperature of about 3.8 Kelvin, the pairing amplitude (i.e., Δ_a) is of the order 1 meV. The minimal value of the bulk gap is two to three orders of magnitude smaller than the pairing amplitude. It is also interesting to notice that, when the smaller Fermi momentum of the corrugated cylindrical Fermi surface along k_x [e.g., k_xia≃0.039π for Fig.2(b)] is smaller than the Fermi momentum of the spheroidal Fermi surface along k_x [e.g., k_x1a≃0.077π for Fig.2(a)], the minimum superconducting gap for the corrugated cylindrical Fermi surface can be smaller than the minimum superconducting gap for the spheroidal Fermi surface. From a practical point of view, and for both spheroidal and corrugated cylindrical Fermi surfaces, the minimum of the superconducting gap acts effectively as point node of the bulk spectrum. The above picture holds as long as R_1 is not extremely (e.g., two to three orders of magnitude) larger than the value used. For the parameter set of `Cylindrical 2' in table I, a two orders of magnitude larger R_1 (i.e., 20 eV) is needed to increase the minimum of the bulk gap along k_x to the same order of magnitude to that along k_y. A further tenfold enhancement in R_1 is required to achieve the same increase in the bulk gap for the parameter set of `Cylindrical 1' in table I. These large values of R_1 are not only inconsistent with the magnitudes of other parameters but also will distort strongly the Fermi surface and thus deviate qualitatively from experiments. Therefore, the bulk spectrum for Δ_4a(𝐤) should be nodal-like for realistic parameters.












For Δ_4b(𝐤)=Δ_bσ_2⊗ s_3, the effective pairing is

    Δ̃_4b(𝐤)   =   Δ_b/E_𝐤[(B_0c_z(𝐤)+R_2d_2(𝐤))ρ_1
    -R_1d_1(𝐤)ρ_2
          -A_0c_x(𝐤)ρ_3]iρ_2.

The minimum of the superconducting gap lies along the k_x=0 plane, where d_1(𝐤)=c_x(𝐤)=0. Along the intersection contour of the Fermi surface with the k_x=0 plane, E_𝐤=μ-ϵ(𝐤) is a smooth function of the wave vector. The variation of the gap is thus determined by B_0c_z(𝐤)+R_2d_2(𝐤). For R_2=0, we reproduce the known result that Δ_4b(𝐤) has bulk point nodes determined by c_z(𝐤)=k_x=0. The number of the point nodes is two (four) for spheroidal (corrugated cylindrical) Fermi surface. For R_20, the above point nodes are gapped out, with the magnitude of the gap proportional to |Δ_bR_2d_2(0,k_Fy,0)|/E_𝐤 [and also |Δ_bR_2d_2(0,k'_Fy,π)|/E_𝐤 for corrugated cylindrical Fermi surface], where k_Fy (k'_Fy) is the k_y component of the Fermi momentum. The original point nodes do not simply vanish. Instead, they are tilted away from the (0,k_y,0) axis [and also the (0,k_y,π) axis for the case with corrugated cylindrical Fermi surface] into the k_yk_z plane. If we have |R_2/B_0|≫1, a fully-gapped bulk spectrum can be obtained. However, for realistic parameters, the point nodes are still present. For the parameters considered in Table I, the point nodes of Δ_4b(𝐤) are in fact still very close to the point nodes for R_2=0.






§ SURFACE ANDREEV BOUND STATES



Since Δ_4b(𝐤) is nodal for practical model parameters, we will focus on the fully-gapped Δ_4a in what follows. For spheroidal Fermi surface, Δ_4a was known to support a peculiar surface Andreev bound states (SABSs) on the xy surface of a sample, which is (almost) flat along the k_x direction of the surface BZ <cit.>. It was argued in later works that, when the Fermi surface becomes two-dimensional-like with copper intercalation, the SABSs for Δ_4a would disappear <cit.>. This conclusion is natural if the Fermi surface is purely cylindrical with no dispersion along k_z, because the existence of SABSs on the xy surface is associated with a sign change in the pairing term upon reflection from the surface, which requires on one hand a finite dispersion along k_z and on the other hand a pairing component that changes sign with the reversal of k_z.



However, there seems to be no reason why the Fermi surface can turn from three-dimensional to completely two-dimensional with copper intercalation. On one hand, one experiment reports corrugated cylindrical rather than completely cylindrical Fermi surface <cit.>. On the other hand, Cu_xBi_2Se_3 is known to be superconducting in a wide range of x values <cit.>. It is natural to expect that the Fermi surface evolves continuously from spheroidal to corrugated cylindrical as x increases. Finally, a completely two-dimensional Fermi surface is inconsistent with the existence of the TSSs, observed experimentally <cit.>. Therefore, compared to completely cylindrical Fermi surface, corrugated cylindrical Fermi surfaces with different degrees of corrugation [e.g., those shown in Fig.1] are better descriptions of the actual Fermi surface. For these Fermi surfaces, the bulk gap for the Δ_4a pairing has minima along the k_x direction, which is explained above and illustrated in Fig.2. Since the conduction band has a finite dispersion along k_z, and the gap is dominated by the c_z(𝐤) term which is odd in k_z, it seems natural to expect the prevalent existence of SABSs. As shown below, this is true. We first give some analytical analysis, which is then followed by numerical results. Finally, we study the stability of the SABSs against surface and bulk nonmagnetic impurities.






 §.§ Analytical analysis for a simplified model



To gain a qualitative understanding of the SABSs for a Fermi surface in the shape of a spheroid or a corrugated cylinder, we ignore the hexagonal warping in the Fermi surface and consider a band with the following dispersion

    ξ_𝐤=k_x^2+k_y^2/2m_1^∗-t_zcos k_z-μ,

where m_1^∗ is the effective mass in the k_xk_y plane, and the dispersion along k_z is determined by t_z. ħ=1 is assumed. m_1^∗ is assumed to be small so that the Fermi momenta along directions in the k_xk_y plane are small, i.e. √(2m_1^∗μ) is small (m_1^∗>0 and μ>0). The Fermi surface is spheroidal (corrugated cylindrical) when |t_z/μ|≥1 (0<|t_z/μ|<1).




For the pairing term, it is convenient to replace the factor E_𝐤 with a constant. Eq.(17) is thus reduced to

    Δ̃_4a(𝐤)   =   Δ̃_0[-R_1d_1(𝐤)ρ_1 -(B_0c_z(𝐤)+R_2d_2(𝐤))ρ_2
          +A_0c_y(𝐤)ρ_3]iρ_2,

where Δ̃_0 is the pairing amplitude divided by the constant representing E_𝐤.




There are several available approaches that we can use to derive the SABSs. Here, we follow the approach of mapping the surface problem by an equivalent junction problem. Namely, we consider a junction at z=0 between the surfaces of two bulk samples; one is extended from z=0 to ∞, and the other  is from z=0 to -∞. The two bulk samples are both described by Eq.(21). The problem of scattering off the surface is mapped to a sign change in the components of the pairing term odd in k_z. The z<0 part of the junction is described simply by Eqs.(21) and (22). The z>0 part of the junction is described by Eq.(21) and Eq.(22) with the sign of the term proportional to c_z(𝐤) reversed. The model on either side of the junction is a 4×4 model with two 2×2 block diagonals for the bare bands and two 2×2 off-diagonal blocks representing the pairing term.



To proceed, we adopt the quasiclassical approximation to the ansatz for the wave function of the SABSs, to separate the fast and slow degrees of freedom. The 4×4 eigenvector of the SABSs is thus taken as

    φ(𝐤_F,𝐫)=[ u(𝐫); v(𝐫) ]=
    e^i𝐤_F·𝐫[ f(𝐫); g(𝐫) ].

In the same spirit, we expand the bulk band around the Fermi momentum 𝐤_F as

    ξ_𝐤≃𝐯_F(𝐤_F)·(-i∇-𝐤_F),

where the Fermi velocity is defined as

    𝐯_F(𝐤_F)=∇_𝐤ξ_𝐤|_𝐤=𝐤_F.

We assume a perfect junction in which the translational invariance within the junction plane is preserved. The problem is thus reduced to a one-dimensional scattering problem along the z direction. That is, the dependencies in the x and y coordinates occur only through the exponential pre-factor of Eq.(23). f and g depend only on z. Consistent with this assumption on the wave function, we replace in the pairing term the k_x and k_y components of the wave vectors with k_Fx and k_Fy. The c_z(𝐤) term is then expanded to linear order of k_z-k_Fz.



The interface localized states are solved by imposing the following boundary conditions to the wave function

    φ(𝐤_F,x,y,z=0^+)=φ(𝐤_F,x,y,z=0^-),


    φ(𝐤_F,x,y,z=-∞)=φ(𝐤_F,x,y,z=+∞)=0.

0^+ and 0^- are positive and negative infinitesimals. Focusing on the direction of k_Fy=0, where the minima of the superconducting gap are attained, we indeed find solutions satisfying the above boundary conditions, with energies

    |E(k_Fx,0,k_Fz)|=|Δ̃_0R_1d_1(k_Fx,0,k_Fz)|.

From the discussions on the bulk superconducting gap in the previous section, the energy of the above bound states are well below the bulk gap in a large part of the bulk gap. Therefore, they are well-defined in-gap states.






 §.§ Numerical results for clean system




Inspired by the existing experiments, we consider the three sets of parameters in Table I, which result separately a spheroidal Fermi surface and two corrugated cylindrical Fermi surfaces with different degrees of corrugation. Besides the shape of the Fermi surface, these parameters allow the simultaneous presence of bulk conduction band and the topological surface states (TSSs) at the Fermi level. As shown in Figure 3 are the surface spectral functions for the three parameter sets, for clean systems. The surface spectral function for a wave vector in the surface BZ are defined as summation over the imaginary part of the particle Green's function, which are obtained in terms of standard iterative Green's function method <cit.>. From Fig.3, the SABSs exist for both spheroidal and corrugated cylindrical Fermi surfaces. One essential feature is the existence of a nearly flat band of Andreev bound states along k_x at the center of the SC gap. The magnitude of the bulk gap depends sensitively on the model parameters and can be vanishingly small for the parameters with (corrugated) cylindrical Fermi surface, consistent with the previous section. A qualitative difference from previous results obtained for simplified models without the hexagonal warping terms is that <cit.>, the SABSs are not exactly flat along (k_x,0), which becomes increasingly clear as the size of the (larger) gap minimum increases from Fig.1(a) to Fig.1(c) and Fig.1(e).



The surface spectral functions can be probed by ARPES <cit.>. The integrated surface spectral function, the surface density of states (SDOS), can be probed by tunneling spectroscopy <cit.>. In Fig.4, we have shown the SDOS together with the corresponding bulk density of states (BDOS). A common characteristic of the results for all three parameter sets is the appearance of a prominent zero-energy peak corresponding to the mid-gap Majorana bound states plus a continuum of low-energy states filling up the bulk superconducting gap. As is explained in Sec.IV, the minimum of the bulk gap scales linearly with R_1. We have made test calculations by increasing R_1 artificially to 2 eV and 20 eV, for the third set of parameters (`Cylindrical 2' in Table I). The minimal size of the bulk gap along k_x increases linearly with R_1, and the BDOS becomes increasingly U shaped with a flat bottom of zero DOS. The SABSs, on the other hand, persist and traverse the bulk gap for all parameters considered. The zero-energy TSSs also persist. As a result, the SDOS is still featured by the existence of in-gap states with a peak at or close to zero energy.
















 §.§ Stability against nonmagnetic impurities



Having verified the presence of SABSs and seen their peculiar dispersions in clean systems, we proceed to test their stability against various imperfections. From a practical point of view relevant to Cu_xBi_2Se_3, it is plausible to focus on the effect of short-range nonmagnetic impurities. We study the effect of the nonmagnetic impurities at two levels, impurities uniformly distributed in the whole sample and separate point-like impurities situating on the surface.



Firstly, we consider the effect of impurities distributed uniformly throughout the whole sample. We first obtain the self-energy correction to the bulk Green's functions in terms of the self-consistent T-matrix approximation <cit.>. We consider the simplest case of uniformly distributed short-range nonmagnetic impurities, V(𝐫)=V_0δ(𝐫-𝐫_0). In this case, the self-energy is 𝐤-independent and is determined by a set of three self-consistent equations: (1) G(𝐤,ω)=[ω+iη-H(𝐤)-Σ(ω)]^-1, (2) Σ(ω)=n_imp[T(ω)-Ṽ], and (3) T(ω)=[I_8-Ṽ/N∑_𝐤G(𝐤,ω)]^-1Ṽ. Here, n_imp is the concentration of the nonmagnetic impurity, Ṽ=V_0τ_3⊗σ_0⊗ s_0, 𝐤 denotes a wave vector in the 3D BZ, and N is the number of wave vectors in the 3D BZ. After obtaining the self-energy Σ(ω) from the above self-consistency loop, we add it as an energy correction to H(𝐤) and obtain the surface Green's function in terms of the iterative Green's function method <cit.>. The resulting surface Green's function is then the proper Green's function for the surface layer in the presence of short-range nonmagnetic impurities uniformly distributed in the bulk. The bulk density of states (BDOS) and surface density of states (SDOS) obtained by this method are as shown in Figure 5. Three concentrations of the impurities (n_imp=0.001, 0.01 and 0.02) are considered. As the bulk superconducting gap of the cases with corrugated cylindrical (spheroidal) Fermi surface is filled up for n_imp=0.01 (n_imp=0.02), the fine structures in the SDOS beyond E=0 disappear, but the zero-energy surface states are still quite robust and manifest as a single zero-energy peak in the SDOS.










Secondly, we consider the effect of individual point-like impurities on the surface of an otherwise clean sample. This is achieved by keeping the bulk of the material clean, and adding impurities only to the surface layer in a manner that different impurities are far away from each other. For this case, we study the changes in the surface Green's functions for the clean system induced by the surface impurities. The effect of the impurities is taken in to account in terms of the T-matrix approximation <cit.>. For a single short-range nonmagnetic impurity, V(𝐫)=V_0δ(𝐫-𝐫_0), the T-matrix is 𝐤-independent

    T(ω)=[I_8-Ṽ/N_xy∑_𝐤G_0(𝐤,ω)]^-1Ṽ,

where N_xy is the number of unit cells (wave vectors) in the xy plane (surface BZ), G_0(𝐤,ω) is the retarded surface Green's function obtained in terms of the iterative Green's function method for a clean system <cit.>. Ṽ=V_0τ_3⊗σ_0⊗ s_0 is the impurity potential in the Nambu space. In terms of the unperturbed Green's function G_0(𝐤,ω) and the T-matrix, the perturbed Green's function is obtained in term of the T-matrix approximation as

    G(𝐫,𝐫',ω)=G_0(𝐫-𝐫',ω)+G_0(𝐫-𝐫_0,ω)T(ω)G_0(𝐫_0-𝐫',ω),

where G_0(𝐫,ω) is the Fourier transformation of G_0(𝐤,ω). The local Green's function at 𝐫 under the influence of the nonmagnetic impurity at 𝐫_0 is defined as G(𝐫,𝐫,ω). We consider the effect of a strong unitary impurity and take V_0=1000 eV. The results of the local density of states (SLDOS) at the impurity site and its six nearest-neighbor sites for the three sets of parameters in Table I are shown in Fig. 6. The SLDOS on the nearest-neighbor site 𝐫-𝐫_0=δ_α is equal to the SLDOS on the nearest-neighbor site 𝐫-𝐫_0=-δ_α (α=1,2,3), so only results for three of the six nearest-neighbor sites are shown. Albeit quantitative differences from the SDOS on Fig.4 and Fig.5, the in-gap surface states persist. On the other hand, the SLDOS for 𝐫-𝐫_0=δ_1 and 𝐫-𝐫_0=δ_2 are identical but are different from the SLDOS for 𝐫-𝐫_0=δ_3, which is consistent with the two-fold anisotropy of the superconducting pairing between the x and y directions. It is also noted that the zero-energy surface states are also very robust at and near the impurity site for the three types of Fermi surfaces (see Fig. 6).











§ CONCLUSION



In conclusion, extensive analyses for the bulk and surface spectra have been made for the Δ_4a(𝐤) nematic pairing, suggested to be the correct SC pairing for the Cu_xBi_2Se_3 compound based on recent Knight shift and field-angle dependent specific heat measurements<cit.>. The purpose of the present work is to explore the consequences deduced from this type of pairing, taking into account the evolution of the Fermi surface from spheroidal to corrugated cylindrical <cit.>. We show that Cu_xBi_2Se_3 with Δ_4a(𝐤) pairing should be a topological superconductor with topological surface states and the surface Andreev bound states (SABSs), even if the Fermi surface has changed from spheroidal <cit.> to the corrugated cylindrical case studied in the present work. The bulk SC spectrum, while fully gapped, show prominent twofold anisotropy with vanishingly small gap minima along k_x. One of the essential features of the SABSs is the exhibition of the zero-energy Majorana bound states regardless the shape of the Fermi surface. The SABSs are shown to be robust against short-range bulk and surface nonmagnetic impurities. This is consistent with recent works which show that the surface states of class DIII topological superconductor are stable against weak disorder and interaction <cit.>. On the other hand, there are experimental controversies on the existence of SABSs <cit.>. As to whether the SABSs could be observed experimentally may depend critically on the condition of the sample surfaces. For instance, the excessive magnetic Cu^2+ (3d^9) ions or Cu (3d^104S^1) atoms on the surfaces could very much suppress the SABSs. For clean and perfect surface, on the other hand, the SABSs should be detectable. The present pairing model, an odd-parity spin-triplet pairing, yields a very anisotropic bulk SC gap while the STM experiment detected an almost isotropic s-wave like bulk gap<cit.>. It appears that the present nematic pairing model, which has successfully explained recent Knight shift and field-angle dependent specific heat measurements, are apparently having difficulty to account for other experimental measurements like the SC density of states <cit.>. Therefore it is necessary to develop a revised pairing model which includes the essential physics as discussed here and also being able to explain other experiments. This will definitely constitute a challenging topic for future study.




L.H. thanks Ting-Kuo Lee and Sungkit Yip for many helpful discussions. This work was supported in part by the Robert A. Welch Foundation under Grant No. E-1146 and the Texas Center for Superconductivity at the University of Houston. L.H. is also supported by NSFC.11204035 and would also like to acknowledge the support from the China Scholarship Council. Part of the calculations were performed in the Center for Advanced Computing and Data Systems in University of Houston.







hor10 Y. S. Hor, A. J. Williams, J. G. Checkelsky, P. Roushan,
J. Seo, Q. Xu, H. W. Zandbergen, A. Yazdani, N. P. Ong,
and R. J. Cava, Phys. Rev. Lett. 104, 057001 (2010).


wray10 L. Andrew Wray, Su-Yang Xu, Yuqi Xia, Yew San Hor,
Dong Qian, Alexei V. Fedorov, Hsin Lin, Arun Bansil,
Robert J. Cava and M. Zahid Hasan, Nature Phys. 1762
(2010).


sasaki11 Satoshi Sasaki, M. Kriener, Kouji Segawa, Keiji Yada,
Yukio Tanaka, Masatoshi Sato, and Yoichi Ando, Phys.
Rev. Lett. 107, 217001 (2011).

koren11 G. Koren, T. Kirzhner, E. Lahoud, K. B. Chashka, and
A. Kanigel, Phys. Rev. B 84, 224521 (2011).

kirzhner12 T. Kirzhner, E. Lahoud, K. B. Chaska, Z. Salman, and
A. Kanigel, Phys. Rev. B 86, 064517 (2012).

kriener11 M. Kriener, Kouji Segawa, Zhi Ren, Satoshi Sasaki, and
Yoichi Ando, Phys. Rev. Lett. 106, 127004 (2011).

das11 Pradip Das, Yusuke Suzuki, Masashi Tachiki, and Kazuo
Kadowaki, Phys. Rev. B 83, 220513(R) (2011).

bay12 T.V. Bay, T. Naka, Y. K. Huang, H. Luigjes, M. S.
Golden, and A. de Visser, Phys. Rev. Lett. 108, 057001
(2012).


levy12 N. Levy, T. Zhang, J. Ha, F. Sharifi, A. A. Talin, Y. Kuk, and
J. A. Stroscio, Phys. Rev. Lett. 110, 117001 (2013).


peng13 Haibing Peng, Debtanu De, Bing Lv, Fengyan Wei, and Ching-Wu Chu, Phys. Rev. B 88, 024515 (2013).


matano16 K. Matano, M. Kriener, K. Segawa, Y. Ando, and Guo-qing Zheng, Nat. Phys. 12, 852 (2016).

fu16 Liang Fu, Nat. Phys. 12, 822 (2016).


zhao14 Lukas Zhao, Haiming Deng, Inna Korzhovska,	Milan Begliarbekov,	Zhiyi Chen,	Erick Andrade, Ethan Rosenthal,	Abhay Pasupathy, Vadim Oganesyan, Lia Krusin-Elbaum, Nature Commun. 6, 8279 (2014).


pan16 Y. Pan, A. M. Nikitin, G. K. Araizi, Y. K. Huang, Y. Matsushita, T. Naka, and A. de Visser, Sci. Rep. 6, 28632 (2016).


du17 Guan Du, YuFeng Li, J. Schneeloch, R. D. Zhong, GenDa Gu, Huan Yang, Hai Lin, and Hai-Hu Wen, Sci. China-Phys. Mech. Astron. 60, 037411 (2017)


qiu15 Y. Qiu, K. N. Sanders, J. Dai, J. E. Medvedeva, W. Wu, P. Ghaemi, T. Vojta, and Y. S. Hor, arXiv:1512.03519.


asaba17 Tomoya Asaba, B. J. Lawson, Colin Tinsman, Lu Chen, Paul Corbae, Gang Li, Y. Qiu, Y. S. Hor, Liang Fu, and Lu Li, Phys. Rev. X 7, 011009 (2017).


lawson16 B. J. Lawson, Paul Corbae, Gang Li, Fan Yu, Tomoya Asaba, Colin Tinsman, Y. Qiu, J. E. Medvedeva, Y. S. Hor, and Lu Li, Phys. Rev. B 94, 041114(R) (2016).



wang16 Zhiwei Wang, A. A. Taskin, Tobias Frölich, Markus Braden, and Yoichi Ando, Chem. Mater. 28, 779 (2016).




yonezawa17 Shingo Yonezawa, Kengo Tajiri, Suguru Nakata, Yuki Nagai, ZhiweiWang, Kouji Segawa,
Yoichi Ando, and Yoshiteru Maeno, Nature Phys. 13, 123 (2017).


fu10 Liang Fu and Erez Berg, Phys. Rev. Lett. 105, 097001
(2010).

fu14 Liang Fu, Phys. Rev. B 90, 100509 (2014).



yip13 S.-K. Yip, Phys. Rev. B 87, 104505 (2013).


hao14 Lei Hao and Ting-Kuo Lee, J. Phys.: Condens. Matter 27, 105701 (2015); ibid arXiv:1407.3329v2.



hao11 Lei Hao and T. K. Lee, Phys. Rev. B 83, 134516 (2011).


yamakage12 Ai Yamakage, Keiji Yada, Masatoshi Sato, and Yukio Tanaka, Phys. Rev. B 85, 180509(R) (2012).


hao15 Lei Hao and Jun Wang, J. Phys.: Condens. Matter 27, 255701 (2015).



lahoud13 E. Lahoud, E. Maniv, M. S. Petrushevsky, M. Naamneh, A. Ribak, S. Wiedmann, L. Petaccia, Z. Salman, K. B. Chashka, Y. Dagan, and A. Kanigel, Phys. Rev. B 88, 195107 (2013).


mizushima14 Takeshi Mizushima, Ai Yamakage, Masatoshi Sato, and Yukio Tanaka, Phys. Rev. B 90, 184516 (2014).



liu15 Zhongheng Liu, Xiong Yao, Jifeng Shao, Ming Zuo, Li Pi, Shun Tan, Changjin Zhang,
and Yuheng Zhang, J. Am. Chem. Soc. 137, 10512 (2015).


nagai15 Yuki Nagai, Phys. Rev. B 91, 060502(R) (2015).



liu10 Chao-Xing Liu, Xiao-Liang Qi, Haijun Zhang,
    Xi Dai, Zhong Fang, and Shou-Cheng Zhang, Phys. Rev. B 82,
    045122 (2010).


fu09 Liang Fu, Phys. Rev. Lett. 103, 266801
    (2009).


zhang09 Haijun Zhang, Chao-Xing Liu, Xiao-Liang Qi,
    Xi Dai, Zhong Fang, and Shou-Cheng Zhang, Nature Phys.
5, 438 (2009).


wang10 Qiang-Hua Wang, Da Wang, and Fu-Chun Zhang,
    Phys. Rev. B 81, 035104 (2010).


fu07 Liang Fu and C. L. Kane, Phys. Rev. B 76, 045302 (2007).



hao13 Lei Hao, Gui-Ling Wang, Ting-Kuo Lee, Jun Wang, Wei-Feng Tsai, and Yong-Hong Yang, Phys. Rev. B 89, 214505 (2014).



acparameters  P. Larson, V. A. Greanya, W. C. Tonjes, Rong Liu, S. D. Mahanti, C. G. Olson,
     Phys. Rev. B 65, 085108 (2002).



yip14 Sungkit Yip, Annu. Rev. Cond. Matter Physics 5, 15 (2014).


yip16 Sungkit Yip, arXiv:1609.04152.


zocher13 Björn Zocher and Bernd Rosenow, Phys. Rev. B 87, 155138 (2013).


hashimoto13 Tatsuki Hashimoto, Keiji Yada, Ai Yamakage, Masatoshi Sato, and Yukio Tanaka, J. Phys. Soc. Jpn. 82, 044704 (2013).


nagai14 Yuki Nagai, Hiroki Nakamura, and Masahiko Machida, J. Phys. Soc. Jpn. 83, 053705 (2014).


takami14 Shota Takami, Keiji Yada, Ai Yamakage, Masatoshi Sato, and Yukio Tanaka, J. Phys. Soc. Jpn. 83, 064705 (2014).



fu15 Liang Fu, Phys. Rev. Lett. 115, 026401 (2015).


kozii15 Vladyslav Kozii and Liang Fu, Phys. Rev. Lett. 115, 207002 (2015).


venderbos16a Jörn W. F. Venderbos, Vladyslav Kozii, and Liang Fu, Phys. Rev. B 94, 094522 (2016).


venderbos16b Jörn W. F. Venderbos, Vladyslav Kozii, and Liang Fu, Phys. Rev. B 94, 180504(R) (2016).



note1 Qualitatively the same effective pairing can be obtained in terms of a simplified method of deriving the low-energy effective model <cit.>. The only change is to substitute μ+M(𝐤)-ϵ(𝐤) for E_𝐤 in the denominator of Eq.(17). The simplified method thus gives the same picture about the nature of the pairing, at the expense of less quantitative accuracy <cit.>.



balatsky06 A. V. Balatsky, I. Vekhter, and Jian-Xin Zhu, Rev. Mod. Phys. 78, 373 (2006).



foster14 Matthew S. Foster, Hong-Yi Xie, and Yang-Zhi Chou, Phys. Rev. B 89, 155140 (2014).


xie15 Hong-Yi Xie, Yang-Zhi Chou, and Matthew S. Foster, Phys. Rev. B 91, 024203 (2015).




