Introduction. In the Bardeen-Cooper-Schrieffer (BCS) theory, when electrons experience attractive interactions mediated by phonons, a superconductor (SC) pairing instability can arise, leading to a gap opening at the electron Fermi surface [1]. Typically, the Fermi surface becomes fully gapped to lower the SC ground state energy, resulting in nodeless swave superconductivity. In contrast, nodal superconductivity can emerge when the SC pairing is of an "unconventional" nature [2]. The most notable example is the high T c superconductivity in cuprates [3], where the d-wave nature of SC pairing leads to nodes at the Fermi surface. The occurrence of nodal superconductivity usually requires unconventional SC pairing mechanism or wavefunction "topology," beyond the phonon mediated superconductivity in normal metals in the BCS theory, as the nodes at the Fermi surface are not energetically favored. However, this scenario might be changed for superconductivity in multiple-flat-band systems when the band width and energy separation between different flat bands are smaller than the attractive interactions between electrons (strong coupling limit), even for the phonon-mediated SC systems, which is the main topic of this work. This topic could also be relevant for the Fermi surfaces with topological obstructions in the presence of electron-phonon interaction.

Flat-band superconductivity has attracted great research interests due to the discovery of superconductivity in Moiré SC materials, e.g., twisted bilayer/trilayer graphene [4,5]. While the microscopic origin of superconductivity in these Moiré SC materials still remains elusive , the experimental studies of the SC gap in twisted bilayer graphene (TBG) suggest an unconventional nature of SC pairing, beyond the conventional s-wave spin-singlet pairing in the BCS theory. The scanning tunneling microscope (STM) [39] and electronic thermal conductance measurements [40] on TBG suggested a nodal SC pairing gap; evidence of nematicity has also been found in the SC state via magnetotransport measurements [41]. On the other hand, recent Angle-resolved photoemission spectroscopy (ARPES) experiments with high spatial resolutions found the replicas of flat bands that are expected to result from the strong electron-K-phonon interaction [42]. This raises the academic questions if phonon-mediated superconductivity can reveal the unconventional features such as nodes and nematicity in the SC state, which can be more stable than the fully gaped SC state, in flat-band land. Note that we are not claiming this as the TBG superconductivity mechanism due to the strong Coulomb interaction, which needs to be screened for SC to appear.

In Ref. [43], we have shown that the electron-phonon interaction between the bare bands of TBG and the K phonons [35,44] can induce two types of pairing channels, the intra-Chern-band and inter-Chern-band pairings, which share a very close T c , based on the self-consistent gap equation in the heavy fermion model without Coulomb interaction [45][46][47]. A recent study in Ref. [33] suggested that the Kondo screening due to the Coulomb interaction can help the intra-Chern-band pairing channel to overcome the inter-Chern-band s-wave pairing channel. The intra-Chern-band pairing is d wave and spin singlet, and belongs to the 2D E 2 irreducible representation (irrep) of the C 6v group. Strikingly, a Euler obstructed pairing that belongs to the intra-Chern-band pairing has recently been proposed to understand the nematicity and nodal superconductivity in TBG [18,33]. Even though the SC nature in TBG remains elusive, the academic question of the nematic and nodal SC is relevant as the evidence of these features have been found in TBG experiments [39][40][41]. Due to the 2D irrep nature, other intra-Chern-band pairing channels, e.g., chiral d-wave pairing, can also exist and share the same T c as the Euler obstructed pairing. Naively, one may anticipate chiral d-wave pairing is more stable at zero temperature as it is fully gapped while the nematic Euler obstructed pairing has nodes. Such situation occurs in the doped graphene model [48,49]. Thus, it is interesting to theoretically understand the stability of nematic Euler obstructed pairing with nodes.

In this work, we study the low temperature behavior of the 2D E 2 -irrep intra-Chern-band pairing induced by the electron-K-phonon interaction in the TBG heavy fermion model [45,46]. Our numerical results of the condensation energy suggest the nematic Euler obstructed pairing has a lower energy than the chiral d-wave pairing and has nodes for a wide range of parameters in the chiral limit, although this limit is away from experimental situations. A detailed theoretical analysis shows the nematic Euler obstructed pairing is at least locally stable. The underlying mechanism is that although gapping out the nodes in Euler obstructed pairing can lower energy around the nodal regions in the momentum space, strong inter-eigenband pairing potential can lead to an energy saving for other momenta away from the nodal regions, which is the dominant contribution. This mechanism is proposed to represent a counter example to the known fact that gapping gapless points leads to lower energy.

Bogoliubov-de Gennes Hamiltonian for the flat bands of TBG. Our previous work [43] has shown that in absence of Coulomb interaction, the coupling between the flat bands and the K phonons in the TBG model can potentially favor two pairing channels, the s-wave spin-singlet inter-Chern-band pairings with 1D A 1 irreducible representation (irrep) of C 6v group and the d-wave spin-singlet intra-Chern-band pairings with the 2D E 2 irrep. These two channels are decoupled in the chiral limit and their T c 's are close (for the bare flat bands). We focus on the d-wave intra-Chern-band channel below. Due to the 2D E 2 irrep nature, all different forms of the d-wave intra-Chern-band pairings share the same T c . The degeneracy of different 2D E 2 irrep pairing channels will be broken down below T c , which can be extracted from the full gap equation.

To study that, we consider the Bogoliubov-de Gennes (BdG) Hamiltonian of TBG with spin singlet pairing,

where η = ± is the valley index and s 0 and s label the identity and Pauli matrices for spin. The basis functions are

, where γ † k,e Y ,η,s is the electron creation operator for the flat bands in TBG with the Chern-band index e Y = ± and spin s. Here the Chern-band basis is related to the eigenstate basis by

for the projected single-particle Hamiltonian, where γ † k,η,n,s is the electron creation operator for the eigen-state n = ± of two flat bands per valley per spin. We use the Chern-band basis to write down the BdG Hamiltonian because the inter-Chern-band and intra-Chern-band pairings are decoupled in the chiral limit [43,50,51]. We note that the realistic TBG is away from the chiral limit, and thus our discussion is of academic interest. On the Chern-band basis, the projected single-particle Hamiltonian on the flat bands at the valley η reads

where μ is the chemical potential, ζ 0 and ζ x,y,z represent the identity matrix and Pauli matrices for the Chern-band basis,

, and n,η (k) are the single-particle eigenenergies that can be extracted from the Bistritzer-MacDonald (BM) model [43,52]. Time reversal (TR) symmetry requires h -η (-k) = h η (k). We always consider the chiral limit with d 0,η (k) = 0 below. The intra-Chern-band pairing has the form

where k,e Y labels the gap function on the Chern-band channel e Y = ±. The intra-Chern-band pairing gap function ( k,+ , k,-) belongs to the 2D E 2 irrep and has a d-wave nature, which can be represented by [43] k

with k -e Y = k xie Y k y . In particular, ( + , -) = ( 0 , 0) or (0, 0 ) gives the TR-breaking chiral d wave while ( + , -) = ( 0 , * 0 ) represents the TR-preserving nematic Euler obstructed pairing. H - BdG (k) can be related to H + BdG (k) by SC particle-hole symmetry (See Sec. A in the Supplemental Material [53]), so we focus on H + BdG (k), which has a block diagonal form in the spin space [see Eq. ( 1)] and we can further decompose it into

BdG a 4×4 matrix Hamiltonian (See Eq. (A10) in Sec. A in the Supplemental Material for H +,λ BdG [53]). The eigenenergy of H +,λ=+ BdG can be solved analytically as

where 2 ) = 0. The general conditions for the nodes are discussed in Ref. [43] and Sec. A in the Supplemental Material [53]. The eigenstates of H +,+ BdG can be encoded in the unitary matrix U k defined by where

Gap equation and superconductor phase diagram. To derive the self-consistent gap equation below T c , we consider the attractive electron-electron interaction mediated by K phonons [43],

where

, and V 0 = 0.4 meV and b = 0.185k θ are material dependent parameters given in Ref. [43] (also see Sec. A in the Supplemental Material [53]).

With the above form of V interaction and gap function ansatz in Eq. ( 5), the self-consistent gap equation is 9) with N M the number of Moiré unit cells. Using the gap equation ( 9) at zero temperature, the condensation energy is given by

where ξ - k,n = -1 2 |d 0,k -μ + nd x,k | is the eigenenergy of the single-particle Hamiltonian.

The SC phase diagram can be extracted by solving the self-consistent gap equation (9). Figures 1(a To investigate the stability of the nematic Euler obstructed pairing, we consider a general form of the gap function,

where 0 is real, the parameter α tuning the gap function form in 2D E 2 irrep, and φ is the relative phase. The nematic Euler obstructed pairing corresponds to α = π 4 while the chiral d-wave pairing corresponds to α = 0 or π 2 . The BdG spectrum is generally gapped unless α = π/4. E c has a very weak dependence on φ (See Fig. S3(B) and the corresponding discussion in Sec. B in the Supplemental Material [53]), and thus we choose φ = 0 as an example to show the 0 and α dependence of E c in Fig. 2(a), in which the red dashed line shows the minimal value of E c located at 0 ≈ 0.22 meV. This value of 0 is consistent with that obtained from the self-consistent gap equation ( 9)

shows E c as a function of α at 0 ≈ 0.22 meV, in which the minima of E c appears at α = π/4. These numerical results suggest the Euler obstructed pairing has a lower energy than the chiral d-wave pairing with α = 0 and is at least locally stable. An argument about the stability of Euler obstructed pairing in the flat band limit is provided in Sec. C in the Supplemental Material [53] (also see Ref. [43]). We will next provide a 12) and ( 13). (d) The energy gap of the BdG spectrum is plotted as a function of k for α = π/4. The energy unit here is meV. theoretical understanding of this local stability when miniband dispersion is considered.

Locally stable nematic Euler obstructed pairing. To understand this local stability, we consider the gap function slightly away from the Euler obstructed pairing and study the variation of the condensation energy at zero temperature. Particularly, we choose α = π 4 + δ 0 with δ 0 1 in the gap function ansatz (11) and expand the condensation energy E c in Eq. ( 10) as

up to the δ 2 0 order. Our main result below is to demonstrate γ = 0 and β > 0 so that Euler obstructed pairing becomes locally stable. Furthermore, we define

and the momentum distribution of β k is shown in Fig. 2(c), from which one can see that the negative value of β k around the nodal regions is compensated by the positive β k away from the nodal regions. Now let us discuss different terms of E c in Eq. (10). The term

is subtle due to the existence of the nodes in the BdG spectrum. The location of the nodes in the momentum space is denoted as k g , whose values are determined by the conditions cos(

We divide the Moiré Brillouin zone (MBZ) into two regions: one region is around the point node k g , denoted as k g , namely k = k g + q with |q| < ε, where ε serves as a momentum cutoff, and the other is the momentum region away from the point nodes, namely k ∈ k g . In the region k ∈ k g , we can expand the BdG eigenenergy for the n =branch around the nodal points and find

for small q and δ 0 , where q 1 = 1 √ 2 (q x + q y ), q 2 = 1 √ 2 (q xq y ), and v 1,2 and B 3 are the parameters extracted from the perturbation expansion (see Sec. C in the Supplemental Material [53]). In Eq. ( 15), E - k,-is gapless at q = 0, δ 0 = 0. Although a small but finite δ 0 will give a gap that is linearly proportional to |δ 0 |, we find its contribution to the condensation energy E c requires the momentum q integral of Eq. ( 15) in the region k g , and this integral gives the term of order δ 2 0 (see Sec. C in the Supplemental Material [53]) under the condition

where

is the average velocity around the point nodes. In our numerical calculations, we first choose the cutoff ε so that the energy dispersion (15) is valid in k g and then choose δ 0 to satisfy Eq. ( 16). The correction to β in the region k g is approximately given by (see Sec. C in the Supplemental Material [53])

where S M is the area of the Moiré unit cell. The negative sign here suggests that this contribution will lower E c for a finite δ 0 .

In numerical calculations, we choose ε ≈ 0.01k θ , so that the area ratio between k g and MBZ, represented by the factor πε 2 S M (2π ) 2 in Eq. ( 17), is only around 10 -4 , and the energy dispersion within k g is well described by Eq. (15). With this value of ε, we then set δ 0 = 0.01 × π 2 ≈ 0.0157 to satisfy Eq. ( 16) as ε ≈ 0.01k θ |B 3 δ 0 | |v| ≈ 0.0014k θ with v ≈ 0.4 meV/k θ and B 3 ≈ 0.036 meV. The contribution to β from the k g region is estimated as k∈ kg β k ≈ -1.6 × 10 -4 meV.

For the n = + branch in the whole MBZ or the n =branch in the momentum region k ∈ k g , the BdG spectrum has a gap, so we can perform the perturbation expansion for E - k,n for a small δ 0 ,

up to the order of δ 2 0 , where

and

for n = ±. From the above derivation, one can see that the lowest order contribution to E c is of order δ 2 0 , no matter inside or outside k g , so γ = 0 in Eq. (12). The full form of β in Eq. ( 12) is given by

in which the first and second terms are the contributions inside and outside k g , respectively. We numerically find β k ∈ kg of order 0.1 meV, several orders larger than the contribution inside k g (∼10 -4 meV). Thus, it is safe to neglect the first term in β. We numerically evaluate β k from the second term in Eq. ( 21), which is shown in Fig. 2(c). By comparing to the BdG energy gap shown in 2(d), we find a negative β k around the point nodes, while β k for a large momentum is always positive. All these features of β k can be understood from Eq. ( 21) (see Sec. C in the Supplemental Material [53]). After summing over the momenta for β k , our numerical results show that the positive β k for large momenta is dominant and β ≈ 0.44 meV > 0. This concludes the Euler obstructed pairing for α = π/4 is locally stable. The locally stable nematic Euler obstructed pairing in our model is in sharp contrast to the case of the doped graphene model [48,49], in which the chiral d-wave pairing with a full gap is found to be more stable. Unlike the dispersive bands in doped graphene, both the bandwidth and the energy splitting between different bands are small for two flat bands (per valley per spin) in the TBG model. Thus, the pairing potential not only exists within one eigenband (intra-eigenband pairing), but also can strongly mix different eigenbands (intereigenband pairing). This inter-eigenband mixing due to the pairing term is much smaller than the energy splitting between two bands in doped graphene. To see the influence of the intereigenband pairing, we manually turn off the inter-eigenband pairing term and perform a similar analysis of E c expanded around α = π/4 + δ 0 (δ 0 1) in Sec. D in the Supplemental Material [53]. Our numerical calculation gives a negative β in the δ 2 0 term in Eq. ( 12), which implies that the Euler obstructed pairing becomes unstable. Figure S6 in Sec. D in the Supplemental Material [53] shows the minima of E c in this case is deviated from the Euler obstructed pairing with α = π/4, which reveals the possibility of a phase transition from time-reversal preserving Euler SC phase to timereversal breaking SC phase via reducing the inter-eigenband pairing.

Conclusion. This work suggests that the gapless nematic Euler obstructed pairing can be locally stable in the nearly degenerate flat band systems as the inter-eigenband pairing potential, which is comparable to the bandwidth of each band and band gap between two flat bands, can induce a strong interband mixing. We emphasize that the pairing mechanism for the superconductivity in TBG remains elusive. Our work begins with an attractive interaction that supports the intra-Chern-band pairing, and once this pairing channel is energetically favored, a similar theoretical framework can be applied, making a nematic superconducting phase possible regardless of whether the exact pairing mechanism arises from electron-phonon interactions or other types of interactions [37]. The mechanism discussed in this work is different from other mechanisms for nematic SC phase in TBG, e.g., strong density wave fluctuation mechanism [30] and the pair mixing mechanism [15], and it may also be applicable to other multiple-flat-band SC systems. Due to the existence of bulk nodes of the Euler obstructed pairing, we expect the emergence of Majorana flat bands located at the 1D boundary [55,56], which can generally be revealed as a zero-bias peak at the sample boundary in the STM measurements.