Magnetic-impurity resonance states for different pairing symmetries in twisted bilayer graphene
MMagnetic-impurity resonance states for different pairing symmetries in twisted bilayergraphene
Liang Chen, ∗ Hui-Zhen Li, and Rong-Sheng Han School of Mathematics and Physics, North China Electric Power University, Beijing, 102206, China (Dated: September 5, 2018)In this work, we study the magnetic-impurity resonance states in superconducting phase of twistedbilayer graphene for different pairing symmetries. Using the two-orbital model proposed by Yuanand Fu in [
Phys. Rev. B , 045103 (2018)], we find that when the impurity is located at one siteof the emergent honeycomb lattice, the spacial distributions of the resonance states will break boththe threefold and twofold rotation symmetries of D group for pairing symmetries which belong tothe irreducible representations of this point group. When the magnetic impurity is located at thecenter of the emergent honeycomb lattice, the appearance of resonance peak at the position closeto the impurity can be considered as a strong evidence of non- s -wave pairing. PACS numbers: 74.20.Rp, 74.25.Ha, 74.78.Fk
I. INTRODUCTION
Recent experiments [1, 2] on bilayer graphene with asmall twisted angle have attracted great attention. The-oretical investigations [3–7] show that when the twistedangle is close to the ‘magic’ angles, the bilayer grapheneexhibits some low energy flat bands in the mini Brillouinzone (MBZ) of the moir´e patten superlattice. Exper-imental measurement [2] demonstrates that these flatbands display the insulating phase at half-filling, i.e.,the charge density is tuned to ± e per unit cell of thesuperlattice via gate voltage. When the charge den-sity is doped slightly away from these correlated insu-lating phases, superconducting phases are observed. Fora record-low carrier density of a few 10 cm − , a remark-able high critical temperature T c = 1 .
7K is observed [1].The strong correlation features in this material are verysimilar to the high- T c cuprate superconductors. Whichmakes graphene moir´e superlattice to be an experimen-tal tunable platform for unconventional superconductiv-ity with completely carbon atoms.The nature of the superconductivity phases and themechanism for the relative high critical temperature in‘magic angle’ twisted bilayer graphene (TwBG) are notclear so far. Theoretical investigations from strong cor-relation physics show that the unconvertional pairingsare preferred in this material, i.e., the topological non-trivial d + id pairing is repeatedly proposed in theoret-ical works[8–13], the p -wave pairing [14], s ± and s ++ pairing [15] are also candidates dependent on differentconditions. Some other theoretical investigations showthat the TwBG exhibits strong electron-phonon couplingat certain twisted angles and electron densities, whichleads to a remarkable high critical temperature of mag-nitude 1 ∼
10K [16–18]. As demonstrated in these works,the pairing potential mediated by electron-phonon cou-pling reveals the s -wave symmetry. Further investiga- ∗ Corresponding Email: [email protected] tions show that, if the electron-electron interaction issufficient large, the d -wave symmetric pairing potentialis also possible for the electron-phonon coupling mecha-nism [16]. Among the comprehensive investigations onthe nature of superconductivity in this material, the de-termination of pairing symmetry is of fundamental im-portance.In this work, we study the magnetic impurity reso-nance states for typical pairing symmetries in the super-conductivity phases of TwBG. Like the standard phase-sensitive tetracrystal measurements [19, 20] and quasi-particle interference (QPI) experiments [21–24], the lo-cal density of states (LDOS) of the resonance states neara magnetic (and nonmagnetic) impurity is an impor-tant method to uncover the symmetries of pairing poten-tials in unconventional superconductors. This method iswidely used in the high- T c cuprate superconductors [25–34], iron-based superconductors [35–37], chiral p -wave su-perconductors and topological superconductors [38–41],etc.The paper is organized as follows: In Sec. II, wepropose the model Hamiltonian with four different pair-ing symmetries and give a theoretical derivation of theLDOS. In Sec. III, we analysis the numerical results forthe pairing symmetries. According to different represen-tations of the D group, the LDOS for three kinds ofdifferent hybridizations are presented. A conclusion isgiven in Sec. IV. II. THEORETICAL MODEL
We use the two-orbital tight-binding model proposedby Yuan and Fu [42] with a few modifications to describethe low-energy physics of TwBG, which is expressed in a r X i v : . [ c ond - m a t . s t r- e l ] S e p the following form, H = (cid:88) k ψ † k h ( k ) ψ k , (1) h ( k ) = (cid:88) i,j (cid:15) i,j ( k ) τ i ⊗ χ j , (2)where i, j = 0 , x, y, z , τ and χ are the 2 × { p x , p y } orbital space, respectively. τ x,y,z and χ x,y,z are the corresponding Pauli matrices in these spaces. ψ k = ( c AB ,p x , k , c AB ,p y , k , c BA ,p x , k , c BA ,p y , k ) T , c α,o, k is theannihilation operator of electron state for specified sub-lattice α = AB , BA, orbital o = p x , p y , and wave-vector k (the spin index is suppressed here). The superscript T means matrix transpose. As the Hermitian conjugate of ψ k , ψ † k is the creation operator of corresponding electronstate. The summation in Eq. (1) is taken over the MBZof TwBG. Fig. 1(a) shows the primitive vectors, unitcell, and p x,y -orbitals of the emergent superlattice. Theexplicit expressions of (cid:15) i,j ( k ) are given in the followingforms, (cid:15) x, ( k ) = t (cid:18) cos k x √ k x √ k y (cid:19) , (3) (cid:15) y, ( k ) = t (cid:18) sin k x √ − k x √ k y (cid:19) , (4) (cid:15) , ( k ) = (2 t + t (cid:48) ) (cid:32) √ k x k y k y (cid:33) − µ, (5) (cid:15) x,z ( k ) = t (cid:48) (cid:18) cos k x √ − cos k x √ k y (cid:19) , (6) (cid:15) y,z ( k ) = t (cid:48) (cid:18) sin k x √ k x √ k y (cid:19) , (7) (cid:15) x,x ( k ) = −√ t (cid:48) sin k x √ k y , (8) (cid:15) y,x ( k ) = −√ t (cid:48) cos k x √ k y , (9) (cid:15) ,z ( k ) = t (cid:48) (cid:32) cos √ k x k y − cos k y (cid:33) , (10) (cid:15) ,x ( k ) = −√ t (cid:48) sin √ k x k y , (11) (cid:15) ,y ( k ) = − t (cid:48)(cid:48) sin √ k x √ k x − k y √ k x + 3 k y , (12)where µ is the Fermi energy, t is the real hopping ampli-tude between nearest neighbors within different sublat-tices, t denotes the real hopping amplitude between nextnearest neighbors within the same sublattice. t (cid:48)(cid:48) rep-resents the inter-sublattice hopping amplitude betweenfifth nearest neighbors with definite chirality. This termbreaks the emergent SU (4) symmetry and hence the four-fold degeneracy (orbital and spin) along the ΓM line in FIG. 1. (color online) (a) Schematic representation of thetwo-orbital tight-binding model. The two primitive vectorsare chosen as A , = A ( √ / , ± / A is the super-lattice constant, which is set to be unit length in this work.The red dots and blue circles represent the two sublatticescorresponding to the AB and BA spots of the moir´e patten.On each site of the emergent honeycomb lattice, there is aset of p x - and p y -like Wannier orbitals which belongs to the E representation of the D point group (see Tab I for the irre-ducible representations of D group). The original point is setto be located at the impurity site. (b) The nearest neighborpairing. For the extended s-wave pairing, the potential is setto be ∆ = ∆ = ∆ = ∆ s (cid:48) . For the p -wave pairing, the po-tential is set to be ∆ = ∆ p , ∆ = e iφ ∆ p , ∆ = e iφ ∆ p with φ = 2 π/
3. (c) The next nearest neighbor pairing. For the d + id -wave pairing, the potential is set as ∆ = ∆ = ∆ d + id ,∆ = ∆ = e iφ ∆ d + id , ∆ = ∆ = e iφ ∆ d + id . the MBZ is split. The terms proportional to t (cid:48) and t (cid:48) refer to the nearest neighbor and next to nearest neigh-bor hopping amplitudes when the inter-orbital scatteringis included. These terms further break the orbital U (1)symmetry and the fourfold degeneracy along the KΓ andMK lines in band structure. It is easy to check thatthis Hamiltonian preserves all the symmetries proposedin Ref. [42]. Tab. I shows the group elements and irre-ducible representations of the point group D . For thecurrent Hamiltonian, the specified representations of thegroup elements are present in the caption of Tab. I. Fig.2 shows the band structure and density of states (DOS)of the model Hamiltonian with proper parameters. Thethree dashed lines in Fig. 2(b) from top to bottom de-note the typical charge density +2 e , 0 and − e in unitcell of the superlattice, respectively. One can find thatfilling ratios 1 / − / ± /
2) and thenear-by Van Hove singularities for the given parameters.When the relative strong electron-electron interactionsare presented, TwBG can be driven into the Mott insu-late phase at these two half-filling ratios. The experiment[1] shows that the superconducting phase appears whenthe electron concentration is doped slightly away fromthese half-filling ratios.Generally, the pairing potential term can be writtenas, H pairing = (cid:88) k ,s,s (cid:48) ψ T − k ,s ∆ s,s (cid:48) ( k ) ψ k ,s (cid:48) + h.c., (13) D E C C (cid:48) linear quadratic cubic A +1 +1 +1 - x + y x ( x − y ) A +1 +1 − y (3 x − y ) E +2 − x, y ) ( x − y , xy ) ( x + xy , x y + y )TABLE I. Character tab for point group D . A , A and E in the first column represent the three irreducible represen-tations. E in the first row refers to the identity element ofthe D point group. 2 C means the two group elements re-lated to threefold rotation operations along the z -axis, underthe basis shown in the context below Eq. (2), these groupelements are given by, C z = τ ⊗ ( − χ − i √ χ y ) and thesquare of C z , which is also the inverse of C z , denoted as C z .3 C (cid:48) in the first row means the three group elements relatedto the twofold rotation operations along the three symmetricaxes, under the basis shown in the context below Eq. (2),these group elements are given by C (cid:48) y = − τ x ⊗ χ z , C z C (cid:48) y C − z ,and C z C (cid:48) y C − z . The last three columns show the exemplifiedpolynomial representations.FIG. 2. (color online) (a) Band structure of the tight-bindingmodel Hamiltonian (1) for t = 1, t = t / t (cid:48) = t / t (cid:48) = t / t (cid:48)(cid:48) = 0 and µ = 0. (b) The corresponding DOS.The dashed lines show the filling ratio, the specified fillingratios 1 / − / where s, s (cid:48) = ↑ , ↓ refer to the spins of electron states, h.c. means Hermitian conjugate, ∆ s,s (cid:48) ( k ) is the pairing po-tential matrix. In this work, we consider only the spin-singlet pairing (potential spin-triplet p -wave and d + id pairing are discussed in Refs. [13, 43]), such that thesummation over spin indices contains only one situation, s = ↑ and s (cid:48) = ↓ . Without causing confusion, we suppressthe superscript of pairing potential. In this case, ∆( k )is a 4 × s -wave pairing,∆ s ( k ) = ∆ s τ ⊗ χ , (14)the second one is the nearest neighbor (extended) s -wave FIG. 3. (color online) Fermi surface for different chemical po-tentials. (a) µ = − . − / µ = 1 . / µ = − . µ = 1 . pairing,∆ s (cid:48) ( k ) = ∆ s (cid:48) (cid:20)(cid:18) cos k x √ k x √ k y (cid:19) τ x ⊗ χ + (cid:18) sin k x √ − k y k x √ (cid:19) τ y ⊗ χ (cid:21) , (15)the third one is the nearest neighbor p -wave pairing,∆ p ( k ) = ∆ p (cid:20)(cid:18) cos k x √ − cos k x √ k y i √ k x √ k y (cid:19) τ x ⊗ χ + (cid:18) sin k x √
3+ sin k x √ k y i √ k x √ k y (cid:19) τ y ⊗ χ (cid:21) , (16)the last one is the d + id pairing between next to nearestneighbor sites,∆ d + id ( k ) = 2∆ d + id (cid:32) cos k y − cos k y √ k x i √ k y √ k x (cid:33) τ ⊗ χ , (17)where the constants ∆ s , ∆ s (cid:48) , ∆ p and ∆ d + id representthe strength of the pairing potentials for these four dif-ferent symmetries. It is easy to check that the s -waveand extended s -wave pairing potentials belong to theone-dimensional A representation of the D group, the p -wave and d + id pairing potentials belong to the two-dimensional representation (the E representation shownin Tab. I) of the D group. Figs. 1(b) and 1(c) showthe nearest neighbor and the next to nearest neighborpairing potentials in superlattice space.The total Hamiltonian for a magnetic impurity coupledto the superconducting TwBG is consisted of four terms, H = H imp + H hyb + H + H pairing , (18) H imp = (cid:88) s (cid:15) d d † s d s + U d †↑ d ↑ d †↓ d ↓ , (19) H hyb = 1 √ N (cid:88) k ,s (cid:16) ψ † k ,s ¯ V k d s + h.c. (cid:17) , (20)where (cid:15) d is the impurity energy, U is the on-site Coulombinteraction, d † s and d s are the creation and annihilationoperators of spin- s impurity state, respectively. N is thetotal number of wave-vectors in the MBZ. ¯ V k is the 4 × k . This term has to beconsidered very carefully. Generally, the impurity Hamil-tonian (19) belongs to the A representation of the D group, and the two Wannier orbitals belong to the two-dimensional representation of the D group. If the im-purity is coupled to only one site of the emergent hon-eycomb lattice, i.e., shown by Eq. (23) in the followingcontext, the hybridization term (20) will break both thethreefold and twofold rotation symmetries. Firstly, weconsider the C rotation symmetric hybridizations, weconsider that the impurity is located at the center of theemergent honeycomb lattice, i.e., the center of an AAspot of the moir´e pattern. The impurity is hybridized tothe six nearest neighbors symmetrically, detailed analysisshow that there are only two different hybridizations be-long to the two-dimensional representation of D group,Fig. 4 shows the form of these hybridizations in realspace. Their explicit expressions in wave-vector spaceare given by,¯ V (1) k = V (cid:16) e i kx √ − e − i kx √ cos k y (cid:17) √ ie − i kx √ sin k y − (cid:16) e − i kx √ − e i kx √ cos k y (cid:17) √ ie i kx √ sin k y , (21)¯ V (2) k = V √ ie − i kx √ sin k y − (cid:16) e i kx √ − e − i kx √ cos k y (cid:17) − √ ie i kx √ sin k y − (cid:16) e − i kx √ − e i kx √ cos k y (cid:17) , (22)where V represents the strength of the hybridizations.When the impurity is located at the center of an AB spotand hybridized to the p x and p y orbitals at this emergentsite equally, the hybridization matrix in Eq. (20) can be FIG. 4. (color online) The two different hybridizations belongto the two-dimensional representation of the D group. Themagnetic impurity (labeled with black dot) is located at thecenter of the emergent honeycomb lattice. The orientationsof the orbital labels refer to proper combinations of the p x,y -orbitals at specified sites. (a) The real space representationof the 1st hybridization ¯ V (1) k given in Eq. (21). (b) The realspace representation of the 2nd hybridization ¯ V (2) k given inEq. (22). written as, ¯ V (3) k = ( V , V , , T . (23)We need to emphasize that, in this work, the originalpoint is set to be located at the impurity site, whichmeans that, for the hybridizations ¯ V (1) k and ¯ V (2) k , theoriginal point is located at the center of the emergenthoneycomb lattice, for the hybridization ¯ V (3) k , the origi-nal point is located at the center of one of the AB spotof the moir´e pattern.Now we consider the solution of the total Hamil-tonian (18), in the strong Coulomb interaction limit, U → ∞ , the double occupation state on the impu-rity site can be excluded, this condition may be repre-sented by introducing the slave-boson operators b and b † , d s = f s b † , d † s = bf † s . The extra degrees of freedomcan be ruled out by the no-double occupation condition, Q = b † b + (cid:80) s f † s f s = I . In the mean-field approxima-tion, the slave-boson operators can be replaced by theexpectation value (cid:104) b (cid:105) = (cid:104) b † (cid:105) = b , the constraint condi-tion can be approximated by introducing a Lagrangianmultiplier term to the Hamiltonian, λ ( b + (cid:80) s f † s f s − b and λ can be determined by minimiz-ing the free energy of the mean-field Hamiltonian. Thismean-field Hamiltonian is given by, H MF = H MFimp + H MFhyb + H + H pairing + λ ( b − , (24) H MFimp = (cid:88) s ˜ (cid:15) d f † s f s , (25) H MFhyb = 1 √ N (cid:88) k ,s ψ † k ,s ˜¯ V k d s + h.c., (26)where ˜ (cid:15) d = (cid:15) d + λ is the renormalized impurity energy,˜¯ V k = b ¯ V k is the renormalized hybridization matrix. Inthe Bogoliubov-de Gennes (BdG) formalism, the mean-field Hamiltonian can be recast as, H BdG = λ ( b −
1) + ˜ (cid:15) d + Φ † ΛΦ + (cid:88) k Ψ † k h BdG ( k )Ψ k + 1 √ N (cid:88) k (cid:104) Ψ † k V k Φ + h.c. (cid:105) (27)where the Nambu spinor is defined as, Ψ † k =( ψ † k , ↑ , ψ T − k , ↓ ), Φ = ( f ↑ , f †↓ ) T , Λ = ˜ (cid:15) d τ ⊗ χ ⊗ ς z , h BdG ( k ) = (cid:18) h ( k ) ∆( k )∆ † ( k ) − h ( k ) (cid:19) , V k = (cid:32) ˜¯ V k − ˜¯ V k (cid:33) , (28) ς z is the third Pauli matrix in the Nambu spinor space.Using the standard functional integration method inquantum field theory, we find that the finite temperaturefree energy of the magnetic impurity can be written as, F = λ ( b −
1) + ˜ (cid:15) d + k B T (cid:88) n Tr ln[ G f ( iω n )] , (29)where k B is the Boltzmann constant and T refers totemperature, ω n = (2 n + 1) πk B T is the Matsubara fre-quency of fermion, G f ( iω n ) = [ iω n − Λ − Σ f ( iω )] − is the Green’s function of impurity state expressed inimaginary frequency representation, Σ f ( iω n ) is the self-energy of the impurity state. It is given by, Σ f ( iω n ) = N (cid:80) k V † k G (0) ψ ( iω n , k ) V k , where G (0) ψ ( iω n , k ) = [ iω n − h BdG ( k )] − is the unperturbed Green’s function of theconduction states in superconducting TwBG. The pa-rameters λ and b are determined by minimizing thefree energy (29), which gives, b + k B T (cid:88) n Tr [ G f ( iω n ) ς z ] = 0 , (30) λ b + k B T (cid:88) n Tr [ G f ( iω n )Σ f ( iω n )] = 0 , (31)The LDOS near the magnetic impurity is obtained byanalytic continuation of the imaginary-time Green’s func-tion, iω n → E + i + , ρ ψ ( E, R ) = − π ImTr (cid:20) G ψ ( E ; R , R ) 1 + ς z ⊗ τ ⊗ χ (cid:21) , (32)where G ψ ( E ; R , R ) is the full Green’s function of theconduction states, which is given by [28], G ψ ( E ; R , R ) = 1 N (cid:88) k , k (cid:48) e i ( k − k (cid:48) ) · R G ψ ( E ; k , k (cid:48) ) , (33) G ψ ( E ; k , k (cid:48) ) = G (0) ψ ( E, k ) (cid:104) δ k , k (cid:48) + T k , k (cid:48) ( E ) G (0) ψ ( E, k (cid:48) ) (cid:105) , (34)where T k , k (cid:48) ( E ) = V k G f ( E ) V † k (cid:48) /N is the T -matrix. III. NUMERICAL RESULTS
Before the detailed discussion of the numerical results,we need to clarify that, as demonstrated in previousworks [28, 44–49], there are two distinct phases for mag-netic doping in superconductor and marginal fermi liq-uid. When the hybridization between the impurity andhost material is not strong enough, the impurity be-haviors as an isolated impurity and decoupled to thehost material, as a consequence, Eqs. (30)-(31) haveno solution. For the other situation, the impurity andthe host material are strong correlated, the in-gap Yu-Shiba-Rusinov [50–52] resonance states present. Here-after, the hybridization and impurity energy are chosenin the strong correlation regime such that we can get thein-gap resonance states. In the practical calculation, amesh size of 3 × × . × − .Now we consider the LDOS for charge density equalsto − e in a unit cell of the moir´e pattern superlattice, µ = − .
7. Figs. 5(a)-5(d) show the LDOSs as func-tions of energy at some specified positions in real spacefor hybridization ¯ V (1) k and different symmetries. In thiscase, the impurity (the original point) is located at thecenter of the emergent honeycomb lattice. The red lineslabeled with ‘ × ’ and solid blue lines present the LDOSsfor R = (0 ,
0) and R = ( − / √ , + ’ refer to the clean LDOSs withoutimpurity. For the s -wave pairing (Fig. 5(a)) and ex-tended s -wave pairing (Fig. 5(b)), one can find that thein-gap resonance peak at R = (0 ,
0) is vanishing, we needto emphasize that this is not because V is too small tocouple the impurity and the superconductor, actually,for the given parameters V = ∆ s and (cid:15) d = − ∆ s / b and λ . Furthermore, detailed numerical analysis showthat the red line labeled with ‘ × ’ and the gray line la-beled with ‘+’ in Fig. 5(a) (and also Fig. 5(b)) co-incide to each other exactly. These observations leadto a central conclusion in this work. Here we demon-strate from symmetry analysis that the non-vanishingresonance peak at the center of the AA spot (the center ofthe emergent honeycomb lattice) can be considered as astrong evidence of the unconventional pairings belong tothe two-dimensional representation or the A represen-tation of the D group . Generally, the in-gap resonancepeak is induced by the impurity scattering term in Eq.(34). If the pairing potential belongs to the A repre-sentation (the identity representation), the BdG Hamil-tonian (27) will be invariant under the D group oper-ations, i.e., C H BdG C − = H BdG and C (cid:48) H BdG C (cid:48)− = H BdG , so that C G (0) ψ ( E, k ) C − = G (0) ψ ( E, C − k C ), C (cid:48) G (0) ψ ( E, k ) C (cid:48)− = G (0) ψ ( E, C (cid:48)− k C (cid:48) ). Here C and C (cid:48) denote the threefold and twofold rotation opera-tions given in Tab I. The hybridizations ¯ V (1) k and ¯ V (2) k FIG. 5. (color online) LDOS for the case with the first hybridization ¯ V (1) k , the impurity is located at the center of the emergenthoneycomb lattice, R = (0 , s -wave, extended s -wave, p -wave and d + id pairing,respectively. The gray lines labeled with ‘+’ in (a)-(d) show the clean LDOS without impurity. The red lines labeled with ‘ × ’in (a)-(d) show the LDOS near the impurity R = (0 , R = ( − / √ , E <
0) and right(
E >
0) in-gap resonance peaks of the blue line in (b), respectively (extended s -wave pairing). (g) and (i) show the real spacedistributions of the LDOS corresponding to the two in-gap resonance peaks of the blue line in (c), respectively ( p -wave pairing).(h) and (j) present the real space distributions of the LDOS of the two in-gap resonance peaks for d + id pairing shown in (d).The silvery lines in (e)-(j) show the emergent honeycomb lattice. The parameters are chosen as follows: ∆ s = ∆ s (cid:48) = ∆ p = 0 . d + id = 0 . V = ∆ s , (cid:15) d = − ∆ s / preserve the desired symmetries, C ¯ V (1 , k = ¯ V (1 , C − k C , C (cid:48) ¯ V (1 , k = ¯ V (1 , C (cid:48)− k C (cid:48) , so we get, C (cid:88) k G (0) ψ ( E, k ) V k C − = (cid:88) k G (0) ψ ( E, k ) V k , (35) C (cid:48) (cid:88) k G (0) ψ ( E, k ) V k C (cid:48)− = (cid:88) k G (0) ψ ( E, k ) V k , (36)it is easy to check that the only solution of these equa-tions is (cid:80) k G (0) ψ ( E, k ) V k = 0. Substituting this re-sult into Eqs. (34) and (33), one can find that, at theoriginal point, R = (0 , s -wavepairing and extended s -wave pairing which belong tothe A representation. For the topological nontrivial p -wave pairing and d + id pairing, one can check that C z ∆ p ( k ) C − z = e − πi/ ∆ p ( C − z k C z ), C z ∆ d + id ( k ) C − z = e πi/ ∆ d + id ( C − z k C z ), such that the in-gap resonancepeaks may be non-vanishing, as evident in Figs. 5(c)and 5(d), there does exhibit one resonance peak belowthe fermi energy for each case (see the red lines labeledwith ‘ × ’).Previous theoretical investigations show that, both thequantum fluxes [53] and magnetic impurities [54] can in-duce Majorana bound states in topological superconduc-tors. We need to emphasize that the unique in-gap reso-nance peak can not be regarded as the Majorana boundstate because the resonance peak is not located at E = 0.In addition, the numerical results show that the other res- onance peak with opposite energy appears at R (cid:54) = (0 , R = ( − / √ , R = (0 ,
0) does notmean that the resonance state is absent. The blue line inFig. 5(b) shows LDOS versus energy at R = ( − / √ , p -wave pairing potential and the d + id pairing potential,respectively. One can find that all of these patterns revealthe sixfold rotation symmetry. The maximum intensityis located at R = (0 ,
0) in Figs. 5(g) and 5(h), whichshow the spacial distributions of LDOSs for the negativeresonance energies. For the other cases, the maximumintensity is located at the bonds linking the nearest andnext to the nearest neighbors of the impurity.Fig. 6 show the results for the second hybridizationgiven in Eq. (22). For both these two kinds of hybridiza-tions, one can find that the results for p -wave pairing and FIG. 6. (color online) LDOS for the case with the second hybridization ¯ V (2) k , the impurity is located at the center of theemergent honeycomb lattice, R = (0 , s -wave, extended s -wave, p -wave and d + id pairing, respectively. The gray lines labeled with ‘+’ in (a)-(d) show the clean LDOS without impurity. The red lines labeledwith ‘ × ’ in (a)-(d) show the LDOS near the impurity R = (0 , R = ( − / √ , p -wave pairing).(f) and (h) present the real space distributions of the LDOS of the two in-gap resonance peaks for d + id pairing shown in (d).The silvery lines in (e)-(j) show the emergent honeycomb lattice. The parameters are chosen as follows: ∆ s = ∆ s (cid:48) = ∆ p = 0 . d + id = 0 . (cid:15) d = − ∆ s / V = ∆ s for the s -wave (a) and extended s -wave (b) pairing, V = ∆ s / p -wave pairing,and V = 2∆ s for the d + id pairing (d). d + id pairing exhibit only quantitative differences. Herewe present these similarities. At the position close tothe impurity, R = (0 , × ’ inthese figures), at the nearest site close to the impurity, R = ( − / √ , p -wave and d + id ), the two kinds of hybridizations shown in Figs. 5and 6 may be distinguishable, i.e., the maximum inten-sity of the negative resonance energy is located at theimpurity site for the first hybridization, it is located onthe bonds link the nearest neighbors of the impurity forthe second hybridization.Now we analysis the results for the third hybridizationgiven in Eq. (23). The first column in Fig. 7 showsthe LDOSs versus energy for different pairing symme-tries and at different positions. One can find that, inthis case, all of the four kinds of pairing symmetries ex-hibit two in-gap resonance peaks at the original point, R = (0 ,
0) (see the red lines labeled with ‘ × ’ in Figs.7(a)-7(d) and inserts therein). Another important dif-ference is that, as shown in Figs. 7(e)-7(k), the spacialdistributions of the LDOSs break both the threefold ro-tation symmetry along the z -axis and the twofold rota-tion symmetry along the y -axis. Like the first two cases,the spacial distributions for the p -wave pairing and the d + id pairing are very similar to each other, which demon-strates that the p -wave pairing and the d + id pairing aredifficult to be distinguished via the LDOS measurementof magnetic impurity resonance states. However, the spa-cial distributions of p -wave pairing and the d + id pairingat the negative resonance energy is significant differentfrom the s -wave pairing and the extended s -wave pair-ing, i.e., the maximum intensity points for the latter arelocated at R = (0 , R ≈ (0 . , − . µ = − . − / µ = 1 . / µ = 1 . µ = − . µ = − .
7, andthe results for µ = 1 . µ = 1 .
5. Which demonstrate that the Van Hovesingularities do not have significant influence on LDOSof resonance states induced by the magnetic impurity.Fig. 9 gives an integrated presentation of LDOSs for µ = 1 . µ = − . FIG. 7. (color online) LDOS for the case with the third hy-bridization ¯ V (3) k , the impurity is located at the center of oneof the AB spot, R = (0 , s -wave, extended s -wave, p -wave and d + id pair-ing, respectively. The gray lines labeled with ‘+’ in (a)-(d)show the clean LDOS without impurity. The red lines la-beled with ‘ × ’ in (a)-(d) show the LDOS near the impurity R = (0 , R = ( − / √ , s = ∆ s (cid:48) = ∆ p = 0 .
08, ∆ d + id = 0 . (cid:15) d = − ∆ s / V = ∆ s for the s -wave (a) and extended s -wave(b) pairing, V = ∆ s / p -wave pairing, and V = 2∆ s for the d + id pairing (d). the two orbitals at one site of the emergent honeycomblattice, the spacial distributions shown in the last tworows in Fig. 9 break both the threefold rotation symme-try along the z -axis and the twofold rotation symmetryalong the y -axis, when the impurity is located at the cen-ter of the emergent honeycomb lattice and hybridized tothe nearest neighbors symmetrically, the spacial distri-butions of the resonance states exhibit sixfold rotationsymmetry, 3) for the impurity location given in 2), theLDOS versus energy for s -wave pairing and extended s -wave pairing at R = (0 ,
0) do not have in-gap resonancepeak, and 4) for the impurity location given in 2), theLDOS versus energy for p -wave pairing and d + id pairing FIG. 8. (color online) LDOS versus energy E for the casewith the third hybridization ¯ V (3) k . (a) s -wave pairing, (b)extended s -wave pairing, (c) p -wave pairing, (d) d + id pair-ing. The green lines, red lines labeled with ‘+’, and theblue dashed lines present the LDOS for R = ( − / √ , R = (1 / √ , / R = (1 / √ , − / at R = (0 ,
0) have only one in-gap resonance peak be-low the Fermi energy. The maximum point in the spacialdistributions of the resonance states need to be empha-sised here. For the third hybridization, ¯ V (3) k given in Eq.(23), the spacial distributions of the LDOSs show simi-lar characteristics for different chemical potentials, see,i.e., Fig. 7 and the last two rows of Fig. 9 for details.For the other two hybridizations, however, the resultsare different. Comparing all the contour plots in Figs.5 and 6, one can find that there are only two subfigureswhere the maximum points in the spacial distributionsof the resonance states are located at R = (0 , V (1) k hybridization for p -wave and d + id pairing symmetries. When the chem-ical potential is doped to close to the other half fillingratio 1 /
2, i.e., µ = 1 .
308 shown in Fig. 9, we find thatin which the maximum points located at R = (0 ,
0) areFigs. 9(42) and 9(45). Both of them correspond to the¯ V (2) k hybridization. Experimentally, when the impurityis located at the center of the emergent honeycomb lat-tice, the maximum point located at R = (0 ,
0) is a strongevidence for unconventional pairing symmetry, however,which one of the two hybridizations is preferred for aspecified magnetic impurity is difficult to discriminate.Our calculations suggest that, if the pairing symmetry is p -wave or d + id , by tuning the chemical potential close tothe two half filling ratios ± /
2, there must be a situationwhere the maximum point in the spacial distributions
FIG. 9. (color online) Integrated presentation of the LDOSs for µ = 1 . / V (1) k . The two rows in the middle show the results forthe second hybridization ¯ V (2) k . The last two rows show the results for the third hybridization ¯ V (3) k . (11), (31) and (51) showLDOSs versus energy for the s -wave pairing with different hybridizations. (14), (34) and (54) show LDOSs versus energy forthe extended s -wave pairing with different hybridizations. (21), (41) and (61) show LDOSs for the p -wave pairing. (24), (44)and (64) show the results for d + id pairing. In each line plot, three situations are presented. The gray lines labeled with ‘+’show the clean LDOSs without impurity. The red lines labeled with ‘ × ’ show LDOSs at R = (0 , R = ( − / √ , E = − . E = 0 . of the resonance states is located at R = (0 ,
0) for thenegative resonance energy.Now we consider the influence of the term proportionalto t (cid:48)(cid:48) given in Eq. (12). We calculate the LDOS versusenergy at R = (0 ,
0) for (1) three different values of t (cid:48)(cid:48) , t (cid:48)(cid:48) = 0, t , and 2 t , (2) two different topological nontriv-ial pairing, the p -wave pairing and d + id pairing, (3) the three kind of hybridizations, and (4) the four differentchemical potentials: filling ratio = ± / t (cid:48)(cid:48) , so we suspect that theeffect of the t (cid:48)(cid:48) term is limited.0 IV. CONCLUSION
Based on the two orbital model proposed in Ref. [42],we give a systematic study of the resonance states near amagnetic impurity in superconducting TwBG for typicalpairing symmetries. We need to emphasize that thoughthe parameters chosen in this work is different from thosegiven in Ref. [42] (especially for the case t (cid:48)(cid:48) (cid:54) = 0 and t (cid:48) = t (cid:48) = 0 used to fit the band structure given in Ref. [55]),our results are general and available for the parametersgiven in Ref. [42]. Here are our main results from bothsymmetry analysis and numerical calculations:1). When the impurity is hybridized to the two orbitalsat one site of the emergent honeycomb lattice, i.e., theimpurity is located at the center of the AB or BA spot,for any pairing symmetry belongs to the irreducible rep-resentations of D point group, the spacial distributionof the resonance states will break both the threefold rota-tion symmetry around the z -axis and the twofold rotationsymmetry around the y -axis.2). When the impurity is located at the center of theemergent honeycomb lattice, i.e., at the center of theAA spot, and hybridized to the six nearest neighborssymmetrically, there are only two kinds of hybridizationswhich belong to the two-dimensional irreducible repre-sentations of the D point group. For each hybridizationand the four pairing symmetries studied, the spacial dis-tributions of the resonance states reveal sixfold rotationsymmetry.3). For the hybridizations given in point 2), the in-gap resonance peak at the position close to the impuritymust be vanishing for the s -wave pairing, extended s -wave pairing, and any other pairing symmetry belongsto the A representation of D group.4). For the hybridizations given in point 2), the uniqueresonance peak with negative resonance resonance energyat the position close to the impurity indicate that thepairing potential is topological nontrivial p -wave pairingor d + id pairing. Here we give a discussion about these results, focus-ing on further investigations. Firstly, these conclusions,especially point 1) and point 2), are essentially depen-dent on the assumption given in Ref. [42] that the twoorbitals belong to the two-dimensional irreducible rep-resentation of D group. Further investigations basedon other models [3, 4, 56–59] will be important refer-ence to find the pairing symmetry of superconductingmagic angle TwBG. Secondly, the pairing symmetrieswe considered here are not complete. When the inter-orbital pairing is considered, there will be lots of possi-ble pairing potentials. For the on-site pairing, it is easyto check that there are 2 kinds of pairing potentials foreach irreducible representation. For the nearest neigh-bor pairing and next to nearest neighbor pairing, thereare numerous pairing potentials for each irreducible rep-resentation, a complete analysis for each pairing poten-tial will be straightforward. Thirdly, as shown in theinserts of Figs. 6(c), 6(d), 9(21) and 9(24), even if thepairing potential is p -wave or d + id symmetric, the in-gap resonance peak at R = (0 ,
0) may be too weak tobe detected. Fortunately, the in-gap resonance peak issignificant for the charge density close to the other half-filling, see the inserts of Figs. 9(41), 9(44), 5(c) and5(d), correspondingly. This is experimental accessible bytuning the gate voltage. Fourthly, when the impurity islocated at the center of the emergent honeycomb lattice,the two hybridizations are obtained from symmetry anal-ysis, for a practical adatom like manganese or chromium,the hybridization may break the rotation symmetries be-forehand. Systematic theoretical investigations based onfirst-principle calculations may help to find the properadatoms with desired hybridizations.
ACKNOWLEDGMENT
We appreciate the support from the NSFC underGrants No. 11504106 and No. 11447167 and the Funda-mental Research Funds for the Central Universities underGrant No. 2018MS049. [1] Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi,E. Kaxiras, and P. Jarillo-Herrero, Nature , 43(2018).[2] Y. Cao, V. Fatemi, A. Demir, S. Fang, S. L. Tomarken,J. Y. Luo, J. D. Sanchez-Yamagishi, K. Watanabe,T. Taniguchi, E. Kaxiras, R. C. Ashoori, and P. Jarillo-Herrero, Nature , 80 (2018).[3] J. M. B. Lopes dos Santos, N. M. R. Peres, and A. H.Castro Neto, Phys. Rev. Lett. , 256802 (2007).[4] R. Bistritzer and A. H. MacDonald, PNAS , 12233(2011).[5] E. Su´arez Morell, J. D. Correa, P. Vargas, M. Pacheco,and Z. Barticevic, Phys. Rev. B , 121407 (2010). [6] G. Trambly de Laissardi`ere, D. Mayou, and L. Magaud,Phys. Rev. B , 125413 (2012).[7] S. Fang and E. Kaxiras, Phys. Rev. B , 235153 (2016).[8] C. Xu and L. Balents, arXiv:1803.08057 (2018).[9] H. Guo, X. Zhu, S. Feng, and R. T. Scslettar,arXiv:1804.00159 (2018).[10] T. Huang, L. Zhang, and T. Ma, arXiv:1804.06096(2018).[11] C.-C. Liu, L.-D. Zhang, W.-Q. Chen, and F. Yang,arXiv:1804.10009 (2018).[12] D. Kennes, J. Lischner, and C. Karrasch,arXiv:1805.06310 (2018).[13] M. Fidrysiak, M. Zegrodnik, and J. Spa(cid:32)lek,arXiv:1805.01179 (2018). [14] B. Roy and V. Juriˇci´c, arXiv:1807.05524 (2018).[15] Y. Sherkunov and J. J. Betouras, arXiv:1807.05524(2018).[16] F. Wu, A. H. MacDonald, and I. Martin,arXiv:1805.08735 (2018).[17] B. Lian, Z. Wang, and B. A. Bernevig, arXiv:1807.04382(2018).[18] T. J. Peltonen, R. Ojaj¨arvi, and T. T. Heikkil¨a,arXiv:1805.01039 (2018).[19] C. C. Tsuei, J. R. Kirtley, C. C. Chi, L. S. Yu-Jahnes,A. Gupta, T. Shaw, J. Z. Sun, and M. B. Ketchen, Phys.Rev. Lett. , 593 (1994).[20] C. C. Tsuei, J. R. Kirtley, Z. F. Ren, J. H. Wang,H. Raffy, and Z. Z. Li, Nature , 481 (1997).[21] J. E. Hoffman, K. McElroy, D.-H. Lee, K. M. Lang,H. Eisaki, S. Uchida, and J. C. Davis, Science , 1148(2002).[22] Q.-H. Wang and D.-H. Lee, Phys. Rev. B , 020511(2003).[23] T. Hanaguri, Y. Kohsaka, J. C. Davis, C. Lupien, I. Ya-mada, M. Azuma, M. Takano, K. Ohishi, M. Ono, andH. Takagi, Nat. Phys. , 865 (2007).[24] T. Hanaguri, Y. Kohsaka, M. Ono, M. Maltseva, P. Cole-man, I. Yamada, M. Azuma, M. Takano, K. Ohishi, andH. Takagi, Science , 923 (2009).[25] S. H. Pan, E. W. Hudson, K. M. Lang, H. Eisaki,S. Uchida, and J. C. Davis, Nature , 746 (2000).[26] E. W. Hudson, K. M. Lang, V. Madhavan, S. H. Pan,H. Eisaki, S. Uchida, and J. C. Davis, Nature , 920(2001).[27] J. Bobroff, H. Alloul, W. A. MacFarlane, P. Mendels,N. Blanchard, G. Collin, and J.-F. Marucco, Phys. Rev.Lett. , 4116 (2001).[28] G.-M. Zhang, H. Hu, and L. Yu, Phys. Rev. Lett. ,704 (2001).[29] A. Polkovnikov, S. Sachdev, and M. Vojta, Phys. Rev.Lett. , 296 (2001).[30] J.-X. Zhu and C. S. Ting, Phys. Rev. B , 020506(2000).[31] M. Vojta, R. Zitzler, R. Bulla, and T. Pruschke, Phys.Rev. B , 134527 (2002).[32] A. Polkovnikov, Phys. Rev. B , 064503 (2002).[33] X. Dai and Z. Wang, Phys. Rev. B , 180507 (2003).[34] S. Baar, N. Momono, K. Kawamura, Y. Kobayashi,S. Iwasaki, T. Sakawaki, Y. Amakai, H. Takano, T. Kuro-sawa, M. Oda, and M. Ido, J. Supercond. Nov. Magn. , 659 (2016).[35] W.-F. Tsai, Y.-Y. Zhang, C. Fang, and J. Hu, Phys.Rev. B , 064513 (2009).[36] Y. Bang, H.-Y. Choi, and H. Won, Phys. Rev. B ,054529 (2009).[37] A. Akbari, I. Eremin, and P. Thalmeier, Phys. Rev. B , 014524 (2010).[38] J. D. Sau and E. Demler, Phys. Rev. B , 205402 (2013).[39] Z.-G. Fu, P. Zhang, Z. Wang, and S.-S. Li, J. Phys.:Condens. Matter , 145502 (2012).[40] G.-Q. Zha and Y.-Y. Jin, Europhys. Lett. , 27002(2017).[41] Y.-W. Guo, W. Li, and Y. Chen, Front. Phys. , 127403(2017).[42] N. F. Q. Yuan and L. Fu, Phys. Rev. B , 045103 (2018).[43] H. Isobe, N. F. Q. Yuan, and L. Fu, arXiv:1805.06449(2018).[44] D. Withoff and E. Fradkin, Phys. Rev. Lett. , 1835(1990).[45] C. R. Cassanello and E. Fradkin, Phys. Rev. B , 15079(1996).[46] C. R. Cassanello and E. Fradkin, Phys. Rev. B , 11246(1997).[47] C. Gonzalez-Buxton and K. Ingersent, Phys. Rev. B ,14254 (1998).[48] H.-B. Zhuang, Q.-F. Sun, and X. C. Xie, Europhys. Lett. , 58004 (2009).[49] B. Uchoa, T. G. Rappoport, and A. H. Castro Neto,Phys. Rev. Lett. , 016801 (2011).[50] L. Yu, Acta Phys. Sin. , 75 (1965).[51] H. Shiba, Prog. Theor. Phys. , 435 (1968).[52] A. I. Rusinov, JETP Lett. , 85 (1969).[53] D. A. Ivanov, Phys. Rev. Lett. , 268 (2001).[54] H. Hu, L. Jiang, H. Pu, Y. Chen, and X.-J. Liu, Phys.Rev. Lett. , 020401 (2013).[55] N. N. T. Nam and M. Koshino, Phys. Rev. B , 075311(2017).[56] J. M. B. Lopes dos Santos, N. M. R. Peres, and A. H.Castro Neto, Phys. Rev. B , 155449 (2012).[57] P. Moon and M. Koshino, Phys. Rev. B87