Possible two-component pairings in electron-doped Bi 2 Se 3 based on a tight-binding model
aa r X i v : . [ c ond - m a t . s up r- c on ] D ec Possible two-component pairings in electron-doped Bi Se based on a tight-bindingmodel Lei Hao and C. S. Ting School of Physics, Southeast University, Nanjing 211189, China Department of Physics and Texas Center for Superconductivity,University of Houston, Houston, Texas 77204, USA (Dated: December 20, 2018)Recent experiments show the spontaneous breaking of rotational symmetry in the superconductingtopological insulators M x Bi Se (M represents Cu, Sr, or Nd), suggesting that the pairing belongsto a two-dimensional representation of the D d symmetry group of the crystal. Motivated by theseprogresses, we construct an exhaustive list of possible two-component pairings of the M x Bi Se superconductors, both for the odd-parity E u representation and for the even-parity E g representa-tion. Starting from a tight-binding model for the normal phase of Bi Se and M x Bi Se , we firstlyconstruct the pairing channels in the spin-orbital basis, up to second-nearest-neighbor pairing cor-relations in the basal plane. We then infer the properties of these pairings by transforming themto the band (pseudospin) basis for the conduction band. A comparison with the key experimentalconsensuses on M x Bi Se superconductors shows that the true pairings should also be multichan-nel. Besides a nematic and time-reversal symmetric pairing combination, the other pairings thatwe have identified are chiral and nematic at the same time, which may be nonunitary and have aspontaneous magnetization. A complementary set of experiments are proposed to identify the truepairing symmetries of these superconductors and their evolution with the doping concentration x . I. INTRODUCTION
The nature of the pairing of superconducting topolog-ical insulators, to be abbreviated as M x Bi Se with Mrepresenting a metallic element that might be Cu, Sror Nd, has been mysterious since their discovery [1–6].Various experiments made conflicting implications, al-luding the pairing to be topologically nontrivial [7–11]or trivial [12, 13]. Recently, a series of new experimentshave shown convincingly that the pairing of this fam-ily of superconductors is unconventional (see Yonezawa[14] for a recent review). For Cu x Bi Se [15, 16] andSr x Bi Se [17, 18], more than one experiments have re-vealed that the superconducting state breaks the three-fold rotational symmetry of the normal phase to two-foldrotational symmetry, which is possible only if the pair-ing belongs to a two-dimensional representation of theunderlying D d point group. The pairing was therebycalled nematic superconductivity [19]. For Nd x Bi Se ,in addition to the broken three-fold in-plane rotationalsymmetry [20, 21], the time-reversal symmetry appearsto be also broken [6].Multiple theoretical analysis have been made to ac-count for these new experimental findings [22–31], whichfocus mostly on the odd-parity E u representation of the D d point group. While the studies can account for thequalitative features of various experiments, a satisfac-tory explanation of all crucial experimental features interms of known pairings appears to be difficult. For ex-ample, while recent experiments deny the presence of in-gap states on the surface of superconducting Cu x Bi Se [12, 13], the most studied E u pairing is shown to have ro-bust low-energy surface states [31]. In addition, the pro-posed E u pairing has not been stabilized as the leadingpairing instability in any theoretical calculations based on a microscopic pairing mechanism, such as the electron-phonon interaction [32–35] or the electron-electron inter-action [36]. In view of the difficulty in first-principlespredictions of the pairing symmetry on one hand andthe extensive experimental observations accumulated upto date on the other hand, a promising approach is toconstruct an exhaustive list of possible two-componentpairings and compare them with the available experi-mental consensuses. From this comparison, we may seeto what extent the existing experiments have constrainedthe pairing symmetry, and what further experiments arenecessary to figure out definitely the genuine pairing sym-metries of the M x Bi Se superconductors.Motivated by the above considerations, we constructin this work a complete list of pairing channels belongingto the two-dimensional irreducible representations of the D d point group. Starting from a tight-binding modelfor Bi Se and the normal phase of M x Bi Se , we con-struct pairings belonging to the odd-parity E u represen-tation and pairings in the even-parity E g representation.In consistency with the tight-binding model which is upto second-nearest-neighbor (2NN) in-plane hoppings, theconstructed pairing channels are restricted to 2NN in-plane pairing correlations. By transforming from thespin-orbital basis to the band (pseudospin) basis and re-taining only the conduction band which contributes tothe Fermi surface, we analyze the general properties ofvarious interesting and typical pairing channels.After a comprehensive review over the experimentalconsensuses on these superconductors, including theirbulk spectrum, surface spectrum, and magnetic prop-erties, we infer the constraints of these experiments onthe pairing symmetry. The major conclusion of the com-parison is that the pairing has to be multichannel, inaddition to having two components. A purely nematicpairing combination in the E u representation can givea fully-gapped and two-fold symmetric bulk spectrum,fully-gapped surface spectrum, and two-fold symmetricelectronic spin susceptibility. In addition, we find sev-eral chiral and nematic pairing combinations that canexplain the key experimental results, in both the E g rep-resentation and the E u representation. The chiral andnematic pairing of the E u representation, besides break-ing the time-reversal symmetry, may also be nonunitaryand having a spontaneous magnetization. We finally dis-cuss the implications of the present work to the futureexperiments, which are highly desirable to determine thetrue pairing symmetry of the M x Bi Se superconductors. II. MODEL IN THE SPIN-ORBITAL ANDPSEUDOSPIN BASIS
The low-energy band structures of Bi Se andM x Bi Se (M denotes Cu, Sr, or Nd) can be describedby the following two-orbital tight-binding model, definedon a quasi-two-dimensional hexagonal lattice [31, 36–39], H ( k ) = ǫ ( k ) I + M ( k )Γ + B c z ( k )Γ + A [ c y ( k )Γ − c x ( k )Γ ] + R d ( k )Γ + R d ( k )Γ . (1)The basis operator is taken as φ † k = [ a † k ↑ , a † k ↓ , b † k ↑ , b † k ↓ ],where the a and b orbitals separately correspond to the p z orbitals of the top and bottom Se layers of the Bi Se quintuple units, with a certain amount of hybridizationwith the p z orbitals of the neighboring Bi layers [37–39]. I is the 4 × = σ ⊗ s , Γ = σ ⊗ s ,Γ = σ ⊗ s , Γ = − σ ⊗ s , and Γ = σ ⊗ s [32, 36–41]. s i and σ i ( i = 1 , ,
3) are Pauli matrices for thespin and orbital degrees of freedom, s and σ are thecorresponding unit matrices. With the parity operator P = σ ⊗ s , it is easy to verify that the model has theinversion symmetry P H ( k ) P − = H ( − k ).The above model was obtained previously [36] basedon symmetry analysis and comparison with a k · p modeldefined near k x = k y = k z = 0 [37]. The lattice of Bi Se and M x Bi Se , which belong to the D d space group, ismapped 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 ˚A and 3 c =28.64 ˚A [42]. ǫ ( k ) = C + 2 C [1 − cos( k · δ )]+ C [3 − cos( k · δ ) − cos( k · δ ) − cos( k · δ )]. M ( k ) is obtained from ǫ ( k ) by making thesubstitutions C i → M i ( i = 0 , , c x ( k ) = √ [sin( k · δ ) − sin( k · δ )], c y ( k ) = [sin( k · δ )+sin( k · δ ) − k · δ )], and c z ( k ) = sin( k · δ ). Finally, d ( k ) = − √ [sin( k · a ) + sin( k · a ) + sin( k · a )] and d ( k ) = − k · δ ) +sin( k · δ )+sin( k · δ )]. Here, the four NN bond vectors ofthe hexagonal lattice are δ = ( √ a, a, δ = ( − √ a , a, δ = (0 , − a, δ = (0 , , c ). The threein-plane 2NN bond vectors in d ( k ) are a = δ − δ , a = δ − δ , and a = δ − δ . The last and second last terms of H ( k ) induce hexagonal warping of the Fermisurface and the topological surface states [37, 38]. Wemention in passing that Nd x Bi Se was reported to havemultiple Fermi surfaces, with possible contributions fromthe d orbitals of the Nd dopants [43]. We will neglect thiscomplexity and work with the above model for all threesuperconductors [26, 27].The dopants of M x Bi Se (M is Cu, Sr, or Nd) dopeelectrons to the system, so that only the conduction bandcontributes to the Fermi surface. In studying the bulkproperties of a superconductor, it is easier to work withthe band (pseudospin) basis and retain only the states ofthe conduction band [44–50]. For this purpose, we turn tothe pseudospin basis for the conduction band. BecauseM x Bi Se in the normal state has both inversion sym-metry and time-reversal symmetry, the conduction bandis twofold degenerate (the Kramers degeneracy) at eachwave vector k . The operator for the inversion symmetryis the parity operator P = σ ⊗ s . The time reversal op-erator is taken as T = − iσ ⊗ s K , where K denotes thecomplex conjugation. We define the pseudospin basis forthe conduction band states on the northern hemisphere(i.e., k z >
0) of the three-dimensional Brillouin zone (BZ)as [ | k , α i , | k , β i ] = [ | k , α ′ i , | k , β ′ i ] u k , (2)where α and β are the two pseudospin degrees of free-dom. The Kramers degeneracy relates the two bases via | k , β i = P T | k , α i . The two auxiliary bases are definedas [31] | k , α ′ i = 1˜ D k N k (cid:18) ˜ E k ˜ M − ( k ) (cid:19) (cid:18) A c + ( k ) D − ( k ) (cid:19) , (3)and | k , β ′ i = P T | k , α ′ i = 1˜ D k N k (cid:18) ˜ M + ( k )˜ E k (cid:19) (cid:18) − D − ( k ) A c − ( k ) (cid:19) , (4)where the first and second two-component vectors areseparately spinors in the subspaces of the original orbitaland spin degrees of freedom. The unitary matrix con-necting the two basis sets is [31] u k = (cid:18) ie i ( ϕ k + φ k ) cos θ k − e iφ k sin θ k e − iφ k sin θ k − ie − i ( ϕ k + φ k ) cos θ k (cid:19) . (5)For notational simplicity, we have introduced the fol-lowing abbreviations in Eqs. (3)-(5): c ± ( k ) = c y ( k ) ± ic x ( k ), ˜ M ± ( k ) = M ( k ) ± i [ B c z ( k ) + R d ( k )], D k = q A [ c x ( k ) + c y ( k )] + R d ( k ), E k = q | ˜ M ± ( k ) | + D k ,˜ E k = E k + D k , N k = p E k ˜ E k , D ± ( k ) = D k ± R d ( k ),˜ D k = p D k D − ( k ), and c + ( k ) = i q c x ( k ) + c y ( k ) e − iϕ k = ic ( k ) e − iϕ k , (6) W k = ˜ E k + ˜ M + ( k ) √ N k = | W k | e iφ k , (7) R d ( k ) + iA c ( k ) = D k e iθ k . (8)The above formulae define the pseudospin basis for theconduction band states on the northern hemisphere of theBZ ( k z > k z < | k , α i = P |− k , α i = − T | − k , β i and | k , β i = P | − k , β i = T | − k , α i . We in-troduce the new Pauli matrices ̺ i ( i = 1 , ,
3) and thecorresponding unit matrix ̺ in the subspace of the twopseudospin bases. The reduced model containing onlystates of the conduction band is simply h ( k ) = E ( k ) ̺ , (9)where the dispersion of the conduction band is E ( k ) = ǫ ( k ) + E k . III. LISTS AND GENERAL PROPERTIES OFTWO-COMPONENT PAIRINGS
We now construct the full lists of the basis functionsfor the E g and E u representations of the D d group, upto 2NN in-plane pairing correlations, consistent with thetight-binding model which is up to 2NN in-plane hop-ping terms [36]. Corresponding to the two basis of themodel defined in the previous section, there are two waysof classifying the possible pairing channels in M x Bi Se .The first approach focuses on the low-energy states closeto the Fermi surface [23, 24, 28]. For the x > x Bi Se super-conductor in the bulk, or if the normal phase is topo-logically trivial, the two approaches give essentially thesame results. However, if we are also interested in thetopological aspect of the system inherited from the topo-logically nontrivial normal phase, such as the coexistenceof the topological surface states with the Fermi surface[2, 51, 52], then it is advantageous, if not imperative, towork with the second approach.We will first construct in the spin-orbital basis the fulllists of pairing channels in both the E g and the E u repre-sentations, up to 2NN in-plane pairing correlations. Basisfunctions for the irreducible representations of the D d symmetry group can be constructed in terms of the Γmatrices or the symmetrized Fourier functions. Here, wedefine the ‘symmetrized Fourier functions’ as linear com-binations of the trigonometric functions cos( k · l ) andsin( k · l ), where k is the wave vector and l represents anNN or 2NN bond vectors defined in section II. Full listsof these basis functions exist in previous works [36, 37]. To be self-contained, we include them in Table I and Ta-ble II. For each representation, there are two sets of basisfunctions in terms of the Γ matrices. Up to a constantnumber of unit module, the two basis sets differ by a fac-tor of Γ . This is easy to understand from the fact thatΓ belongs to the A g representation, which respects thefull symmetry of the crystal and maps an existing basisset to a new basis set belonging to the same representa-tion. For the E g and E u representations, the two basissets in Table I transform in the same manner under the D d group [36]. The symmetrized Fourier functions inTable II and their expansions in the limit of small in-plane wave vectors (i.e., k x a ≃ k y a ≃
0) are ϕ ( k ) = 13 [cos k · δ + cos k · δ + cos k · δ ] ≃ −
16 ( k x + k y ) a , (10) ϕ ( k ) = 12 [cos k · δ − cos k · δ ] ≃ − √ k x k y a , (11) ϕ ( k ) = [cos k · δ + cos k · δ − k · δ ]2 √ ≃ − √
38 ( k x − k y ) a , (12) ϕ ( k ) = d ( k ) ≃ ( k x − k x k y ) a = 12 ( k + k − ) a , (13) ϕ ( k ) = d ( k ) ≃ (3 k y k x − k y ) a = 12 i ( k − k − ) a , (14) ϕ ( k ) = c x ( k ) ≃ k x a, (15)and ϕ ( k ) = c y ( k ) ≃ k y a. (16)We have introduced the abbreviation k ± = k x ± ik y . Ifwe extend to include the inter-quintuple-layer pairings,we can replace ϕ ( k ) with ϕ ′ ( k ) = cos k · δ and replace ϕ ( k ) with ϕ ′ ( k ) = sin k · δ . However, we will focus onthe intra-quintuple-layer pairings in this work.By multiplying the basis functions in Table I and thosein Table II, we can get various product representations ofthe D d group. These product representations can be de-composed into the irreducible representations accordingto the group theory [53]. For example, A u ⊗ E g = E u and E g ⊗ E u = A u + A u + E u . In such a manner, we canidentify all the realizations of the possible irreducible rep-resentations. When taken as a part of the model Hamilto-nian, they are subject to further constraints. If a term istaken as a part of the model for the electronic structuresin the normal state, this term should belong to the A g representation and has to be Hermitian [36, 37]. These TABLE I: Basis functions in terms of the Γ matrices. Thesymbols in the brackets of the first column are anothercommonly used name for the corresponding representation[36, 37]. The semicolons in the second column separate dif-ferent basis sets of the same representation.Representation Basis A g (˜Γ +1 ) I ; Γ A g (˜Γ +2 ) Γ ; Γ E g (˜Γ +3 ) { Γ ,Γ } ; { Γ ,Γ } A u (˜Γ − ) Γ ; Γ A u (˜Γ − ) Γ ; Γ E u (˜Γ − ) { Γ ,Γ } ; { Γ ,Γ } TABLE II: Basis functions in terms of the symmetrizedFourier functions. The symmetrized Fourier functions are de-fined as linear combinations of cos( k · l ) and sin( k · l ), where k is the wave vector and l represents an NN or 2NN bondvectors defined below Eq.(1). The symbols in the bracketsof the first column are another commonly used name for thecorresponding representation [36, 37]. The semicolon in thesecond line of the second column separates different basis setsof the A g representation.Representation Basis A g (˜Γ +1 ) 1; ϕ ( k ) A g (˜Γ +2 ) none E g (˜Γ +3 ) { ϕ ( k ), ϕ ( k ) } A u (˜Γ − ) ϕ ( k ) A u (˜Γ − ) ϕ ( k ) E u (˜Γ − ) { - ϕ ( k ), ϕ ( k ) } constraints, together with the time-reversal symmetry ofthe materials, eliminate many combinations.If a basis set is taken as the superconducting pairingterm, then it has to satisfy the Fermi exchange statis-tics. Another important aspect about the symmetry ofthe pairing term is related to the peculiar transformationproperty of the pairing term under the time-reversal op-eration [54, 55]. In the spin-orbital basis φ † k , the pairingterm has the following general expression φ † k ∆( k )[ φ †− k ] T + H.c. , (17)where the superscript ‘T’ means taking the transpose,and H.c. means the Hermitian conjugate of the first term.According to Blount [54, 55], the time-reversed creationoperator transforms under the symmetry operation justlike the corresponding annihilation operator. That is, T [ φ † k ] T = − iσ ⊗ s [ φ †− k ] T (18)transforms in the same manner as φ k under the actionof the D d symmetry group and the time-reversal opera-tion [54]. This means that, by transforming the creation operator part of the pairing term to φ † k ∆( k )[ φ †− k ] T = φ † k ∆( k ) iσ ⊗ s ( − iσ ⊗ s )[ φ †− k ] T , (19)the matrix ∆( k ) iσ ⊗ s (20)has the same transformation property as the terms inthe model for the normal state electronic structures. Asa result, we can construct the basis functions according tothe general procedure for the normal state [37, 53]. Sincethe obtained basis is of the form of Eq.(20), we multiply − iσ ⊗ s from the right and get the basis functions of thepairing term. Then we single out from the results thoseobeying the Fermi exchange statistics, namely ∆ T ( − k ) = − ∆( k ). The above discussions correspond to taking theNambu basis as [ φ † k , φ T − k ]. If we take the Nambu basis as[ φ † k , φ T − k ( iσ ⊗ s )] instead, then the pairing term will bein the form of Eq.(20) spontaneously [44]. Hereafter, wewill stick to the first Nambu basis.The resulting basis sets (up to 2NN in-plane pairingcorrelations) for the E g and E u representations are sep-arately shown in Table III and Table IV. The new Γ-matrices in the tables with two subindices are defined asΓ µν = i [Γ µ , Γ ν ], where both µ and ν run from 1 to 5.Explicitly, Γ = σ ⊗ s , Γ = − σ ⊗ s , Γ = σ ⊗ s ,Γ = σ ⊗ s , Γ = σ ⊗ s , Γ = σ ⊗ s , Γ = σ ⊗ s ,Γ = σ ⊗ s , Γ = σ ⊗ s , and Γ = σ ⊗ s . Notethat, we have multiplied a factor of σ ⊗ is to each com-ponent of the basis sets. To get the final expressions forthe pairing components, we have to multiply back a fac-tor of σ ⊗ ( − is ) to each component listed in the tables[44, 54–56]. Also notice that, each basis set can be mul-tiplied by a factor of an arbitrary linear combination of ϕ ( k ) and a constant. TABLE III: Basis functions for the even-parity two-dimensional representation E g , expressed as linear combina-tions of products between the Γ matrices and the symmetrizedFourier functions. The first column is the numbering of thevarious pairing channels. The second and third columns areseparately the two components of the corresponding basissets. E ( n ) g ψ ( n )1 ( k )( σ ⊗ is ) ψ ( n )2 ( k )( σ ⊗ is ) n = 1 - I ϕ ( k ) I ϕ ( k ) n = 2 -Γ ϕ ( k ) Γ ϕ ( k ) n = 3 -Γ ϕ ( k ) Γ ϕ ( k ) n = 4 Γ ϕ ( k ) -Γ ϕ ( k ) n = 5 Γ ϕ ( k ) Γ ϕ ( k ) n = 6 Γ ϕ ( k ) Γ ϕ ( k ) n = 7 Γ ϕ ( k )+Γ ϕ ( k ) -Γ ϕ ( k )+Γ ϕ ( k ) The symmetry channels listed in Table III and TableIV are one central result of this work. Two featuresof the Tables are apparent. Firstly, only the two basis
TABLE IV: Basis functions for the odd-parity two-dimensional representation E u , expressed as linear combina-tions of products between the Γ matrices and the symmetrizedFourier functions. The first column is the numbering of thevarious pairing channels. The second and third columns areseparately the two components of the corresponding basissets. E ( n ) u ˜ ψ ( n )1 ( k )( σ ⊗ is ) ˜ ψ ( n )2 ( k )( σ ⊗ is ) n = 1 Γ Γ n = 2 Γ ϕ ( k ) Γ ϕ ( k ) n = 3 Γ ϕ ( k ) -Γ ϕ ( k ) n = 4 Γ ϕ ( k )+Γ ϕ ( k ) Γ ϕ ( k )-Γ ϕ ( k ) n = 5 Γ ϕ ( k ) Γ ϕ ( k ) n = 6 Γ ϕ ( k ) Γ ϕ ( k ) n = 7 Γ ϕ ( k ) Γ ϕ ( k ) n = 8 -Γ ϕ ( k ) Γ ϕ ( k ) n = 9 Γ ϕ ( k ) -Γ ϕ ( k ) n = 10 Γ ϕ ( k ) Γ ϕ ( k ) n = 11 − Γ ϕ ( k )+Γ ϕ ( k ) Γ ϕ ( k )+Γ ϕ ( k ) n = 12 Γ ϕ ( k )+Γ ϕ ( k ) Γ ϕ ( k )-Γ ϕ ( k ) functions of E (1) u are completely k -independent. This isthe pairing channel that has attracted the most atten-tion as a promising candidate for the nematic pairing inCu x Bi Se [19, 22]. Secondly, among the listed basis setsin the E g channel, only the leading three channels arespin-singlet in the spin-orbital basis. Among the twelve E u channels, only E (3) u is spin-singlet in the spin-orbitalbasis.According to a theorem by Yip and Garg, the mostgeneral pairing can be written as a linear combination ofall independent basis sets of the representation [57]. Wetherefore write the general expression of the pairing termin the E g representation as∆ g ( k ) = X α =1 ∆ α [ η ψ ( α )1 ( k ) + η ψ ( α )2 ( k )] . (21)( η , η ) is the same vector for all seven E g pairing chan-nels, so that the order parameter transforms as a well-defined vector in the subspace of the E g representation.Similarly, the general pairing in the E u representation iswritten as∆ u ( k ) = X α =1 ˜∆ α [˜ η ˜ ψ ( α )1 ( k ) + ˜ η ˜ ψ ( α )2 ( k )] . (22)Again, (˜ η , ˜ η ) are the same vector for all twelve E u pairing channels. Notice that, we have restricted thepairing to the same (i.e., E u or E g ) irreducible repre-sentation of the symmetry group. The superconductingstate with mixed even-parity and odd-parity components,which was found to stabilize under suitable circumstances[28, 58], will not be considered in the present work. Among the pairing combinations contained in Eqs.(21)and (22), we are particularly interested in those pairingsthat have been stabilized as the ground state in previ-ous studies based on a microscopic or phenomenologicalpairing mechanism, and those that are possibly consis-tent with more than one key experimental consensuseson the M x Bi Se superconductors. The theoretical stud-ies motivated by the recent experiments focus on pairingsin the odd-parity E u channel. Besides E (1) u [19, 31], other E u pairings have been studied in previous works. Directcomparison shows that E (6) u and E (11) u were studied byYuan et al [26], E (6) u and E (9) u (if we replace ϕ ( k ) with ϕ ′ ( k ) = sin k · δ ) were studied by Chirolli et al [27].These pairing channels exhaust the three kinds of pair-ings in Table I and Table II, as regards the differencebetween the two pairing components: For E (1) u and E (9) u the difference comes from the two different Γ-matrices,Γ versus Γ for E (1) u and Γ versus Γ for E (9) u ; for E (6) u the difference between the two components comesfrom the ϕ ( k ) and ϕ ( k ) symmetry factors; for E (11) u the distinction comes from a combination of the differ-ence between Γ and Γ and the difference between ϕ ( k ) and ϕ ( k ).The even-parity E g channels have attracted much lessattention. The only theoretical paper focusing on theeven-parity pairings studied the leading pairing insta-bilities resulting from the purely repulsive short-rangeelectron-electron interactions [36, 56]. The six pairingsfound in that paper could be identified as the six pair-ing components of E (1) g , E (2) g , and E (3) g . In all these three E g channels, the difference between the two basis compo-nents comes completely from the k -dependent symmetryfactors.The nature of the pairing defined by Eq.(21) [Eq.(22)]depends both on the pairing strengths ∆ α ( ˜∆ α ) and onthe two-component vector ( η , η ) [(˜ η , ˜ η )]. For sim-plicity, we will assume in the following analysis that ∆ α ( α = 1 , ...,
6) and ˜∆ α ( α = 1 , ...,
12) are all real numbers.By this convention, we neglect pairings analogous to thesingle-component chiral pairings, like the s + is ′ , d + id ′ ,and p + ip ′ pairings [59–61]. The nonzero ∆ α or ˜∆ α indicate the pairing channels that contribute to the su-perconducting order parameter. The relative magnitudesof the ∆ α or ˜∆ α parameters characterize the contribu-tions of different pairing channels. Then, depending onthe ( η , η ) or the (˜ η , ˜ η ) vector, we may define thechirality and nematicity of the two-component supercon-ducting order parameters [19, 23, 24]: The pairing is ne-matic if at least one of | η | − | η | ( | ˜ η | − | ˜ η | ) and η η ∗ + η ∗ η (˜ η ˜ η ∗ + ˜ η ∗ ˜ η ) is nonzero. The pairing is chi-ral if η η ∗ − η ∗ η = 0 (˜ η ˜ η ∗ − ˜ η ∗ ˜ η = 0), or equivalently η /η (˜ η / ˜ η ) is a nonzero and finite complex number.To understand the properties of various pairings, weturn from the spin-orbital basis to the band (pseudospin)basis [44–50]. Defining U k = [ | k , α i , | k , β i ] and U − k =[ | − k , α i , | − k , β i ], for k z >
0. A pairing wave func-tion (i.e., superconducting order parameter) expressedas ∆( k ) in the original spin-orbital basis transforms to˜∆( k ) = U † k ∆( k ) U ∗− k (23)in the pseudospin basis [31]. As the properties of thesix symmetry factors ϕ i ( k ) ( i = 1 , ...,
6) are known fromEqs.(11)-(16), the remaining task is to calculate the basistransformations for the sixteen 4 × A. The E g pairings First consider the E g representation. Define I ′ = I ( − σ ⊗ is ) and Γ ′ i = Γ i ( − σ ⊗ is ) ( i = 1 , ..., ̺ i ( i = 0 , , , I ′ ( k ) = − i̺ , (24)˜Γ ′ ( k ) = − A c y ( k ) E k i̺ , (25)˜Γ ′ ( k ) = A c x ( k ) E k i̺ , (26)˜Γ ′ ( k ) = − R d ( k ) E k i̺ , (27)˜Γ ′ ( k ) = − B c z ( k ) + R d ( k ) E k i̺ , (28)˜Γ ′ ( k ) = − M ( k ) E k i̺ . (29)While only three (i.e., E (1) g , E (2) g , and E (3) g ) out ofthe seven channels listed in Table III are spin-singlet inthe original spin-orbital basis, all seven E g channels arepseudospin singlets in the pseudospin basis. In particu-lar, although E (7) g is very different from E (1) g and E (2) g in the spin-orbital basis, the symmetry factors of the twobasis components of E (7) g behave like k x − k y and − k x k y in the band (pseudospin) basis, qualitatively the same asthe corresponding basis components of E (1) g and E (2) g , ifthe slight anisotropy introduced by E k and M ( k ) are ne-glected. In addition, E (4) g is identical to E (6) g in the pseu-dospin basis, up to a k -independent constant factor. Fi-nally, we point out that the seemingly different E (3) g and E (5) g are closely related. In fact, if we replace ϕ ( k ) = d ( k ) in E (5) g by ϕ ( k ) + B R ϕ ′ ( k ) = d ( k ) + B R c z ( k ),which have the same symmetry as that of ϕ ( k ) under D d , then E (3) g and E (5) g are identical in the pseudospinbasis up to a constant factor. Notice that, the six spin-singlet pairings identified in aprevious study combine to E (1) g , E (2) g , and E (3) g [36, 56].The two components of E (3) g were incorrectly identifiedas belonging to the E u representation [36, 56], because ofneglecting the different transformation properties of thepairing term compared to the model for the normal stateelectronic structures, as we explained in Eqs.(17)-(20).To understand the qualitative properties of the vari-ous pairing channels, we consider each pairing channelseparately. That is, we assume that only one amongthe seven ∆ α parameters in Eq.(21) is nonzero. Fromthe transformations in Eqs.(24)-(29) and the expansionsin Eqs.(11)-(16), it is easy to see that all the fourteenpairing components, ψ ( n ) m ( n = 1 , ..., m = 1 , ψ (1)1 and ψ (1)2 can give a fully-gapped bulk spectrum. Thesame is true for the two components of E (2) g and the twocomponents of E (7) g . If the Fermi surface is spheroidal,the chiral combination of the two components of E (1 , , g give a bulk spectrum with two point nodes at k x = k y = 0of the Fermi surface. The remaining four pairing chan-nels, E (3) g to E (6) g , have line nodes for arbitrary ( η , η ),for both spheroidal and corrugated cylindrical Fermi sur-faces. One common set of line nodes for E (3) g comesby setting the B c z ( k ) + R d ( k ) factor of Eq.(28) tozero. Six line nodes persist for E (4) g and E (6) g , which bothcome from ϕ ( k ) = d ( k ) = 0. E (5) g also has six preva-lent line nodes, which come from ϕ ( k ) = d ( k ) = 0.For a spheroidal Fermi surface, the six line nodes from d ( k ) = d ( k ) = 0 connect at the two points of the Fermisurface with k x = k y = 0.We can also estimate the magnitude of the supercon-ducting gaps. For M x Bi Se (M is Cu, Sr, or Nd), thechemical potential µ > k x a and k y a are all very small for wavevectors on the Fermi surface. The superconducting gapof a certain pairing channel can therefore be character-ized in terms of its power in ka = q k x + k y a . For stateslying on the Fermi surface, we have E k + ǫ ( k ) = µ . ǫ ( k )and M ( k ) vary only slightly over states on the Fermi sur-face, and so does E k [31]. As an approximation, we treat E k = µ − ǫ ( k ) and M ( k ) as constants. Under these con-ditions, we see that the gap of E (1) g is of the order ( ka ) . E (2) g and E (7) g also open superconducting gaps in the or-der of ( ka ) , but reduced by a factor of M ( k ) /E k and A /E k compared to E (1) g . The superconducting gaps of E (4 , , g are all in the order of ( ka ) and are two powerssmaller than ( ka ) . For c z ( k ) = 0, the superconductinggap for E (3) g is also in the order of ( ka ) . However, thetwo components of E (3) g behave more like ( k x a )( k z c ) and( k y a )( k z c ), and are more efficient than E (4 , , g in open-ing the superconducting gap. The angular dependenceof the pairing amplitudes on the k x k y plane, which wasneglected in the above analysis in terms of the ka factor,can be obtained from Eqs.(11)-(16).Among the seven pairing channels in Table III, onlythe two components of E (3) g can have a sign change in thepseudospin basis, when we substitute − k z for k z . Thisimplies that only E (3) g can give SABSs on the natural xy surface of the M x Bi Se superconductors. As regards thetopological surface states of the normal phase, accordingto a previous theoretical study [62], all the pairing chan-nels in Table III except for E (3) g can open a gap in thetopological surface states. Therefore, among all seven E g pairings, E (3) g is special as regards the surface properties.The surface states for the E (3) g pairings were studied ina previous work [36]. On the other hand, since all seven E g pairings are pseudospin singlets, they are expected tohave trivial isotropic electronic spin susceptibility in the xy plane. B. The E u pairings We next study the pairings belonging to the E u rep-resentation. We define Γ ′ µν = Γ µν ( − iσ ⊗ s ) for µ, ν =1 , ..., µ < ν . As have been noticed in previousstudies, the E (1) u channel, which is k -independent in thespin-orbital basis, has a complicated k -dependence in theband (pseudospin) basis [19, 31, 44]. This is generallytrue for all the 12 E u channels listed in Table II. Besidesthe E (1) u channel, only the symmetry factors of E (2) u and E (3) u are even functions of k . In addition, E (2) u is a directproduct of Γ which belongs to the A u representation(Table I) and { ϕ ( k ) , ϕ ( k ) } which belongs to the E g representation (Table II). Γ ′ was a chief candidate ofthe pairing for Cu x Bi Se in early theoretical discussions[32]. The relevant basis transformations for these threechannels are˜Γ ′ ( k ) = [( B c z + R d ) ̺ − R d ̺ − A c x ̺ ] i̺ E k , (30)˜Γ ′ ( k ) = [ R d ̺ + ( B c z + R d ) ̺ − A c y ̺ ] i̺ E k , (31)˜Γ ′ ( k ) = [ A c x ̺ + A c y ̺ + ( B c z + R d ) ̺ ] i̺ E k , (32)˜Γ ′ ( k ) = [ − A c y ̺ + A c x ̺ − R d ̺ ] i̺ E k . (33)The k -dependence of the terms in the results are sup-pressed to simplify the notations. In contrast to thesalient two-fold anisotropy in the spin structure factorsin the two bases for E (1) u , the spin structure factors for E (2) u and E (3) u are fairly symmetric in the xy -plane. The symmetrized Fourier functions for the remaining E u channels are all odd functions of k . The propertiesof these pairing channels can also be understood by com-bining the symmetry factors and the expressions of theremaining six Γ-matrices in the pseudospin basis. Bystraightforward applications of Eq.(23), we get the fol-lowing results:˜Γ ′ ( k ) = {− [ A c x ( B c z + R d ) + A c y R d ] ̺ − [ A c y ( B c z + R d ) − A c x R d ] ̺ − [ E k ( E k + M ) − A c ] ̺ } i̺ E k ( E k + M ) , (34)˜Γ ′ ( k ) = { [ A c x ( B c z + R d ) + A c y R d ] ̺ +[ A c y ( B c z + R d ) − A c x R d ] ̺ (35) − [ M ( E k + M ) + A c ] ̺ } i̺ E k ( E k + M ) . ˜Γ ′ ( k ) = { [ R d ( B c z + R d ) − A c x c y ] ̺ +[ E k ( E k + M ) − A c y − R d ] ̺ (36) − [ A c x R d + A c y ( B c z + R d )] ̺ } i̺ E k ( E k + M ) , ˜Γ ′ ( k ) = {− [ E k ( E k + M ) − A c x − R d ] ̺ +[ A c x c y + R d ( B c z + R d )] ̺ (37)+[ A c x ( B c z + R d ) − A c y R d ] ̺ } i̺ E k ( E k + M ) , ˜Γ ′ ( k ) = {− [ E k ( E k + M ) − A c y − ( B c z + R d ) ] ̺ − [ A c x c y + R d ( B c z + R d )] ̺ (38) − [ A c x ( B c z + R d ) − A c y R d ] ̺ } i̺ E k ( E k + M ) , ˜Γ ′ ( k ) = {− [ A c x c y − R d ( B c z + R d )] ̺ +[ A c x + ( B c z + R d ) − E k ( E k + M )] ̺ (39) − [ A c x R d + A c y ( B c z + R d )] ̺ } i̺ E k ( E k + M ) . Again, the k -dependence of the terms in the results aresuppressed. The expression multiplying ̺ in the thirdline of Eq.(35) can be rewritten as M ( E k + M )+ A c = E k ( E k + M ) − R d − ( B c z + R d ) . The transformations in Eqs.(30)-(39) and those inEqs.(24)-(29) are another central result of the presentwork.There is an interesting relation among the six basistransformations in Eqs. (34)-(39). According to Table I,the six primed Γ-matrices separate into three pairs: Γ ′ and Γ ′ = Γ ′ Γ , Γ ′ and Γ ′ = Γ ′ Γ , Γ ′ and Γ ′ =Γ ′ Γ . The connection between the two components ofeach pair is revealed by summing over the correspondingexpressions in the pseudospin basis, which gives˜Γ ′ ( k ) + ˜Γ ′ ( k ) = − E k + ME k ̺ ( i̺ ) , (40)˜Γ ′ ( k ) + ˜Γ ′ ( k ) = − E k + ME k ̺ ( i̺ ) , (41)˜Γ ′ ( k ) + ˜Γ ′ ( k ) = − E k + ME k ̺ ( i̺ ) . (42)Note that I + Γ = I + P , which underlies the abovecombinations, is the projection operator to the even-parity subspace in the spin-orbital basis. Eqs. (40)-(42)have a clear cyclical structure. This should be relatedto the definition of the pseudospin basis, which was cho-sen to make the even-parity component of the magneticmoment operator to transform like a proper axial vector[31]. The combinations corresponding to the projectionoperator I − Γ = I − P , which projects to the odd-parity subspace, are˜Γ ′ ( k ) − ˜Γ ′ ( k )= {− A c x ( B c z + R d ) + A c y R d ] ̺ − A c y ( B c z + R d ) − A c x R d ] ̺ (43)+[ A c − R d − ( B c z + R d ) ] ̺ } i̺ E k ( E k + M ) , ˜Γ ′ ( k ) − ˜Γ ′ ( k )= { [ A ( c x − c y ) + R d − ( B c z + R d ) ] ̺ +2[ A c x c y + R d ( B c z + R d )] ̺ (44)+2[ A c x ( B c z + R d ) − A c y R d ] ̺ } i̺ E k ( E k + M ) , ˜Γ ′ ( k ) − ˜Γ ′ ( k )= { A c x c y − R d ( B c z + R d )] ̺ +[ A ( c y − c x ) + R d − ( B c z + R d ) ] ̺ (45)+2[ A c x R d + A c y ( B c z + R d )] ̺ } i̺ E k ( E k + M ) . Although more complex than Eqs.(40)-(42), Eqs.(43)-(45) are simpler than Eqs.(34)-(39) because the E k ( E k + M ) factors are eliminated from the numerators.An implication of Eqs.(40)-(45) is that, we can rehy-bridize four pairs of the E u channels in Table IV to sim-plify the discussions on the properties of those pairingchannels. Explicitly, the four pairs include { E (5) u , E (6) u } , { E (7) u , E (8) u } , { E (9) u , E (10) u } , and { E (11) u , E (12) u } . The eightnew pairing combinations are defined in Table V. Theexpressions in the pseudospin basis of the eight pairingchannels in Table V follow directly from the definitions in Table IV and Eqs.(40)-(45), and will not be written outexplicitly. In shifting from Table IV to Table V, Eq.(22)is also changed in a straightforward manner. For exam-ple, the pairing components corresponding to E (5) u and E (6) u are changed to E (1 p ) u and E (1 m ) u by X α =5 ˜∆ α [˜ η ˜ ψ ( α )1 ( k ) + ˜ η ˜ ψ ( α )2 ( k )]= m X α ′ =1 p ˜∆ α ′ [˜ η ˜ ψ ( α ′ )1 ( k ) + ˜ η ˜ ψ ( α ′ )2 ( k )] , (46)where the (˜ η , ˜ η ) vector does not change, ˜∆ p = ( ˜∆ +˜∆ ) / m = ( ˜∆ − ˜∆ ) /
2. Hereafter, we consider E (1) u to E (4) u in Table IV and the eight pairings in TableV as the 12 independent E u pairing channels. TABLE V: Redefinitions of the basis functions for the eightchannels of the E u representation in Table IV, from E (5) u to E (12) u . E ( n ′ ) u ˜ ψ ( n ′ )1 ( k ) ˜ ψ ( n ′ )2 ( k ) n ′ = 1 p ˜ ψ (5)1 ( k ) + ˜ ψ (6)1 ( k ) ˜ ψ (5)2 ( k ) + ˜ ψ (6)2 ( k ) n ′ = 1 m ˜ ψ (5)1 ( k ) − ˜ ψ (6)1 ( k ) ˜ ψ (5)2 ( k ) − ˜ ψ (6)2 ( k ) n ′ = 2 p ˜ ψ (7)1 ( k ) + ˜ ψ (8)1 ( k ) ˜ ψ (7)2 ( k ) + ˜ ψ (8)2 ( k ) n ′ = 2 m ˜ ψ (7)1 ( k ) − ˜ ψ (8)1 ( k ) ˜ ψ (7)2 ( k ) − ˜ ψ (8)2 ( k ) n ′ = 3 p ˜ ψ (9)1 ( k ) + ˜ ψ (10)1 ( k ) ˜ ψ (9)2 ( k ) + ˜ ψ (10)2 ( k ) n ′ = 3 m ˜ ψ (9)1 ( k ) − ˜ ψ (10)1 ( k ) ˜ ψ (9)2 ( k ) − ˜ ψ (10)2 ( k ) n ′ = 4 p ˜ ψ (11)1 ( k ) + ˜ ψ (12)1 ( k ) ˜ ψ (11)2 ( k ) + ˜ ψ (12)2 ( k ) n ′ = 4 m ˜ ψ (11)1 ( k ) − ˜ ψ (12)1 ( k ) ˜ ψ (11)2 ( k ) − ˜ ψ (12)2 ( k ) Now we make an order of magnitude estimation overthe superconducting gap amplitudes of the various pair-ing channels in Table IV and Table V, by taking advan-tage of the basis transformations of the Γ-matrices inEqs.(30)-(45) and the expansions of the symmetry fac-tors in Eqs.(11)-(16). As we have explained in the previ-ous subsection, k x a and k y a are all very small for wavevectors on the Fermi surface. The pairing amplitude of apairing channel can be characterized in terms of its powerin ka = q k x + k y a . Also, as an approximation, we maytreat E k = µ − ǫ ( k ) and E k + M ( k ) as constants. Underthe above conditions, the gap of E (1) u is in the order of ka ,the gaps of E (2 , , u are in the order of ( ka ) . While E (1 p ) u and E (4 p ) u open gaps in the the order of ka , E (1 m ) u and E (4 m ) u open gaps in the order of ( ka ) . E (2 p ) u and E (3 p ) u also open gaps in the order of ( ka ) . The gap opened by E (2 m ) u and E (3 m ) u are of the order ( ka ) . The dependenceof the superconducting gap amplitudes on the directionsof the wave vectors on the k x k y plane can be read fromthe expansions in Eqs.(11)-(16). Notice that, if the Fermisurface is corrugated cylindrical, the c z ( k ) factor under-goes drastic variations in scanning over the Fermi surfaceand may dominate the variation of the pairing amplitude[31]. The above analysis is then an estimate of the bulkpairing magnitudes on the k z = 0 and k z = ± π planes ofthe BZ. From the above analysis, E (1) u , E (1 p ) u and E (4 p ) u are the most efficient in opening a large superconductinggap. In contrast, the E (2 m ) u and E (3 m ) u channels are leasteffective in producing a superocnducting gap, and are un-likely to be the leading pairing channels. The remainingseven E u pairing channels have the intermediate abilityin opening a bulk superconducting gap.The surface states on the natural xy surface of super-conducting M x Bi Se may include the SABSs and thetopological surface states. The SABSs are directly re-lated to the superconducting order parameter, and exist(do not exist) if the superconducting pairing wave func-tion undergoes (does not undergo) a sign reversal uponthe reflection changing k z to − k z . By this rule, it iseasy to see that only five pairing channels do not sup-port SABSs on the xy surface, including E (3) u , E (1 p ) u , E (2 p ) u , E (3 p ) u , and E (4 p ) u . The topological surface statesinherited from the normal state may or may not opena gap at the Fermi level in the superconducting phase,depending both on the nature of the bulk pairing and onthe nature of the topological surface states. A full listof the pairings that open a gap in the topological sur-face states of Bi Se was constructed previously [40, 62].According to that list, E (1) u , E (2) u , and E (4) u cannot gapthe topological surface states. On the other hand, all theremaining pairings can at least partially gap the topo-logical surface states. For example, ˜ ψ (3)1 can gap all thetopological surface states at the Fermi level, except forthe crossing points between the topological surface statesat the Fermi level and the lines determined by ϕ ( k ) = 0.To infer the qualitative behaviors of the electronic spinsusceptibility relevant to the Knight shift experiment, wereformulate the pairing defined by Eq.(23) in the stan-dard expression as˜∆( k ) = [ d ( k ) ̺ + d ( k ) · ̺ ] i̺ , (47)where d ( k ) is the pseudospin-singlet component ofthe pairing, d ( k ) = ( d ( k ) , d ( k ) , d ( k )) is a three-component vector for the pseudospin-triplet part of thepairing, and ̺ = ( ̺ , ̺ , ̺ ). From the above results,we have d ( k ) = (0 , ,
0) for the even-parity E g pairingsand d ( k ) = 0 for the odd-parity E u pairings. The d ( k )vector of an E u pairing is perpendicular to the spin direc-tion of the corresponding Cooper pair [63–66]. The spinsusceptibilities show distinct behaviors depending on therelative orientation between the external magnetic field H and the d ( k ) vector [15, 47, 63–66]: For H ⊥ d ( k ), thespin susceptibility barely changes across the supercon-ducting transition temperature T c ; for H k d ( k ), the spinsusceptibility decreases below T c similar to a spin-singletsuperconductor. In this manner, we can understand thequalitative behaviors of the electronic spin susceptibili-ties for the E u pairing channels in Table IV and TableV, in terms of Eqs.(30)-(45). For simplicity and without losing much rigor, we can first set R = R = 0 and workwith the effective pairing in the simplified model. It iseasy to see that both of the two components of the E (1) u channel break the spin rotational symmetry in the xy -plane to two-fold symmetry. It is also clear that both E (2) u and E (3) u respect the spin rotational symmetry in the xy plane. For E (4) u , while the two Γ-matrices that are con-tained in each channel breaks the spin rotational symme-try in the xy -plane, the linear combination restores thissymmetry. E (1 p ) u and E (1 m ) u respect the spin rotationalsymmetry in the xy -plane. The two components of E (2 p ) u and E (3 p ) u both break the three-fold in-plane spin rota-tional symmetry in the xy -plane to two-fold symmetry.The two components of E (2 m ) u and E (3 m ) u only slightlybreak the three-fold in-plane spin rotational symmetryin the xy -plane. The components of E (4 p ) u and E (4 m ) u respect the spin rotational symmetry in the xy -plane. IV. COMPARISON TO EXPERIMENTALRESULTS
In this section, we first make a survey over the mainexperimental consensuses on the three superconductors.Then we combine the results in the previous section toinfer the most probable pairings for the three supercon-ductors. To be clear, we categorize the relevant experi-ments on the M x Bi Se (M is Cu, Sr, or Nd) supercon-ductors into three broad classes, which separately probethe bulk spectrum, the surface spectrum (along the nat-ural xy surface parallel to the basal plane), and the mag-netic properties. Since both the bulk spectrum and thesurface spectrum influence the magnetic properties, thisdivision is by no means absolute.For the bulk spectrum, both Cu x Bi Se [12, 67] andSr x Bi Se [68] are reported to be fully-gapped. However,there is no solid consensus on the momentum-dependenceof the band gap. For Cu x Bi Se , while an STM experi-ment gives a tunneling spectrum consistent with isotropic s -wave pairing [12], a later field-angle-dependent specificheat experiment suggests an energy gap structure withsalient two-fold symmetry in the k x k y plane [16]. ForSr x Bi Se , an STM experiment implies an anisotropic s -wave pairing [68]. Later, an experiment shows with resis-tivity and Laue diffraction measurements that the energygap structure of Sr x Bi Se also has two-fold symmetry inthe k x k y plane [17]. For Nd x Bi Se , an experiment findsevidence for point nodes in its bulk spectrum [69].For the surface spectrum, the situation is rather con-fusing, for all three superconductors. While several earlyexperiments, in particular those based on point-contactspectroscopy, claim to have found evidence of SABSs forCu x Bi Se [7–9], later experiments denied the existenceof SABSs [12, 13]. Similar confusion exists for Sr x Bi Se .An experiment infers the existence of surface states inSr x Bi Se , through the Shubnikov-de Haas oscillations[3]. However, an STM/STS experiment does not see0any in-gap states [68]. For superconducting Nd x Bi Se ,the angle-resolved photoemission spectroscopy experi-ment by Qiu et al shows that the topological surfacestates in the normal state is preserved in the supercon-ducting state [6]. But it is unclear whether or not thetopological surface states are gapped at the Fermi level,and whether or not there are SABSs on the xy surface ofNd x Bi Se .For the magnetic properties, Nd x Bi Se was reportedto have spontaneous magnetization in the superconduct-ing phase [6], which seems to be related to the magnetismof the Nd dopants. The Cu x Bi Se and Sr x Bi Se su-perconductors are usually considered as nonmagnetic inthe absence of external magnetic fields [6]. On the otherhand, for all three superconductors, there are magnetism-related experiments indicating that the three-fold rota-tional symmetry in the normal phase is broken downto two-fold rotational symmetry in the superconduct-ing phase. The Knight shift experiment by Matano et al shows that the electronic spin susceptibilities ofCu x Bi Se is two-fold symmetric in the basal plane andis invariant for an out-of-plane magnetic field [15]. An-other magnetic measurement is the upper critical field, B c . By varying the direction of the magnetic field inthe basal plane of the superconductors, a two-fold sym-metry in the upper critical field is observed for Cu x Bi Se [16], for Sr x Bi Se [18], and also for Nd x Bi Se [21]. ForNd x Bi Se , the two-fold rotational symmetry in the basalplane is also confirmed in a torque magnetometry mea-surement [20].In the light of the above survey, the most well-established experimental feature is the two-fold rota-tional symmetry in the basal plane, for all three super-conductors. This two-fold symmetry has two incarna-tions, (1) in the electronic spin susceptibility, and (2) inthe momentum-dependence of the superconducting gapamplitude. As we explain below, these two aspects of thetwo-fold symmetry are independent of each other.The nematic combination of two components of E (1) u is the most well-known candidate for the two-fold sym-metry in the electronic spin susceptibility [15, 19, 22]. E (2 p ) u and E (3 p ) u can also lead to two-fold symmetry inthe spin susceptibility in the xy plane. Because for E (2 p ) u and E (3 p ) u the d ( k ) vectors in the pseudospin basis do nothave the third component, as shown in Eqs.(41) and (42),the corresponding spin susceptibility will keep invariantfor an out-of-plane magnetic field. On the other hand,the d ( k ) vectors for the two E (1) u pairing components,as shown in Eqs.(30) and (31), have finite third compo-nents. The spin susceptibility for a pairing in the E (1) u channel should therefore decrease slightly for an out-of-plane magnetic field [47]. In this respect, E (2 p ) u and E (3 p ) u fit the Knight shift experiment better than the E (1) u pair-ing [15].The two-fold symmetry in the k -dependence of the su-perconducting gap amplitudes could be explained by the nematic realization of E (1) u or E (3) g . E (2 m ) u and E (3 m ) u alsoresult in two-fold symmetry in the k x k y plane. But the( ka ) dependence of their pairing amplitude implies thatthey are unlikely the leading pairing instability. The re-maining channels in Table III and Table IV do not giveclear two-fold symmetry in the superconducting gap am-plitude in the k x k y plane to account for the experiments.Therefore, as regards the two-fold symmetry in the basalplane, the E (1) u channel is the only channel that naturallyaccount for both of the two aspects of the two-fold rota-tional symmetry. On the other hand, each incarnation ofthe two-fold symmetry has alternative realizations otherthan E (1) u .Now we consider further constraints imposed by thenodal structures of the bulk spectra. For both Cu x Bi Se [12, 67] and Sr x Bi Se [68], the bulk spectra are fullygapped according to relevant experiments. This meansthe absence of true or approximate nodes in the bulkspectrum. ˜ ψ (1)2 ( k ) of E (1) u is the only pairing componentin Table III to Table V that has approximate (point)nodes, for which the gap minima is two to three orders ofmagnitude smaller than the gap maxima [31]. The twocomponents of E (4 p ) u individually opens a full gap on theFermi surface if we have a corrugated cylindrical Fermisurface, with the size of the gap scaling as ka for wavevectors on the Fermi surface. For spheroidal Fermi sur-faces, both of the two components of E (4 p ) u give two pointnodes at k x = k y = 0. All the remaining individual pair-ing components in Table III to Table V (including ˜ ψ (1)1 ( k )of E (1) u ) give true nodes to the bulk spectrum, for bothspheroidal and corrugated cylindrical Fermi surfaces.For several of the pairing channels other than E (4 p ) u ,we can also get a fully-gapped bulk spectrum by form-ing chiral combinations of the two pairing components,for suitable Fermi surface topology. These channels in-clude E (1 , , g of the E g representation, and E (2 , , p ) u ofthe E u representation. When the Fermi surface is a cor-rugated cylinder, the chiral combinations (e.g., η = 1and η = i , or ˜ η = 1 and ˜ η = i ) of the two compo-nents of E (1 , , g or E (2 , , p ) u all lead to a fully-gappedbulk quasiparticle spectrum. For a chiral and nematiccombination (e.g., η = 0 . η = i , or ˜ η = 0 . η = i ) of the two components of E (1 , , g or E (2 , , p ) u , thebulk quasiparticle spectrum are also fully gapped, if wehave a corrugated cylindrical Fermi surface. When theFermi surface is a spheroid, these chiral or chiral and ne-matic pairings have two point nodes at the k x = k y = 0points of the Fermi surface. None of the above pairingchannels, including E (1 , , g and E (2 , , p, p ) u , can accountfor the two-fold symmetries observed in the experiments.The time-reversal symmetry breaking combinations ofthe two components of E (1) u have more complicated be-haviors. For a spheroidal Fermi surface, the purely chiralcombinations of the two components of E (1) u have twopoint nodes at the k x = k y = 0 points of the Fermi sur-1face [24, 26]. These point nodes disappear both whenwe consider a chiral and nematic combination of the twocomponents, and when the Fermi surface turns to a cor-rugated cylinder. On the other hand, all these varioustime-reversal symmetry breaking pairing combinationswithin the E (1) u channel support SABSs on the xy surfaceand keep the topological surface states at the Fermi levelungapped.Summing up the above discussions, for Cu x Bi Se andSr x Bi Se , there is not a single pairing channel thatcan simultaneously explain the two-fold symmetric andfully-gapped bulk spectrum, with a time-reversal sym-metric pairing combination. For Nd x Bi Se , while the E (1) u channel can explain the two-fold symmetry and bulkspectrum with (real or approximate) point nodes, thecomplexity of the Fermi surface of this compound com-plicates the comparison [43]. Focusing on Cu x Bi Se andSr x Bi Se , an implication of the above comparison is thatwe have to consider more than one pairing channels tosimultaneously account for the fully-gapped bulk spec-trum and the two-fold rotational symmetry in the basalplane. In other words, the pairing has to be multichan-nel , in addition to having two components. This is themajor conclusion of our comparison. What follows weexplore various possible multichannel pairing combina-tions. We relax the constraint of time-reversal symmetryand explore all the possible pairing combinations thatgive both the two-fold symmetry in the basal plane andthe fully-gapped bulk quasiparticle spectrum.To get a fully-gapped bulk spectrum with two-foldsymmetry in the superconducting gap amplitude, interms of a multichannel pairing in the E g representation,the ( η , η ) vector in Eq.(21) must be simultaneously chi-ral and nematic. The chirality and nematicity of the pair-ing are separately responsible for the fully-gapped bulkspectrum and the two-fold symmetry in the k x k y plane.The full gap may be achieved by a chiral combinationof E (1) g , or E (2) g , or E (7) g , under the premise that theFermi surface is a corrugated cylinder. Including a finitecontribution of E (3) g , the nematicity of the pairing canaccount for the two-fold symmetry in the gap amplitude.The fully-gapped bulk spectrum requires the strength ofthe E (1) g (or E (2) g , or E (7) g ) channel to be larger than thestrength of the E (3) g channel, namely ∆ (or ∆ , or ∆ )is sufficiently large compared to ∆ . In this case, the su-perconducting order parameter (the wave function of theCooper pairs) does not undergo a sign reversal upon scat-tering off the xy surface, implying the absence of SABSs.In addition, the chiral combination of the dominant E (1) g (or E (2) g , or E (7) g ) pairing component fully-gaps the topo-logical surface states. Therefore, there are no low-energyin-gap states on the xy surface for this chiral and nematicpairing in the E g representation.If we impose further the constraint of two-fold sym-metry in the electronic spin susceptibility [15], we haveto consider the E u channels. A combination of E (1) u (or E (2 p ) u , or E (3 p ) u ) and E (4 p ) u is possible to give a purelynematic pairing that has a fully-gapped bulk spectrum,if ˜∆ p is large compared to ˜∆ (or ˜∆ p , or ˜∆ p ). Letus define (I) the combination of E (1) u with E (4 p ) u , (II) thecombination of E (2 p ) u with E (4 p ) u , and (III) the combina-tion of E (3 p ) u with E (4 p ) u . For the Knight shift experiment,the combinations II and III explain the invariant c -axisspin susceptibility better than combination I. In addi-tion, they can also explain more naturally the absenceof SABSs and the gapped topological surface states onthe xy surfaces of the superconductor. For these experi-mental features, combinations II and III are better alter-natives to combination I. On the other hand, if we con-sider the two-fold symmetry in the superconducting gapamplitude, combination I explains it naturally, whereascombinations II and III do not lead to salient two-foldsymmetry in the gap amplitude. Overall, to explainat least qualitatively the key experimental features ofCu x Bi Se with a nematic and time-reversal symmetricpairing, we have to consider a pairing consisting mainly E (1) u and E (4 p ) u , and possibly supplemented by E (2 p ) u and(or) E (3 p ) u . Finally, following the same analysis as thatfor the above chiral and nematic E g pairing, there are nolow-energy in-gap states on the xy surface. Explicitly,because ˜∆ p is assumed to be large compared to ˜∆ (or˜∆ p , or ˜∆ p ), the superconducting order parameter doesnot undergo sign change upon the substitution of − k z for k z , so that there are no SABSs on the xy surface. Inaddition, since ˜∆ p gaps the topological surface states atthe Fermi level, there are no topological surface statesthat may contribute to the low-energy in-gap states.For other E u combinations, we again have to consider achiral and nematic combination to account for the fully-gapped bulk spectrum and the two-fold in-plane rota-tional symmetry at the same time, similar to the E g case. If both the superconducting gap amplitudes andthe electronic spin susceptibilities are required to be two-fold symmetric in the basal plane, we can choose the chi-ral and nematic combinations of E (1) u with E (2 , u . If onlythe electronic spin susceptibilities are required to be two-fold symmetric, the chiral and nematic combinations of E (2 p, p ) u with E (2 , u are eligible. Again, the pairing am-plitudes for E (2 , u , which give the fully-gapped bulk spec-trum, should be larger than the pairing amplitudes for E (1) u . Because the dominant third component of the d ( k )vector for E (2) u is proportional to B c z + R d ≃ B c z ,and both E (1) u and E (2) u not gap the topological surfacestates at the Fermi level, the chiral and nematic com-bination of E (1) u and E (2) u have SABSs and ungappedtopological surface states on the xy surface. The chi-ral and nematic combination of E (1) u and E (3) u , with themagnitude of ˜∆ much larger than the magnitude of ˜∆ to ensure a fully-gapped bulk spectrum, does not haveSABSs and ungapped topological surface states on the xy surface.A salient feature of the above multichannel pairings is2the ubiquity of the chiral and nematic pairing combina-tions, in both the E g and the E u representations. Besidesthe broken in-plane rotational symmetry, the chiral char-acter of the pairing implies that it can show typical sig-natures in experiments such as muon spin relaxation andoptical Kerr effect which probe the broken time-reversalsymmetry [70–72]. The chiral and nematic E u pairingsare also nonunitary, which can be verified by proving thenonvanishing of the vector product d ( k ) × d ∗ ( k ) = 0for general k [73, 74]. Nonunitary pairings, as a specialcase of pairings with broken time-reversal symmetry, areknown to have a spontaneous moment from the Cooperpairs [73, 74]. Also note that, the chiral or chiral andnematic pairing combinations of the two components of E (1) u are all nonunitary.Another common character of the above pairing com-binations is the requirement of a corrugated cylindricalFermi surface, to have a fully-gapped bulk spectrum. Theundoped Bi Se is known to have a spheroidal Fermi sur-face. As the concentration of the dopants increases, itis natural to expect a continuous evolution of the Fermisurface from spheroidal to corrugated cylindrical. Along-side, the nodal structure of the pairing may also change.While the experiment of Lahoud et al shows that thesuperconducting Cu x Bi Se has a corrugated cylindricalFermi surface in the normal phase [51], there is presentlyno systematic studies on the evolution of the Fermi sur-face for any of the three superconductors.Finally, it is interesting to notice that several pairingchannels in Table III to Table V exhibit prominent four-fold rotational symmetry, including the E (1 , , g channelsof the E g representation and the E (2 , u channels of the E u representation, through the k -dependence of the ϕ ( k )and ϕ ( k ) symmetry factors. The four-fold symmetryalso breaks the three-fold rotational symmetry of the un-derlying crystal lattice in the normal phase. Several ofthe chiral and nematic pairing combinations proposedabove contain these pairing channels. They are there-fore expected to exhibit some characters of the four-fold symmetry. On the experimental side, a four-foldsymmetric component in the superconducting gap am-plitude was indeed observed by Du et al for Sr x Bi Se [68]. These chiral and nematic pairings are however in-consistent with the general belief that the superconduct-ing state of Sr x Bi Se preserves the time-reversal sym-metry. V. IMPLICATIONS FOR FUTUREEXPERIMENTS
From the analysis of the previous section, the avail-able experiments have imposed stringent constraints onthe true pairing symmetries of the M x Bi Se (M is Cu,Sr, or Nd) superconductors. In particular, they pointto the multichannel two-component pairings. These ex-periments, on the other hand, do not provide sufficientand consistent information to unambiguously identify the true pairing symmetry of these superconductors. In par-ticular, for each of the three superconductors, the param-eter x defines a series of superconductors which mighthave qualitatively different Fermi surface topology in thenormal state and different pairing symmetries in the su-perconducting state. This doping dependence has notbeen systematically investigated for any of three super-conductors, although the phase diagram of Cu x Bi Se exists [75]. It is therefore highly desirable to make a sys-tematic experimental study on each member of the serieswith a complementary set of experimental tools.Based on the survey over the experimental consen-suses and the comparison with the possible pairing chan-nels, the relevant properties and the experiments thatmay be performed to probe them include: (1) The evo-lution of the Fermi surface with the doping concentra-tion x . The magnetic oscillation experiments, includingthe Shubnikov-de Haas oscillation and the de Haas-vanAlphen effect, can probe the Fermi surface contour inthe normal phase [43, 51, 76]. From the discussions ofthe above section, the geometry of the Fermi surface sig-nificantly influences the bulk quasiparticle spectrum ofthe superconducting state. It is highly desirable to deter-mine the systematic evolution of the Fermi surface as thedoping concentration x increases. In addition, the angle-resolved photoemission spectroscopy (ARPES), whichmay probe the continuous evolution of the topologicalsurface states with x [2, 51, 52], supplements the mag-netic oscillation experiment which gives the bulk Fermisurface. (2) The bulk quasiparticle spectrum of the su-perconducting state, fully gapped or not. For each dop-ing concentration that turns the material to a supercon-ductor, we may determine whether the superconductingstate is fully gapped or not by combining the scanningtunneling spectroscopy (STS) measurements [12] and themeasurements of the zero-field bulk specific heat [67].The STS experiment, besides probing the bulk quasi-particle spectrum, probes also the in-gap states on thesurface of the superconductors [12, 68]. The spin-latticerelaxation rates [77] and penetration depth [69] can alsobe used to probe the nodal structures of the bulk super-conducting spectrum. (3) The momentum-dependenceof the superconducting gap amplitude. This property onone hand refines the understanding obtained from theabove step, on the other hand shows the presence ornot of the anisotropy in the superconducting gap am-plitudes. The relevant experiments include the field-angle-dependent specific heat experiments [16, 78], thefield-angle-dependent upper critical magnetic field mea-surements [16, 18, 21], the field-angle-dependent resis-tivity measurements [17], and field-angle dependent STS[79]. Besides the above experiments for the magnitude ofthe superconducting gap, the phase of the superconduct-ing gap can be probed with the orientation-dependentJosephson junctions [80]. (4) The electronic spin suscep-tibility, which is the most direct diagnostic tool for thestructure of the d -vectors of the pseudospin-triplet pair-ings. The major relevant experiment is the Knight shift3measurement [15]. (5) The persistence or not of the time-reversal symmetry in the superconducting state. Thetime-reversal symmetry breaking of the superconductorscan be probed by several techniques, such as the muonspin resonance [70], optical Kerr effect [71, 72], and theJosephson effect [66]. The zero-field Hall effect was alsoused by Qiu et al to probe the broken time-reversal sym-metry in Nd x Bi Se [6]. (6) The presence or not of spon-taneous magnetization. The above chiral and nematic E u pairings are mostly nonunitary and should lead to spon-taneous magnetization. The experiment of Qiu et al re-porting the spontaneous magnetization of Nd x Bi Se wasbased on the field-dependent dc magnetization, which isa measurement of the global magnetization [6]. It is de-sirable to study the spontaneous magnetization with alocal measurement, like the polarized neutron diffractionexperiments [66].For each member of the three series of superconduc-tors, a systematic study of the above experiments shouldbe able to identify the genuine pairing symmetry. Fromthe breadth of the experimental tools involved, extensiveexperimental collaborations by sharing the same high-quality samples are highly desirable.Besides the experiments listed above, other experi-ments which may probe further implications of the can-didate pairings would also be of great potential inter-est. For example, as a natural consequence of the two-component nature of the pairing, the domain structureshould exist in the superconducting state. Althoughmost of the multichannel two-component pairings in-ferred from the analysis of the previous section do notsupport SABSs on the xy surface, all the E u pairingsshould give ABSs as domain wall states for domain wallsparallel to the z axis, because the superconducting orderparameters are odd functions of k x and k y in the pseu-dospin basis. Evidence of domains was reported in theexperiment of Yonezawa et al [16]. However, it is unclearwhether or not there are nontrivial domain wall states.Further experiments like those performed for Sr RuO [81] and superfluid He- A [82] are highly desirable. Asanother example, for nonunitary pairings in the E u rep- resentation, there should be collective modes associatedwith the oscillation of the magnetization of the Cooperpairs [73]. It would be interesting to detect this collec-tive mode via experiments like electron spin resonance orultrasound attenuation [73]. VI. SUMMARY
Starting from a tight-binding model for the normalstate electronic structures, we have constructed the fulllists of two-component pairings for M x Bi Se (M is Cu,Sr, or Nd). We then transform the pairings to the pseu-dospin basis, based on which we study their qualitativeproperties. Comparisons to the main experimental con-sensuses on these superconductors show that the truepairing symmetry for them has to be multichannel, inaddition to having two components. Besides a time-reversal symmetric nematic pairing belonging to the E u representation, we identify chiral and nematic pairings inboth the E g and the E u representations. However, for allthree superconductors, the existing experiments are in-sufficient to unambiguously determine the nature of thesuperconducting state. In particular, the studies on thedependence of the Fermi surface and the superconductingproperties on the doping concentration x are inadequatefor all three superconductors. A complementary set ofexperiments is suggested to identify unambiguously thegenuine pairing symmetries of the three superconductingelectron-doped Bi Se . Acknowledgments
L.H. thanks Ting-Kuo Lee, Sungkit Yip, and JunWang for helpful discussions. L.H. was supported byNSFC.11204035. C.S.T. was supported by the RobertA. Welch Foundation under Grant No. E-1146 and theTexas Center for Superconductivity at the University ofHouston. [1] 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. , 057001 (2010).[2] 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, Nat. Phys. , 855(2010).[3] Zhongheng Liu, Xiong Yao, Jifeng Shao, Ming Zuo, LiPi, Shun Tan, Changjin Zhang, and Yuheng Zhang, J.Am. Chem. Soc. , 10512 (2015).[4] Shruti, V. K. Maurya, P. Neha, P. Srivastava, and S.Patnaik, Phys. Rev. B , 020506(R) (2015).[5] C. Q. Han, H. Li, W. J. Chen, Fengfeng Zhu, Meng-YuYao, Z. J. Li, M. Wang, Bo F. Gao, D. D. Guan, Canhua Liu, C. L. Gao, Dong Qian, and Jin-Feng Jia, Appl. Phys.Lett. , 171602 (2015).[6] Y. Qiu, K. N. Sanders, J. Dai, J. E. Medvedeva, W. Wu,P. Ghaemi, T. Vojta, and Y. S. Hor, arXiv:1512.03519.[7] Satoshi Sasaki, M. Kriener, Kouji Segawa, Keiji Yada,Yukio Tanaka, Masatoshi Sato, and Yoichi Ando, Phys.Rev. Lett. , 217001 (2011).[8] T. Kirzhner, E. Lahoud, K. B. Chaska, Z. Salman, andA. Kanigel, Phys. Rev. B , 064517 (2012).[9] X. Chen, C. Huan, Y. S. Hor, C. A. R. S`a de Melo, Z.Jiang, arXiv:1210.6054.[10] T.V. Bay, T. Naka, Y. K. Huang, H. Luigjes, M. S.Golden, and A. de Visser, Phys. Rev. Lett. , 057001(2012). [11] Pradip Das, Yusuke Suzuki, Masashi Tachiki, and KazuoKadowaki, Phys. Rev. B , 220513(R) (2011).[12] N. Levy, T. Zhang, J. Ha, F. Sharifi, A. A. Talin, Y. Kuk,and J. A. Stroscio, Phys. Rev. Lett. , 117001 (2013).[13] Haibing Peng, Debtanu De, Bing Lv, Fengyan Wei, andChing-Wu Chu, Phys. Rev. B , 024515 (2013).[14] Shingo Yonezawa, arXiv:1812.01378.[15] K. Matano, M. Kriener, K. Segawa, Y. Ando, and Guo-qing Zheng, Nat. Phys. , 852 (2016).[16] Shingo Yonezawa, Kengo Tajiri, Suguru Nakata, YukiNagai, ZhiweiWang, Kouji Segawa, Yoichi Ando, andYoshiteru Maeno, Nat. Phys. , 123 (2017).[17] Guan Du, YuFeng Li, J. Schneeloch, R. D. Zhong, GenDaGu, Huan Yang, Hai Lin, and Hai-Hu Wen, Sci. China-Phys. Mech. Astron. , 037411 (2017).[18] Y. Pan, A. M. Nikitin, G. K. Araizi, Y. K. Huang, Y.Matsushita, T. Naka, and A. de Visser, Sci. Rep. , 28632(2016).[19] Liang Fu, Phys. Rev. B , 100509 (2014).[20] Tomoya Asaba, B. J. Lawson, Colin Tinsman, Lu Chen,Paul Corbae, Gang Li, Y. Qiu, Y. S. Hor, Liang Fu, andLu Li, Phys. Rev. X , 011009 (2017).[21] Junying Shen, Wen-Yu He, Noah Fan Qi Yuan, ZengleHuang, Chang-woo Cho, Seng Huat Lee, Yew San Hor,Kam Tuen Law, and Rolf Lortz, npj Quant. Mat. , 59(2017).[22] Liang Fu, Nat. Phys. , 822 (2016).[23] J¨orn W. F. Venderbos, Vladyslav Kozii, and Liang Fu,Phys. Rev. B , 094522 (2016).[24] J¨orn W. F. Venderbos, Vladyslav Kozii, and Liang Fu,Phys. Rev. B , 180504(R) (2016).[25] Yuki Nagai and Yukihiro Ota, Phys. Rev. B , 134516(2016).[26] Noah F. Q. Yuan, Wen-Yu He, and K. T. Law, Phys.Rev. B , 201109(R) (2017).[27] Luca Chirolli, Fernando de Juan, and Francisco Guinea,Phys. Rev. B , 201110(R) (2017).[28] Fengcheng Wu and Ivar Martin, Phys. Rev. B , 144504(2017).[29] Fengcheng Wu and Ivar Martin, Phys. Rev. B , 224503(2017).[30] A. A. Zyuzin, Julien Garaud, and Egor Babaev, Phys.Rev. Lett. , 167001 (2017).[31] Lei Hao and C. S. Ting, Phys. Rev. B , 144512 (2017).[32] Liang Fu and Erez Berg, Phys. Rev. Lett. , 097001(2010).[33] Xiangang Wan, Sergey Y. Savrasov, Nature Commun. ,4144 (2014).[34] P. M. R. Brydon, S. Das Sarma, Hoi-Yin Hui, and JayD. Sau, Phys. Rev. B , 184512 (2014).[35] Xiao-Long Zhang and Wu-ming Liu, Sci. Rep. , 8964(2015).[36] Lei Hao, Gui-Ling Wang, Ting-Kuo Lee, Jun Wang,Wei-Feng Tsai, and Yong-Hong Yang, Phys. Rev. B ,214505 (2014).[37] Chao-Xing Liu, Xiao-Liang Qi, Haijun Zhang, Xi Dai,Zhong Fang, and Shou-Cheng Zhang, Phys. Rev. B ,045122 (2010).[38] Liang Fu, Phys. Rev. Lett. , 266801 (2009).[39] Haijun Zhang, Chao-Xing Liu, Xiao-Liang Qi, Xi Dai,Zhong Fang, and Shou-Cheng Zhang, Nature Phys. ,438 (2009).[40] Lei Hao and T. K. Lee, Phys. Rev. B , 134516 (2011).[41] Qiang-Hua Wang, Da Wang, and Fu-Chun Zhang, Phys. Rev. B , 035104 (2010).[42] P. Larson, V. A. Greanya, W. C. Tonjes, Rong Liu, S. D.Mahanti, C. G. Olson, Phys. Rev. B , 085108 (2002).[43] B. J. Lawson, Paul Corbae, Gang Li, Fan Yu, TomoyaAsaba, Colin Tinsman, Y. Qiu, J. E. Medvedeva, Y. S.Hor, and Lu Li, Phys. Rev. B , 041114(R) (2016).[44] S.-K. Yip, Phys. Rev. B , 104505 (2013).[45] Sungkit Yip, arXiv:1609.04152.[46] Bj¨orn Zocher and Bernd Rosenow, Phys. Rev. B ,155138 (2013).[47] Tatsuki Hashimoto, Keiji Yada, Ai Yamakage, MasatoshiSato, and Yukio Tanaka, J. Phys. Soc. Jpn. , 044704(2013).[48] Yuki Nagai, Hiroki Nakamura, and Masahiko Machida,J. Phys. Soc. Jpn. , 053705 (2014).[49] Shota Takami, Keiji Yada, Ai Yamakage, Masatoshi Sato,and Yukio Tanaka, J. Phys. Soc. Jpn. , 064705 (2014).[50] Lei Hao and Ting-Kuo Lee, J. Phys.: Condens. Matter , 105701 (2015); ibid arXiv:1407.3329v2.[51] E. Lahoud, E. Maniv, M. Shaviv Petrushevsky, M.Naamneh, A. Ribak, S. Wiedmann, L. Petaccia, Z.Salman, K. B. Chashka, Y. Dagan, and A. Kanigel, Phys.Rev. B , 195107 (2013).[52] Y. Tanaka, K. Nakayama, S. Souma, T. Sato, N. Xu,P. Zhang, P. Richard, H. Ding, Y. Suzuki, P. Das, K.Kadowaki, and T. Takahashi, Phys. Rev. B , 125111(2012).[53] M. S. Dresselhaus, G. Dresselhaus, and A. Jorio, GroupTheory-Application to the Physics of Condensed Matter (Springer-Verlag, Berlin, 2008).[54] E. I. Blount, Phys. Rev. B , 2935 (1985).[55] Sungkit Yip, Annu. Rev. Cond. Matter Physics , 15(2014).[56] In a previous work [36], because of not incorporating theextra factor of σ ⊗ ( − is ), E (3) g was mistakenly catego-rized as belonging to the E u representation.[57] Sungkit Yip and Anupam Garg, Phys. Rev. B , 3304(1993).[58] Yuxuan Wang and Liang Fu, Phys. Rev. Lett. ,187003 (2017).[59] Bruno Uchoa and A. H. Castro Neto, Phys. Rev. Lett. , 146801 (2007).[60] Carsten Honerkamp, Phys. Rev. Lett. , 146404(2008).[61] Hong-Yan Lu, Lei Hao, Rui Wang, and C. S. Ting, Phys.Rev. B , 241410(R) (2016).[62] Lei Hao and Jun Wang, J. Phys.: Condens. Matter ,255701 (2015).[63] Anthony J. Leggett, Rev. Mod. Phys. , 331 (1975).[64] Robert Joynt and Louis Taillefer, Rev. Mod. Phys. ,235 (2002).[65] Andrew Peter Mackenzie and Yoshiteru Maeno, Rev.Mod. Phys. , 657 (2003).[66] Yoshiteru Maeno, Shunichiro Kittaka, Takuji Nomura,Shingo Yonezawa, and Kenji Ishida, J. Phys. Soc. Jpn. , 011009 (2012).[67] M. Kriener, Kouji Segawa, Zhi Ren, Satoshi Sasaki, andYoichi Ando, Phys. Rev. Lett. , 127004 (2011).[68] Guan Du, Jifeng Shao, Xiong Yang, Zengyi Du, De-long Fang, Jinghui Wang, Kejing Ran, Jinsheng Wen,Changjin Zhang, Huan Yang, Yuheng Zhang, and Hai-Hu Wen, Nature Commun. , 14466 (2017).[69] M. P. Smylie, H. Claus, U. Welp, W.-K. Kwok, Y. Qiu, Y. S. Hor, and A. Snezhko, Phys. Rev. B , 180510(R)(2016).[70] G. M. Luke, Y. Fudamoto, K. M. Kojima, M. I. Larkin, J.Merrin, B. Nachumi, Y. J. Uemura, Y. Maeno, Z. Q. Mao,Y. Mori, H. Nakamura, and M. Sigrist, Nature (London) , 558 (1998).[71] S. Spielman, J. S. Dodge, L. W. Lombardo, C. B. Eom,M. M. Fejer, T. H. Geballe, and A. Kapitulnik, Phys.Rev. Lett. , 3472 (1992).[72] Jing Xia, Yoshiteru Maeno, Peter T. Beyersdorf, M.M. Fejer, and Aharon Kapitulnik, Phys. Rev. Lett. ,167002 (2006).[73] Tetsuo Ohmi and Kazushige Machida, J. Phys. Soc. Jpn. , 4018 (1996).[74] A. D. Hillier, J. Quintanilla, B. Mazidian, J. F. Annett,and R. Cywinski, Phys. Rev. Lett. , 097001 (2012).[75] M. Kriener, Kouji Segawa, Zhi Ren, Satoshi Sasaki,Shohei Wada, Susumu Kuwabata, and Yoichi Ando,Phys. Rev. B , 054513 (2011). [76] Ben J. Lawson, Y. S. Hor, and Lu Li, Phys. Rev. Lett. , 226406 (2012).[77] Yuki Nagai, Yukihiro Ota, and Masahiko Machida, Phys.Rev. B , 180502(R) (2015); Yuki Nagai and YukihiroOta, Phys. Rev. B , 134516 (2016).[78] Kristin Willa, Roland Willa, Kok Wee Song, G. D. Gu,John A. Schneeloch, Ruidan Zhong, Alexei E. Koshelev,Wai-Kwong Kwok, and Ulrich Welp, Phys. Rev. B ,184509 (2018).[79] Ran Tao, Ya-Jun Yan, Xi Liu, Zhi-Wei Wang, YoichiAndo, Qiang-Hua Wang, Tong Zhang, and Dong-LaiFeng, Phys. Rev. X , 041024 (2018).[80] Yukio Tanaka and Satoshi Kashiwaya, Phys. Rev. B ,892 (1997).[81] Francoise Kidwingira, J. D. Strand, D. J. Van Harlingen,and Yoshiteru Maeno, Science , 1267 (2006).[82] J. Kasai, Y. Okamoto, K. Nishioka, T. Takagi, and Y.Sasaki, Phys. Rev. Lett.120