Valley Jahn-Teller effect in twisted bilayer graphene

TL;DR

This paper explores how specific phonon modes, especially a Jahn-Teller coupled mode, influence the insulating and superconducting states in twisted bilayer graphene, highlighting a phonon-driven mechanism for these phenomena.

Contribution

It identifies a particular phonon mode with strong Jahn-Teller coupling that can induce insulating states and potentially facilitate superconductivity in twisted bilayer graphene.

Findings

Discovery of a special phonon mode with strong Jahn-Teller coupling.

Identification of phonon-induced insulating states at various electronic fillings.

Proposal of phonon-mediated pairing leading to extended s-wave or d+id superconductivity.

Abstract

The surprising insulating and superconducting states of narrow-band graphene twisted bilayers have been mostly discussed so far in terms of strong electron correlation, with little or no attention to phonons and electron-phonon effects. We found that, among the 33492 phonons of a fully relaxed $\theta=1.08^\circ$ twisted bilayer, there are few special, hard and nearly dispersionless modes that resemble global vibrations of the moir\'e supercell, as if it were a single, ultralarge molecule. One of them, doubly degenerate at $\Gamma$ with symmetry $A_1+B_1$, couples very strongly with the valley degrees of freedom, also doubly degenerate, realizing a so-called $\text{E}\otimes\text{e}$ Jahn-Teller (JT) coupling. The JT coupling lifts very efficiently all degeneracies which arise from the valley symmetry, and may lead, for an average atomic displacement as small as $0.5~$mA, to an insulating state at charge neutrality. This insulator possesses a non-trivial topology testified by the odd winding of the Wilson loop. In addition, freezing the same phonon at a zone boundary point brings about insulating states at most integer occupancies of the four ultra-flat electronic bands. Following that line, we further study the properties of the superconducting state that might be stabilized by these modes. Since the JT coupling modulates the hopping between AB and BA stacked regions, pairing occurs in the spin-singlet Cooper channel at the inter-(AB-BA) scale, which may condense a superconducting order parameter in the extended $s$-wave and/or $d\pm id$-wave symmetry.

Figures

Tables

Table 1. Table 1: Nontrivial irreducible representations of the space group D 6 subscript 𝐷 6 D_{6} . Each representation has the degeneracy shown in parenthesis. We also list the action of the symmetry operations for each representation, where ϕ = 2 ​ π / 3 italic-ϕ 2 𝜋 3 \phi=2\pi/3 .

	${\text{C}}_{3z}$	${\text{C}}_{2x}$	${\text{C}}_{2y}$
A${}_1{}^{}\left(1\right)$	+1	+1	+1
A${}_2{}^{}\left(1\right)$	+1	-1	-1
B${}_1{}^{}\left(1\right)$	+1	+1	-1
B${}_2{}^{}\left(1\right)$	+1	-1	+1
E${}_1{}^{}\left(2\right)$
E${}_2{}^{}\left(2\right)$

Equations

$$\Delta N_v$$

$${\Psi }_{\alpha ,\mathbf{R}\sigma }^{}$$

$${\Psi }_{\mathbf{R}\sigma }^{}$$

$$\Delta N_v$$

$$\displaystyle\mathcal{C}_{2x}\Big{(}\Psi^{\phantom{\dagger}}_{\mathbf{R}\sigma}\Big{)}$$

$$\displaystyle\mathcal{C}_{2x}\Big{(}\Phi^{\phantom{\dagger}}_{\mathbf{R}\sigma}\Big{)}$$

$$\displaystyle\mathcal{C}_{2z}\Big{(}\Psi^{\phantom{\dagger}}_{\mathbf{R}\sigma}\Big{)}$$

$$\displaystyle\mathcal{C}_{2z}\Big{(}\Phi^{\phantom{\dagger}}_{\mathbf{R}\sigma}\Big{)}$$

$$\displaystyle\mathcal{C}_{2y}\Big{(}\Psi^{\phantom{\dagger}}_{\mathbf{R}\sigma}\Big{)}$$

$$\displaystyle\mathcal{C}_{2y}\Big{(}\Phi^{\phantom{\dagger}}_{\mathbf{R}\sigma}\Big{)}$$

$$\displaystyle\mathcal{C}_{2y}\Big{(}\Psi^{\phantom{\dagger}}_{\mathbf{k}\sigma}\Big{)}=\sigma_{1}\,\tau_{1}\,\mu_{3}\,\Psi^{\phantom{\dagger}}_{\mathbf{k}\sigma}\,,$$

$$2T_3$$

$$\boldsymbol{\eta}(\mathbf{r})=\sum_{a}\,\sum_{i,j=1}^{3}\sin\big{(}\mathbf{q}_{ij}\cdot\mathbf{r}_{a}\big{)}\,\mathbf{u}_{a,ij}\,\delta\big{(}\mathbf{r}-\mathbf{r}_{a}\big{)}\,,$$

$$H_{\text{el-ph}}$$

$$J_3$$

$$H_{\text{el}}=-\Delta \sum _{\mathbf{R}\sigma }{\Psi }_{\mathbf{R}\sigma }^{†}{\sigma }_0{\tau }_0{\mu }_3{\Psi }_{\mathbf{R}\sigma }^{}+T_{\text{el}},$$

$$H_{\text{ph}}=\frac{\displaystyle\omega}{\displaystyle 2}\,\sum_{\mathbf{R}}\,\Big{(}\,\mathbf{p}_{\mathbf{R}}\cdot\mathbf{p}_{\mathbf{R}}+\mathbf{q}_{\mathbf{R}}\cdot\mathbf{q}_{\mathbf{R}}\,\Big{)}\,,$$

$$\gamma (\mathbf{k})$$

$${\gamma }_{+1}(\mathbf{k})$$

$${\gamma }_{-1}(\mathbf{k})$$

$$\gamma_{\pm 1}\Big{(}\text{C}_{3z}(\mathbf{k})\Big{)}=\omega^{\pm 1}\,\gamma_{\pm 1}(\mathbf{k})\,.$$

$$\left(\begin{array}{c}{\Phi }_{+1,\mathbf{k}\sigma }\\ {\Phi }_{-1,\mathbf{k}\sigma }\end{array}\right)=\frac{1}{\sqrt{2}}\left(\begin{array}{cc}1& -i\\ 1& +i\end{array}\right)\left(\begin{array}{c}{\Phi }_{s,\mathbf{k}\sigma }\\ {\Phi }_{p,\mathbf{k}\sigma }\end{array}\right),$$

$$\sum _{\sigma }\sigma {\Phi }_{+1,AB,\mathbf{k}\sigma }^{†}{\tau }_1{\Phi }_{+1,BA,-\mathbf{k}-\sigma }^{†}$$

$$\sum _{\sigma }\sigma {\Phi }_{-1,AB,\mathbf{k}\sigma }^{†}{\tau }_1{\Phi }_{-1,BA,-\mathbf{k}-\sigma }^{†}$$

$$\displaystyle\frac{\displaystyle 1}{\displaystyle\sqrt{2}}\sum_{\sigma}\sigma\Big{(}\Phi^{\dagger}_{+1,AB,\mathbf{k}\sigma}\,\tau_{1}\,\Phi^{\dagger}_{-1,BA,-\mathbf{k}-\sigma}$$

$$\displaystyle\qquad+\!(-)\Phi^{\dagger}_{-1,AB,\mathbf{k}\sigma}\,\tau_{1}\,\Phi^{\dagger}_{+1,BA,-\mathbf{k}-\sigma}\Big{)}$$

$$\frac{1}{\sqrt{2}}\sum _{\sigma }\sigma {\Psi }_{AB,\mathbf{k}\sigma }^{†}{\tau }_1{\mu }_{0(3)}{\Psi }_{BA,-\mathbf{k}-\sigma }^{†}$$

$$H_{\text{el-el}}$$

$$⟨{\text{A}}_{1\mathbf{k}}^{†}⟩$$

$$+\gamma (\mathbf{k})⟨A_{1,\mathbf{k}}⟩+\gamma (\mathbf{k})⟨A_{1,\mathbf{k}}^{\prime }⟩,$$

$$\displaystyle t(\mathbf{d})=V_{pp\sigma}(d)\bigg{[}\frac{\textbf{d}\cdot\textbf{e}_{z}}{d}\bigg{]}^{2}\!\!\!+V_{pp\pi}(d)\bigg{[}1-\Big{(}\frac{\textbf{d}\cdot\textbf{e}_{z}}{d}\Big{)}^{2}\bigg{]}\,,$$

$$V_{pp\sigma }(x)=V_{pp\sigma }^0{\text{e}}^{-\frac{x-d_0}{r_0}}{V}_{pp\pi }(x)=V_{pp\pi }^0{\text{e}}^{-\frac{x-a_0}{r_0}}$$

$$C_{\alpha \beta }(il,js)=\frac{{\partial }^{2}U}{\partial R_{\alpha il}\partial R_{\beta js}}$$

$$D_{\alpha \beta ij}(\mathbf{q})=\frac{1}{M_C}\sum _l{C}_{\alpha \beta }(il,j0)e^{-i\mathbf{q}{R}_l}$$

$$\sum _{j\beta }{D}_{\alpha \beta ij}(\mathbf{q}){ϵ}_{j\beta }(\mathbf{q})={\omega }_{\mathbf{q}}^2{ϵ}_{i\alpha }(\mathbf{q})$$

Full text

I Introduction

The recent discovery of superconductivity in magic angle ($\theta\approx 1.1^{\circ}$) twisted bilayer graphene (tBLG) Cao et al. (2018a, b); Yankowitz et al. (2018); Lu et al. (2019) has stimulated an intense theoretical and experimental research activity. This unexpected phenomenon occurs upon slightly doping insulating states found at fractional fillings, the latter contradicting the metallic behaviour predicted by band structure calculations. Despite the huge size of the unit cell, containing more than $\approx 11000$ atoms, the band structure of tBLG at the first magic angle has been computed with a variety of methods, including tight-binding Trambly de Laissardiére et al. (2010); Trambly de Laissardière et al. (2012); Shallcross et al. (2010); Sboychakov et al. (2015); Suárez Morell et al. (2010); Angeli et al. (2018), continuum models Bistritzer and MacDonald (2011); Tarnopolsky et al. (2018) and DFT Song et al. (2018); Lucignano et al. (2019). These approaches predict that all the exotic properties mentioned above arise from four extremely flat bands (FBs), located around the charge neutrality point, with a bandwidth of the order of $\approx 10-20$ meV. Owing to the flatness of these bands, the fractional filling insulators found in Cao et al. (2018a, b) are conjectured to be Mott Insulators, even though rather anomalous ones, since they turn frankly metallic above a critical temperature or above a threshold Zeeman splitting in a magnetic field, features not expected from a Mott insulator. Actually, the linear size of the unit cell at the magic angle is as large as $\approx 14$ nm, and the effective on-site Coulomb repulsion, the so-called Hubbard $U$, must be given by the charging energy in this large supercell projected onto the FBs, including screening effects due to the gates and to the other bands. Even neglecting the latter, the estimated $U\sim 9~{}\text{meV}$ is comparable to the bandwidth of the FBs Kang and Vafek (2018). Since the FBs are reproducibly found in experiments Lu et al. (2019); Cao et al. (2018a, b); Yankowitz et al. (2018); Kerelsky et al. (2018); Choi et al. (2019); Polshyn et al. (2019); Cao et al. (2019); Tomarken et al. (2019) to be separated from the other bands by a gap of around $\sim 30-50$ meV, the actual value of $U$ should be significantly smaller, implying that tBLG might not be more correlated than a single graphene sheet Wehling et al. (2011). In turn, this suggests that the insulating behaviour at $\nu=\pm 2$ occupancy might instead be the result of a weak-coupling Stoner or CDW band instability driven by electron-electron and/or electron-phonon interactions, rather than a Mott localization phenomenon.

As pointed out in Po et al. (2018), in order to open an insulating gap the band instability must break the twofold degeneracy at the $\mathbf{K}$-points imposed by the $D_{6}$ space group symmetry of the moiré superlattice, as well as the additional twofold degeneracy due to the so-called valley charge conservation. This conserved quantity is associated with an emergent dynamical $U(1)$ symmetry that appears at small twist angles, a symmetry which, unlike spatial symmetries, is rather subtle and elusive. It is therefore essential to identify a microscopic mechanism that could efficiently break this emergent symmetry, hereafter refereed to as $U_{v}(1)$. The most natural candidate is the Coulomb repulsion Po et al. (2018); Xie and MacDonald (2018); You and Vishwanath (2018), whose Fourier transform decays more slowly than that of electron hopping, possibly introducing a non negligible coupling among the two valleys even at small twist angles. Indeed DFT-based calculations show a tiny valley splitting Song et al. (2018); Lucignano et al. (2019), almost at the limits of accuracy of the method, which is nevertheless too small to explain the insulating states found in the FBs.

Here we uncover another $U_{v}(1)$-breaking mechanism involving instead the lattice degrees of freedom, mostly ignored so far. It must be recalled that ab-initio DFT-based calculations fail to predict well defined FBs separated from other bands unless atomic positions are allowed to relax Nam and Koshino (2017); Carr et al. (2018); Angeli et al. (2018); Lucignano et al. (2019); Carr et al. (2019); Choi and Choi (2018), especially out-of-plane. That alone already demonstrated that the effects of atomic motions in the lattice are not at all negligible in tBLG, further supported by the significant phonon contribution to transport Polshyn et al. (2019); Wu et al. (2018); Yudhistira et al. (2019). We calculate the phonon spectrum of the fully-relaxed bilayer at $1.08^{\circ}$ twist angle, which shows the presence, among the about thirty thousands phonons, of a small set of very special optical modes, with C-C stretching character, very narrow and uniquely coherent over the moiré supercell Brillouin Zone. Among them, we find a doubly degenerate optical mode that couples to the $U_{v}(1)$ symmetry much more efficiently than Coulomb repulsion seems to do in DFT calculations. A subsequent frozen-phonon tight-binding calculation shows that this mode is able to fully lift the valley degeneracy even when its lattice deformation amplitude is extremely small. Remarkably, both electrons and phonons are twofold $U_{v}(1)$ -degenerate, and the coupling of this mode with the electron bands actually realizes an $\text{E}\otimes\text{e}$ Jahn-Teller (JT) effect Englman (1972). This effect is able to stabilize insulating states at integer occupancies of the FBs, both even and odd. Moreover, a surprising and important additional result will be that the electron-phonon coupling magnitude controlling this process is extremely large, and not small as one could generally expect for a very narrow band. We conclude by studying the superconducting state that might be mediated by the Jahn-Teller coupling in a minimal tight-binding model of the FBs that reproduces symmetries and topological properties of the realistic band structure calculations.

The work is organized as follows. In section II we specify the geometry of the tBLG studied and define useful quantities that are used throughout the article. In section III we briefly discuss the band structure obtained by a realistic tight-binding calculation of the tBLG with fully relaxed atomic positions. The phonon spectrum and its properties, especially focusing on special optical modes strongly coupled with the valley $U_{v}(1)$ symmetry, are throughly discussed in section IV. Section VII addresses the properties of the superconducting state that might be stabilized by the particular phonon mode identified in the previous section, through a mean-field calculation using a model tight-binding Hamiltonian of the FBs. Finally, section VIII is devoted to concluding remarks.

II moiré superlattice and symmetries

In Fig. 1 we show the geometry of the tBLG that we shall use hereafter. We start from a AA stacked bilayer and rotate in opposite directions the two layers around the center of a hexagon by an angle $1.08^{\circ}/2$, leading to a unit cell with 11164 carbon atoms. The moiré superlattice, top left panel, with the reference frame defined in the top right panel, possesses a full $D_{6}$ spatial symmetry. Specifically, resolving the action of each symmetry operation in the indices that identify the two layers, 1 and 2, the two sublattices within each layer, A and B, and, finally, the two sublattices AB and BA of the moiré superlattice (see bottom panel of Fig. 1), we have that:

• •

the rotation $\text{C}_{3z}$ by $120^{\circ}$ degrees around the $z$-axis is diagonal in all indices, 1 and 2, A and B, AB and BA;

• •

the $\text{C}_{2x}$ rotation by $180^{\circ}$ degrees around the $x$-axis interchanges 1 with 2, A with B, but is diagonal in AB and BA;

• •

the $\text{C}_{2y}$ rotation by $180^{\circ}$ degrees around the $y$-axis interchanges 1 with 2, AB with BA, but is diagonal in A and B;

• •

finally, the action of a $\text{C}_{2z}$ rotation by $180^{\circ}$ degrees around the $z$-axis is a composite symmetry operation obtained by noting that $\text{C}_{2z}=\text{C}_{2x}\times\text{C}_{2y}$.

In Table 1 we list the irreducible representations (irreps) of the $D_{6}$ space group and the action on each of them of the symmetry transformations $\text{C}_{3z}$, $\text{C}_{2x}$ and $\text{C}_{2y}$.

II.1 $U_{v}(1)$ valley symmetry

The mini Brillouin zone (MBZ) that corresponds to the real space geometry of Fig. 1 and its relationship with the original graphene Brillouin zones are shown in the bottom panels of that figure (for better readability, at a larger angle than the actual $1.08^{\circ}$ which we use). Because of the chosen geometry, the Dirac point $K_{+}$($K^{\prime}_{+}$) of the top layer and $K^{\prime}_{-}$($K_{-}$) of the bottom one fold onto the same point, $K_{1}$ or $K_{2}$, of the MBZ, so that a finite matrix element of the Hamiltonian between them is allowed by symmetry. Nevertheless, as pointed out by Bistritzer and MacDonald (2011), the matrix element of the one-body component of the Hamiltonian is negligibly small at small twist angles, so that the two Dirac points, hereafter named valleys, remain effectively independent of each other. This implies that the operator

[TABLE]

where $N_{1}$ and $N_{2}$ are the occupation numbers of each valley, must commute with the non-interacting Hamiltonian of the tBLG at small angles. That operator is in fact the generator of the $U_{v}(1)$ symmetry. As we will see in the following (see III.1), the interplay between the valley symmetry and $C_{2y}$ is responsible for the additional two-fold degeneracies beyond those of Table 1 in both electron band structure and phonon spectra along all the points in the MBZ invariant under that symmetry.

II.2 Wannier orbitals

We shall assume in the following Wannier orbitals (WOs) that are centered at the Wyckoff positions $2c$ of the moiré triangular superlattice, i.e., at the AB and BA region centers, even though their probability distribution has a very substantial component in other regions such as AA. The large size of the unit cell has originated so far an intense effort to find a minimal model faithfully describing the FBs physics. While a variety of WOs centered at different Wyckoff positions has been proposed Rademaker and Mellado (2018); Zou et al. (2018); Po et al. (2018); Kang and Vafek (2018); Po et al. (2018); Koshino et al. (2018); Yuan and Fu (2018); Song et al. (2018); Carr et al. (2019), our simple assumption is well suited for our purposes. The site symmetry at the Wyckoff positions $2c$ is $D_{3}$, and includes only $\text{C}_{3z}$ and $\text{C}_{2x}$ with irreps $A_{1}$, $A_{2}$ and $E$. Inspired by the symmetries of the Bloch states at the high-symmetry points Angeli et al. (2018), see Sec. III), we consider two $A_{1}$ and $A_{2}$ one-dimensional irreps (1d-irreps), both invariant under $\text{C}_{3z}$ and eigenstates of $\text{C}_{2x}$ with opposite eigenvalues $c_{2x}=\pm 1$. In addition, we consider one two-dimensional irrep (2d-irrep) $E$, which transforms under $\text{C}_{3z}$ and $\text{C}_{2x}$ as the 2d-irreps in Table 1, hence comprises eigenstates of $\text{C}_{2x}$, with opposite eigenvalues $c_{2x}=\pm 1$, which are not invariant under $\text{C}_{3z}$.

We define on sublattice $\alpha=AB,BA$ the spin-$\sigma$ WO annihilation operators $\Psi_{\alpha,\mathbf{R}\sigma}$ and $\Phi_{\alpha,\mathbf{R}\sigma}$ corresponding to the two 1d-irreps ( $\Psi$) or the single 2d-irrep ($\Phi$), respectively

[TABLE]

where the subscript $s$ refers to $c_{2x}=+1$, and $p$ to $c_{2x}=-1$, while the labels 1 and 2 refer to the two valleys. It is implicit that each component is itself a spinor that includes fermionic operators corresponding to different WOs that transform like the same irrep. We shall combine the operators of different sublattices into a single spinor

[TABLE]

We further introduce three different Pauli matrices $\sigma_{a}$ that act in the moiré sublattice space (AB, BA), $\mu_{a}$ in the $c_{2x}=\pm 1$ space ($s$, $p$), and $\tau_{a}$ in the valley space (1, 2), where $a=0,1,2,3$, $a=0$ denoting the identity.

With these definitions, the generator (1) of the valley $U_{v}(1)$ symmetry becomes simply

[TABLE]

It is now worth deriving the expression of the space group symmetry operations in this notation and representation.By definition, the $\text{C}_{2x}$ transformation corresponds to the simple operator

[TABLE]

The $180^{\circ}$ rotation around the $z$-axis that connects sublattice AB with BA of each layer ($C_{2z}$) is not diagonal in the valley indices and can be represented by Song et al. (2018); Zou et al. (2018)

[TABLE]

Finally, since $\text{C}_{2y}=\text{C}_{2z}\times\text{C}_{2x}$, then

[TABLE]

III Lattice relaxation and simmetry analysis of the band structure

Since the tBLG must undergo , relative to the ideal superposition of two rigid graphene layers, a substantial lattice relaxation , whose effects have also been observed in recent experiments Yoo et al. (2018); Zhang and Tadmor (2018); Kerelsky et al. (2018); Choi et al. (2019), we performed lattice relaxations via classical molecular dynamics techniques using state-of-the-art force-fields, allowing for both in-plane and out -of-plane deformations. The details about the relaxation procedure are in Appendix A and are essentially those in Ref.Angeli et al. (2018). It is well known Uchida et al. (2014); Dai et al. (2016); Nam and Koshino (2017); Jain et al. (2017); Gargiulo and Yazyev (2018); Angeli et al. (2018); Lucignano et al. (2019); Choi and Choi (2018); Guinea and Walet (2019) that, after full relaxation, the energetically less favourable AA regions shrink while the Bernal-stacked ones, AB and BA, expand in the $(x,y)$ plane. In addition, the interlayer distance along $z$ of the AA regions increases with respect to that of the AB/BA zones, leading to significant out-of-plane buckling deformations, genuine ’corrugations’ of the graphene layers, which form protruding AA bubbles. The main effect of the layer corrugations is to enhance the bandgaps between the FBs and those above and below Angeli et al. (2018); Lucignano et al. (2019); Liu et al. (2018a), important even if partially hindered by the h-BN encapsulation of the samples during the experiments Cao et al. (2018a). Moreover, since both the AB and BA triangular domains have expanded, the initially broad crossover region between them sharpens into narrow domain walls (DWs) that merge at the AA centers in the moiré superlattice (see 1). The electronic structure shown in Fig. 2 is obtained with standard tight-binding calculations, see A for further details, with the relaxed atomic positions, using hopping amplitudes tuned so to reproduce $ab$-$initio$ calculations Trambly de Laissardiére et al. (2010). The colored circles and triangles at the ${\boldsymbol{\Gamma}}$ and $\mathbf{K}$ points, respectively, indicate the irreps that transform like the corresponding Bloch states. For instance, at ${\boldsymbol{\Gamma}}$ the FBs consist of two doublets, the lower corresponding to the irreps $A_{1}+B_{1}$, and the upper to $A_{2}+B_{2}$. Right above and below the FBs, we find at ${\boldsymbol{\Gamma}}$ two quartets, each transforming like $E_{1}+E_{2}$. At $\mathbf{K}$, the FBs are degenerate and form a quartet $E+E$. Consistently with the $D_{3}$ little group containing $\text{C}_{3z}$ and $\text{C}_{2y}$, at $\mathbf{K}$ we find either quartets, like at the FBs, made of degenerate pairs of doublets, each transforming like $E$, or doublets transforming like $A_{1}+A_{2}$, where $A_{1}$ and $A_{2}$ differ in the parity under $\text{C}_{2y}$. This overall doubling of degeneracies beyond their expected $D_{6}$ space group irreps, reflects the valley $U_{v}(1)$ symmetry of FBs that will be discussed below, and whose eventual breaking will be addressed later in this paper. We end by remarking that the so-called ’fragile’ topology Song et al. (2018); Po et al. (2018, 2018, 2018), diagnosed by the odd winding of the Wilson loop (WL) Song et al. (2018), is actually robust against lattice corrugations Liu et al. (2018a) and relaxation, as shown by panel (c) in Fig. 2.

III.1 SU(2) symmetry and accidental degeneracy along $C_{2y}$ invariant lines

We note that along all directions that are invariant under $\text{C}_{2y}$, which include the diagonals as well as all the edges of the MBZ, the electronic bands show a twofold degeneracy between Bloch states that transform differently under $\text{C}_{2y}$. These lines corresponds to the domain walls (DWs) in real space. This "accidental" degeneracy is a consequence of the interplay between $U_{v}(1)$ and $\text{C}_{2y}$ symmetries. Indeed, along $\Gamma\to K_{1,2}$ and $M\to K_{1,2}$ in the MBZ we have that:

[TABLE]

and similarly for $\Phi^{\phantom{\dagger}}_{\mathbf{k}\sigma}$. It follows that the generator of $U_{v}(1)$, i.e., the operator $\sigma_{0}\,\tau_{3}\,\mu_{0}$, anticommutes with the expression of $\text{C}_{2y}$ along the lines invariant under that same symmetry, namely the operator $\sigma_{1}\,\tau_{1}\,\mu_{3}$, and both commute with the Hamiltonian. Then also their product, $\sigma_{1}\,\tau_{2}\,\mu_{3}$, commutes with the Hamiltonian and anticommutes with the other two. The three operators

[TABLE]

thus realise an $SU(2)$ algebra and all commute with the Hamiltonian $\hat{H}_{\mathbf{k}}$ in momentum space for any $\text{C}_{2y}$-invariant $\mathbf{k}$-point. This emergent $SU(2)$ symmetry is therefore responsible of the degeneracy of eigenstates with opposite parity under $\text{C}_{2y}$. We note that $C_{2x}$ instead commutes with $U_{v}(1)$, so that there is no $SU(2)$ symmetry protection against valley splitting along $C_{2x}$ invariant lines ($\Gamma\to M$).

IV Phonons in twisted bilayer graphene

The $U_{v}(1)$ valley symmetry is an emergent one since, despite the fact that its generator (1) does not commute with the Hamiltonian, the spectrum around charge neutrality is nonetheless $U_{v}(1)$–invariant. It is therefore not obvious to envisage a mechanism that could efficiently break it.

However, since the lattice degrees of freedom play an important role at equilibrium, as discussed in Section III, it is possible that they could offer the means to destroy the $U_{v}(1)$ valley symmetry. In this Section we shall show that they indeed provide such a symmetry-breaking tool.

IV.1 Valley splitting lattice modulation: the key role of the Domain Walls

In Section II.1 we mentioned that the valley symmetry arises because, even though inequivalent Dirac nodes of the two layers should be coupled to each other by the Hamiltonian after being folded onto the same point of the MBZ (see bottom panel in Fig. 1), at small angle these matrix elements are vanishingly small and thus the valleys are effectively decoupled.

That is true if the carbon lattice, although mechanically relaxed, is unperturbed by the presence of the electrons. Once coupling with electrons is considered, we cannot exclude that for example a lattice distortion modulated with the wave vectors connecting the inequivalent Dirac nodes of the two layers, $K_{+}$ with $K^{\prime}_{-}$ and $K^{\prime}_{+}$ with $K_{-}$ in Fig. 1, might instead yield a significant matrix element among the valleys. To investigate that possibility we build an ad-hoc distortion into the bilayer carbon atom positions. We define the vector ${\mathbf{q}_{ij}}={\mathbf{K}_{+,i}}-{\mathbf{K}_{-,j}^{\prime}}$, where $i,j=1,2,3$ run over the three equivalent Dirac points of the BZ of each layer, and the $D_{6}$ conserving displacement field

[TABLE]

where $a$ runs over all atomic positions, and $\mathbf{u}_{a,ij}$ corresponds to a displacement of atom $a$ in direction $\mathbf{q}_{ij}$, whose (tuneable) magnitude is the same for all atoms. In proximity of the AA regions the distortion is locally similar to the graphene breathing mode, the in-plane transverse optical phonon at $\mathbf{K}$ Jorio et al. (2017) with $A_{1}$ symmetry. Since by construction $\mathbf{q}_{ij}$ is a multiple integer of the reciprocal lattice vectors, $\boldsymbol{\eta}(\mathbf{r})$ has the same periodicity of the unit cell, i.e., the distortion is actually at the ${\boldsymbol{\Gamma}}$ point. Moreover, the distortion’s $D_{6}$ invariance implies no change of space group symmetries.

Since the most direct evidence of the $U_{v}(1)$ symmetry is the accidental degeneracy in the band structure along all $\text{C}_{2y}$–invariant lines, corresponding just to the DWs directions in real space, we further assume the action of the displacement field $\boldsymbol{\eta}(\mathbf{r})$ to be restricted to a small region in proximity of the DWs, affecting only $\approx 1\%$ of the atoms in the moiré supercell, see top panels in Fig. 3. The modified FBs in presence of the displacement field $\boldsymbol{\eta}(\mathbf{r})$ with \big{|}\mathbf{u}_{a,ij}\big{|}=20 mÅ are shown in the bottom panel of Fig. 3. Remarkably, despite the minute distortion magnitude and the distortion involving only the minority of carbon atoms in the DWs, the degeneracy along ${\boldsymbol{\Gamma}}\to\mathbf{K}_{1}$ and $\mathbf{M}\to\mathbf{K}_{2}$ is lifted to such an extent that the four bands split into two similar copies. This is remarkable in two aspects. First, we repeat, because there is no space symmetry breaking. The only symmetry affected by the distortion is $U_{v}(1)$, since, by construction, $\boldsymbol{\eta}(\mathbf{r})$ preserve the full space group symmetries. Second, the large splitting magnitude reflects an enormous strength of the effective electron-phonon coupling, whose origin is interesting. Generally speaking, in fact, broad bands involve large hoppings and large absolute electron-phonon couplings, while the opposite is expected for narrow bands. The large electron-phonon couplings which we find for the low energy bands of tBLG suggests a possible broad-band origin of the FBs, as we shall discuss in Sec. VII. A consequence of this e-p coupling magnitude is that all potential phenomena involving lattice distortions, either static or dynamic, should be considered with a much larger priority than done so far.

IV.2 The phonon spectrum

We compute the phonon eigenmodes in the relaxed bilayer structure, by standard methods, see Appendix A for details of the calculation. Fig. 4 shows the phonon density of states $F(\omega)$ for tBLG at three different twist angles in comparison with the Bernal AB-stacked bilayer. As previously reported Choi and Choi (2018); Cocemasov et al. (2013), $F(\omega)$ is almost independent of the twist angle, which only affects the inter-layer van der Waals forces, much weaker than the in-plane ones arising from the stiff $C-C$ bonds. As a consequence, phonons in tBLG are basically those of the Bernal stacked bilayer. This is true except for a small set of special phonon modes, clearly distinguishable in Fig. 5 that depicts the phonon spectrum zoomed in a very narrow energy region $\approx 0.04$ $meV$ around the high frequency graphene $\mathbf{K}$-point peak of the phonon density of states. Specifically, within the large number of energy levels of all other highly dispersive phonon bands, unresolved in the narrow energy window, a set of 10 almost dispersionless modes emerges. We note that these special modes show the same accidental degeneracy doubling along the $C_{2y}$ invariant lines as that of the electronic bands around the charge neutrality point.

Similar to the electronic degeneracy, whose underlying $U_{v}(1)$ symmetry arises from vanishingly weak hopping matrix elements between interlayer K-K’ points in the Hamiltonian, the mechanically weak van der Waals interlayer coupling here leads to an effective $U_{v}(1)$ symmetry for this group of lattice vibrations. Their poor dispersion is connected with a displacement which is non-uniform in the supercell, and is strongly modulated on the moiré length scale, a distinctive feature of these special modes that we shall denote as ‘moiré phonons’. In particular, it is maximum in the center of the AA zones, finite in the DWs, and negligible in the large AB and BA Bernal regions. The overlap between the displacement $\boldsymbol{\eta}(\mathbf{r})$ in (10) and the $33492$ phonons of the $\theta=1.08$ tBLG at the ${\boldsymbol{\Gamma}}$ point is non-negligible only for those moiré phonons. In particular, we find the highest overlap with the doubly degenerate mode marked by a red arrow in Fig. 5, and which transforms like $A_{1}+B_{1}$. In Fig. 6a) we show the real space distortion corresponding to the $A_{1}$ component of the doublet, where the displacement direction is represented by small arrows, while its intensity is encoded in colors. This inspection of the eigenvectors of these modes at the atomistic level reveals a definite underlying single-layer graphene character, specifically that of the $A_{1}$ symmetry transverse optical mode at $\mathbf{K}$. Since the graphene $\mathbf{K}$-point does not fold into the bilayer $\Gamma$-point, their appearance along the whole $\Gamma\to\mathbf{K}$ line in the spectrum of the fully relaxed bilayer must be merely a consequence of relaxation, a relaxation that is particularly strong precisely in the AA and DW regions.

V Insulating state at charge neutrality

We will next focus on the effect on the electronic band structure of a carbon atom displacement corresponding to the two degenerate phonon modes $A_{1}$ and $B_{1}$ at $\Gamma$, which should affect the valley symmetry as the displacement in Fig. 3. In order to verify that, we carried out a frozen phonon calculation of the modified FB electronic structure with increasing intensity of the deformation see supplementary video at http://…… . Remarkably, despite transforming as different irreps ($A_{1}$ or $B_{1}$), both frozen phonon distortions are not only degenerate, but have exactly the same effect on the bands. As soon as the lattice is distorted, see Fig. 7 b), the fourfold degeneracy at $\mathbf{K}_{1}$ and $\mathbf{K}_{2}$, and the twofold one at $\Gamma$ and $\mathbf{M}$ is lifted, and small avoided crossings appear and start to move from $\mathbf{K}_{1(2)}$ towards the $\mathbf{M}$ points. Once they cross $\mathbf{M}$, they keep moving along $\mathbf{M}\to{\boldsymbol{\Gamma}}$, see Fig. 7 c). However, along these directions, the $C_{2x}$ symmetry prevents the avoided crossings, and thus leads to six elliptical Dirac cones. Finally, once they reach ${\boldsymbol{\Gamma}}$ at a threshold value of the distortion, the six Dirac points annihilate so that a gap opens at charge neutrality, see Fig. 7 d). This gap-opening mechanism is very efficient, with large splittings even for small values of the atomic displacement amplitude shown in Fig. 7. For instance, an average displacement as small as $\approx 0.5~{}\text{m\AA}$ per atom is enough to completely separate the four FBs and to open a gap at charge neutrality.

We emphasize that this occurs without breaking any spatial symmetry of the tBLG, just $U_{v}(1)$. As a consequence, the insulator state possesses a non trivial topology, as highlighted by the odd-winding of the Wilson loop of the lowest two bands, shown in the inset of the same figure 7 d). In turn, the non trivial topology of the system implies the existence of edge states within the gap separating the two lower flat bands from the two upper ones. We thus recalculated the band structure freezing the moiré phonon with $A_{1}$ symmetry at ${\boldsymbol{\Gamma}}$ in a ribbon geometry, which is obtained by cutting the tBLG along two parallel domain walls in the $y$-direction at a distance of 7 supercells. The ribbon has therefore translational symmetry along $y$, but it is confined in the $x$-direction. In Fig. 8 we show the single-particle energy levels as function of the momentum $k_{y}$ along the $y$-direction. Edge states within the gap at charge neutrality are clearly visible. In particular, we find for each edge two counter-propagating modes.

As a matter of fact, recent experiments do report the existence of a finite gap also at charge neutrality Lu et al. (2019); Jiang et al. (2019), which is actually bigger than at other non zero integer fillings, and appears without a manifest breakdown of time reversal symmetry T Lu et al. (2019) or $\text{C}_{2z}$ symmetry Jiang et al. (2019), and with the FBs still well separated from other bands Jiang et al. (2019). These evidences seem not to support an interaction driven gap, which would entail either T or $\text{T}\,\text{C}_{2z}$ symmetry breakings Lu et al. (2019); Liu et al. (2019), or else the FBs touching other bands at the ${\boldsymbol{\Gamma}}$ point Lu et al. (2019). Our phonon-driven insulator at charge neutrality breaks instead $U_{v}(1)$ and, eventually, $\text{C}_{2y}$ if the frozen phonon has $B_{1}$ character, both symmetries experimentally elusive. However, the edge states that we predict could be detectable by STM or STS, thus providing support or disproving the mechanism that we uncovered.

We end mentioning that near charge neutrality there are compelling evidences of a substantial breakdown of $\text{C}_{3z}$ symmetry Kerelsky et al. (2018); Choi et al. (2019); Jiang et al. (2019), which, although is not expected to stabilise on its own an insulating gap unless the latter were in fact just a pseudo gap Liu et al. (2019), still it is worth being properly discussed, which we postpone to Section VI.1.

V.1 $\text{E}\otimes\text{e}$ Jahn-Teller effect

The evidence that the two degenerate modes marked by a red arrow in Fig. 5 produce the same band structure in a frozen-phonon distortion is reminiscent of a $\text{E}\otimes\text{e}$ Jahn-Teller effect, i.e., the coupling of a doubly degenerate vibration with a doubly degenerate electronic state Englman (1972).

Let us consider the action of the two $\Gamma$-point $A_{1}$ and $B_{1}$ phonon modes, hereafter denoted as $q_{1}$ and $q_{2}$, along the $\text{C}_{2y}$ invariant lines. We note that the $A_{1}$ mode, $q_{1}$, although invariant under the $D_{6}$ group elements, is able to split the degeneracy along those lines. Therefore it must be coupled to the electrons through a $D_{6}$ invariant operator that does not commute with the $U_{v}(1)$ generator $\tau_{3}$. That in turn cannot but coincide with $\text{C}_{2y}$ itself, which, along the invariant lines, is the operator $T_{1}=\sigma_{1}\,\tau_{1}\,\mu_{3}/2$ in Eq. (9). On the other hand, the $B_{1}$ mode, $q_{2}$, is odd under $\text{C}_{2y}$, and thus it must be associated with an operator that anticommutes with $\text{C}_{2y}$ and does not commute with $\tau_{3}$. The only possibility that still admits a $U(1)$ valley symmetry is the operator $T_{2}=\sigma_{1}\,\tau_{2}\,\mu_{3}/2$ in Eq. (9). Indeed, with such a choice, the electron-phonon Hamiltonian is

[TABLE]

with $g$ the coupling constant. This commutes with the operator

[TABLE]

(where $\mathbf{p}=(p_{1},p_{2})$ the conjugate variable of the displacement $\mathbf{q}=(q_{1},q_{2})$), the generator of a generalized $U_{v}(1)$ symmetry that involves electron and phonon variables. As anticipated, the Hamiltonian (11) describes precisely a $\text{E}\otimes\text{e}$ Jahn-Teller problem Englman (1972); Fabrizio et al. (1996); Fabrizio and Tosatti (1997).

Since the phonon mode $\mathbf{q}$ is almost dispersionless, see Fig. 5, we can think of it as the vibration of a moiré supercell, as if the latter were a single, though very large, molecule, and the tBLG a molecular conductor. In this language, the band structures shown in Fig. 7 would correspond to a static Jahn-Teller distortion. However, since the phonon frequency is substantially larger than the width of the flat bands, we cannot exclude the possibility of a dynamical Jahn-Teller effect that could mediate superconductivity Fabrizio et al. (1996); Capone et al. (2004), or even stabilize a Jahn-Teller Mott insulator Fabrizio and Tosatti (1997) in presence of a strong enough interaction.

The Jahn-Teller nature of the electron-phonon coupling entails a very efficient mechanism to split the accidental degeneracy, linear in the displacement within a frozen phonon calculation. However, this does not explain why an in-plane displacement as small as 0.57 mÅ is able to split the formerly degenerate states at ${\boldsymbol{\Gamma}}$ by an amount as large as $15~{}\text{meV}$, see Fig. 7 d), of the same order as the original width of the flat bands. To clarify that, we note that this displacement would yield a change in the graphene nearest neighbour hopping of around $\delta t\simeq 3.6~{}\text{meV}$, see Eq. (25), which in turn entails a splitting at ${\boldsymbol{\Gamma}}$ of $6\,\delta t\sim 21~{}\text{meV}$, close to what we observe. We believe that such correspondence is not accidental, but indicates that the actual energy scale underneath the flat bands is on the order of the bare graphene bandwidth, rather than the flat bandwidth itself. We shall return on this issue later in Section VII.

VI Insulating states at other

commensurate fillings

The $\Gamma$-point distortion described above can lead to an insulating state at charge neutrality, possibly connected with the insulating state very recently reported Lu et al. (2019); Jiang et al. (2019). The same phonon branch might also stabilize insulating states at other integer occupancies of the mini bands besides charge neutrality. However, this necessarily requires freezing a mode at a high-symmetry $\mathbf{k}$-point different from ${\boldsymbol{\Gamma}}$ in order to get rid of the band touching at the Dirac points, $\mathbf{K}_{1}$ and $\mathbf{K}_{2}$, protected by the $\text{C}_{6z}$ symmetry. We consider one of the two degenerate phonon modes at the $\mathbf{M}$-point marked by a green arrow in Fig. 5. This mode has similar features as the Jahn-Teller one at ${\boldsymbol{\Gamma}}$, even if it belongs to an upper branch due to an avoided crossing along ${\boldsymbol{\Gamma}}\to\mathbf{M}$ (blue dashed line in Fig. 5). In Fig. 9 we depict this mode, which still transforms as the $A_{1}$ of graphene $\mathbf{K}$-point on the microscopic graphene scale, but whose long-wavelength modulation now forms a series of ellipses elongated along some of the DWs, thus macroscopically breaking the $C_{3z}$ symmetry of the moiré superlattice.

In Fig. 10 we show the DOS of the FBs obtained by a frozen-phonon realistic tight-binding calculation. Besides the band gaps at $\nu=\pm 4$, which separates the FBs from the others bands, now small gaps opens at $\nu=\pm 2$ with an average atomic displacement of $1.8~{}\text{m\AA}$ induced by the mode at $\mathbf{M}$.

Finally we also considered a more exotic multi-component distortion induced by a combination of the modes at the three inequivalent $\mathbf{M}$ points, which quadruples the unit cell (see Fig. 11). The resulting DOS of the FBs is shown in Fig. 12, and displays small gaps at the odd integer occupancies $\nu=1$ and $\nu=\pm 3$.

The very qualitative conclusion of this exploration is that frozen phonon distortions with various $\mathbf{k}$-vectors can very effectively yield Peierls-like insulating states at integer hole/electron fillings. Of course all distortions at $\mathbf{k}$-points different from $\Gamma$ also represent super-superlattices, with an enlargement of the unit cell that should be verifiable in such a Peierls state, even if the displacement is a tiny fraction of the equilibrium C–C distance. We note here that the the zone boundary phonons are less effective in opening gaps at non zero integer fillings than the zone center phonons are at the charge neutrality point. The reason is that away from charge neutrality the Jahn-Teller effect alone is no longer sufficient; one needs to invoke zone boundary phonons that enlarge the unit cell, and thus open gaps at the boundaries of the folded Brillouin zone. The efficiency of such gap opening mechanism is evidently lower than that of the Jahn-Teller ${\boldsymbol{\Gamma}}$-point mode.

VI.1 $\text{C}_{3z}$ symmetry breaking

We observe that the multicomponent distortion with the phonons frozen at the inequivalent $\mathbf{M}$ points leads to a width of the FBs five times larger than in the undistorted case, see Fig. 12 as opposed to Fig. 7 a), which is not simply a consequence of the tiny gaps that open at the boundary of the reduced Brillouin zone. Such substantial bandwidth increase suggests that the tBLG may be intrinsically unstable to $\text{C}_{3z}$ symmetry breaking, especially near charge neutrality. Electron-electron interaction treated in mean-field does a very similar job Choi et al. (2019). Essentially, both interaction- and phonon-driven mechanisms act right in the same manner: they move the two van Hove singularities of the fully symmetric band structure away from each other, and split them into two, with the net effect of increasing the bandwidth. As such, those two mechanisms will cooperate to drive the C3z symmetry breaking, or enhance it when explicitly broken by strain, not in disagreement with experiments Kerelsky et al. (2018); Choi et al. (2019); Jiang et al. (2019). The main difference is that the moiré phonons also break $U_{v}(1)$, and thus are able to open gaps at commensurate fillings that $\text{C}_{3z}$ symmetry breaking alone would not do.

VII Phonon mediated superconductivity

Here we indicate how the phenomena described above can connect to a superconducting state mediated by electron-phonon coupling. For that, we need to build a minimal tight-binding model containing a limited set of orbitals. For the sake of simplicity we shall not require the model Hamiltonian to reproduce precisely the shape of all the bands around charge neutrality, especially those above or below the FBs, but only the correct elementary band representation, topology and, obviously, the existence of the four flat bands separated from the all others. Considering for instance only the 32 states at ${\boldsymbol{\Gamma}}$ closest to the charge neutrality point, and maintaining our assumption of WOs centred at the Wyckoff positions $2c$, those states (apart from avoided crossings allowed by symmetry) would evolve from ${\boldsymbol{\Gamma}}$ to $\mathbf{K}_{1}$ in accordance with the $D_{6}$ space group as shown in Fig. 13. Once we allow same-symmetry Bloch states to repel each other along ${\boldsymbol{\Gamma}}\to\mathbf{K}_{1}$, the band representation can look similar to the real one (Fig. 2), including the existence of the four FBs that start at ${\boldsymbol{\Gamma}}$ as two doublets, $A_{1}+B_{1}$ and $A_{2}+B_{2}$, and end at $\mathbf{K}_{1}$ as two degenerate doublets, each transforming as the 2d-irrep $E$, see the two solid black lines in Fig. 13. While this picture looks compatible with the actual band structure, we shall take a further simplification and just consider the thicker red, blue and green bands in Fig. 13, which could still produce flat bands with the correct symmetries. This oversimplification obviously implies giving up the possibility to accurately reproduce the shape of the FBs - but it makes the algebra much simpler.

Within this approximation the components of the spinor operators in Eq. (2) are actually single fermionic operators, so that we limit ourselves to just four WOs for each sublattice, AB or BA, and valley, 1 and 2. Two of such WOs transform like the 1d-irreps, one even, $A_{1}$, and the other odd, $A_{2}$, under $\text{C}_{2x}$. The other two instead transform like the 2d-irrep $E$. We build a minimal tight-binding model

[TABLE]

where $\Delta$ splits the $s$ from the $p$ WOs of the 1d-irreps, and $T_{\text{el}}$ includes first and second neighbour hopping between AB and BA regions in the moiré superlattice compatible with all symmetries, see Appendix B for details. In Fig. 13 we show the resulting FBs as well as their Wilson loop. We emphasise that the FBs arise in this picture from a sequence of avoided crossings between a large set of relatively broad bands that strongly repel each other away from the high symmetry points, rather than from truly localized’ WOs. This mechanism is also compatible with, and in fact behind, the strong electron-phonon coupling strength, and the consequently large effects on the FBs of the Jahn-Teller phonons that, according to Eq. (11), simply modulate the AB–BA hopping.

VII.1 Mean-field superconducting state

Neglecting the extremely small dispersion of the moiré Jahn-Teller phonons, we can write their Hamiltonian simply as

[TABLE]

with $\omega\simeq 207~{}\text{meV}$.

Rather than trying to model more faithfully the Jahn-Teller coupling (11), we shall follow a simplified approach based just on symmetry considerations.

In general, we could integrate out the phonons to obtain a retarded electron-electron attraction that can mediate superconductivity. However, since here the phonon frequency is much larger that the bandwidth of the FBs, where the chemical potential lies, we can safely neglect retardation effects making a BCS-type approximation virtually exact. The attraction thus becomes instantaneous and can be represented as in Fig. 14. The phonon couples electrons in nearest neighbor AB and BA regions, giving rise to an inter-moiré site spin-singlet pairing, a state which we expect to be much less affected by Coulomb repulsion than an on-site one. Therefore, neglecting Coulomb repulsions we can concentrate on the pairing channel between nearest neighbor AB and BA regions. The scattering processes in Fig. 14 imply that the pairing channels are only $\tau_{1}\,\mu_{0}$ and $\tau_{1}\,\mu_{3}$, corresponding to inter valley pairing, as expected because time reversal interchanges the two valleys.

Having assumed pairing between nearest neighbor AB and BA regions, we must identify pair functions in momentum space that connect nearest neighbor unit cells, and transform properly under $\text{C}_{3z}$.

These functions are

[TABLE]

where $\omega=\text{e}^{i2\pi/3}$. Specifically, $\gamma(\mathbf{k})\sim A_{1}$ is invariant under $\text{C}_{3z}$, while

[TABLE]

In other words, \big{(}\gamma_{+1}(\mathbf{k}),\gamma_{-1}(\mathbf{k})\big{)} form a representation of the 2d-irrep E=\big{(}E_{+1},E_{-1}\big{)} in which $\text{C}_{3z}$ is diagonal with eigenvalues $\omega$ and $\omega^{*}$.

Here it is more convenient to transform the spinor $\Phi_{\mathbf{k}\sigma}$

[TABLE]

so that $\Phi^{\phantom{\dagger}}_{\pm 1,\mathbf{k}\sigma}$ is associated with a WO that transforms like $E_{\pm 1}$. Under the assumption of pairing diagonal in the irreps, we can construct the following spin-singlet Cooper pairs:

[TABLE]

which can be combined with the $\mathbf{k}$-dependent functions in (15) to give pair operators that transform like the irreps of $D_{3}$. For instance, multiplying (18) by $\gamma_{+1}(\mathbf{k})$, (19) by $\gamma_{-1}(\mathbf{k})$, (21) with the plus sign by $\gamma(\mathbf{k})$, or (22) with $\mu_{0}$ by $\gamma(\mathbf{k})$, we obtain pair operators that all transform like $A_{1}$. We shall denote their sum as $\textbf{A}^{\dagger}_{1\mathbf{k}}$, and, similarly, all other symmetry combinations as $\textbf{A}^{\dagger}_{2\mathbf{k}}$ and $\textbf{E}^{\dagger}_{\pm 1,\mathbf{k}}$. Evidently, since in our modeling the FBs are made of 1d and 2d irreps, the gap function will in general involve a combination of $\gamma(\mathbf{k})$ and $\gamma_{\pm 1}(\mathbf{k})$, namely, it will be a superposition of $s$ and $d\pm id$ symmetry channels. The dominant channel will depend on which is the prevailing WO character of the Bloch states at the chemical potential, as well as on the strength of the scattering amplitudes in the different pairing channels, $\textbf{A}^{\dagger}_{1\mathbf{k}}$, $\textbf{A}^{\dagger}_{2\mathbf{k}}$ and $\textbf{E}^{\dagger}_{\pm 1,\mathbf{k}}$.

In general, we may expect the totally-symmetric $A_{1}$ channel to have the largest amplitude, thus we assume the following expression of the phonon-mediated attraction

[TABLE]

that involves a single parameter $\lambda\sim g^{2}/2\omega$. We treat the full Hamiltonian (13) plus (23) in mean field allowing for a superconducting solution, which is always stabilized by the attraction provided the density of states is finite at the chemical potential. We find that superconductivity opens a gap everywhere in the Brillouin zone. Since

[TABLE]

the order parameter may have finite components with different symmetries, $E_{\pm 1}$ and $A_{1}$. In the model calculation all components acquire similar magnitude, implying a mixture of $s$ and $d\pm id$ wave symmetries. In Fig. 15 we show $\langle\,A^{\prime}_{1,\mathbf{k}}\,\rangle$ and $\langle\,E_{-1,\mathbf{k}}\,\rangle$ at the Fermi surface corresponding to densities $\nu\approx-1$ and $\nu\approx-2$ with respect to charge neutrality. We conclude by emphasizing that the Cooper pair is made by one electron in AB and one in BA, thus leading, in the spin-singlet channel, to extended $s$ and/or $d\pm id$ symmetries. That is merely a consequence of the phonon mode and electron-phonon properties, hence it does not depend on the above modeling of the FBs.

There are already in literature several proposals about the superconducting states in tBLG. Most of them, however, invoke electron correlations as the element responsible, or strongly effecting the pairing Xu and Balents (2018); Po et al. (2018); Baskaran (2018); Dodaro et al. (2018); Isobe et al. (2018); You and Vishwanath (2018); Kennes et al. (2018); Liu et al. (2018b); Zhang et al. (2019); Fidrysiak et al. (2018); Laksono et al. (2018); Su and Lin (2018); Gu et al. (2019); González and Stauber (2019, 2019); Roy and Juričić (2019); Classen et al. (2019); Guo et al. (2018); Huang et al. (2018); Lin and Nandkishore (2019). There are few exceptions Wu et al. (2018); Lian et al. (2018); Wu et al. (2018) that instead propose, as we do, a purely phonon mediated attraction. In both works Wu et al. (2018) and Lian et al. (2018) the tBLG phonons are assumed to coincide with the single layer graphene ones, as if the interlayer coupling were ineffective in the phonon spectrum, which is not what we find for the special modes discussed above. Moreover, they both discuss the effects of such phonons only in a continuum model for the FBs. In particular, the authors of Wu et al. (2018) consider a few selected graphene modes, among which the transverse optical mode at $\mathbf{K}$ that has the largest weight in the tBLG phonon that we consider. They conclude that such graphene mode mediates $d\pm id$ pairing in the $A_{2}$ channel, $\tau_{1}\,\mu_{3}$ in our language, leading to an order parameter odd upon interchanging the two layers. On the contrary, the authors of Lian et al. (2018) focus just on the acoustic phonons of graphene, and conclude they stabilize an extended $s$-wave order parameter.

VIII Conclusions

In this work, we uncovered a novel and strong electron-phonon coupling mechanism and analysed its potential role in the low temperature physics of magic angle twisted bilayer graphene, with particular emphasis on the insulating and superconducting states. By working out the phonon modes of a fully relaxed tBLG, we found a special group of few of them modulated over the whole moiré supercell and nearly dispersionless, thus showing that the moiré pattern can induce flat bands also in the phonon spectrum. In particular, two of these modes, which are degenerate at any $\mathbf{k}$-point invariant under $180^{\circ}$ rotation around the $y$-axis and its $\text{C}_{3z}$ equivalent directions, are found to couple strongly to the valley $U_{v}(1)$ symmetry, that is responsible for the accidental degeneracy of the band structure at the same $\mathbf{k}$-points where the two modes are degenerate. This particular phonon doublet is strongly Jahn-Teller coupled to the valley degrees of freedom, realizing a so-called $\text{E}\otimes\text{e}$ Jahn-Teller model. This mechanism, if static, would generate a filling-dependent broadening of the flat bands and eventually insulating phases at all the commensurate fillings. Interestingly, freezing the modes at ${\boldsymbol{\Gamma}}$ can stabilize a topologically non-trivial insulator at charge neutrality that sustains edge modes. We also investigate the symmetry properties of a hypothetical superconducting state stabilized by this Jahn-Teller mode. We find that the phonon mediated coupling occurs on the moiré scale, favoring spin-singlet pairing of electrons in different Bernal (AB/BA) regions, which may thus condense with an extended $s$- and/or $d\pm id$-wave order parameter. In a mean-field calculation with a model tight-binding Hamiltonian of the flat bands, we find that the dominant symmetry depends on the orbital character that prevails in the Bloch states at the Fermi energy, as well as on the precise values of scattering amplitudes in the different Cooper channels allowed by the $\text{C}_{3z}$ symmetry. We cannot exclude that a nematic component might arise due to higher order terms not included in our mean-field calculation, as discussed in Sigrist and Ueda (1991) and Kozii et al. (2018). These results herald a role of phonons and of lattice distortions of much larger impact than supposed so far in twisted graphene bilayers, which does not exclude a joint action of electron-phonon and electron-electron interaction. Further experimental and theoretical developments will be called for in order to establish their actual importance and role.

Acknowledgments

The authors are extremely grateful to D. Mandelli for the technical support provided during the lattice relaxation procedure. We acknowledge useful discussions with P. Lucignano, A. Valli, A. Amaricci, M. Capone, E. Kucukbenli, A. Dal Corso and S. de Gironcoli. M. F. acknowledges funding by the European Research Council (ERC) under H2020 Advanced Grant No. 692670 “FIRSTORM”. E. T. acknowledges funding from the European Research Council (ERC) under FP7 Advanced Grant No.320796 “MODPHYSFRICT” , later continuing under Horizon 2020 Advanced Grant No. 824402 ”ULTRADISS”.

Appendix A Details on the lattice relaxation, band structure and phonon calculations

The lattice relaxation and bandstructure calculation procedures are the same as those thoroughly described in Ref. Angeli et al. (2018), with the exception that now the carbon-carbon intralayer interactions are modeled via the Tersoff potential Lindsay and Broido (2010). The interlayer interactions are modelled via the Kolmogorov-Crespi (KC) potential Kolmogorov and Crespi (2005), using the recent parametrization of Ref. Ouyang et al., 2018. Geometric optimizations are performed using the FIRE algorithm Bitzek et al. (2006). The hopping amplitudes of the tight-binding Hamiltonian are:

[TABLE]

where $\mathbf{d}=\mathbf{r}_{i}-\mathbf{r}_{j}$ is the distance between atom $i$ and $j$, $d=|\mathbf{d}|$, and $\textbf{e}_{z}$ is the unit vector in the direction perpendicular to the graphene planes. The out-of-plane ($\sigma$) and in-plane ($\pi$) transfer integrals are:

[TABLE]

where $V_{pp\sigma}^{0}=0.48\;eV$ and $V_{pp\pi}^{0}=-2.8\;eV$ are values chosen to reproduce ab-initio dispersion curves in AA and AB stacked bilayer graphene, $d_{0}=3.344\text{\AA}$ is the starting inter-layer distance, $a_{0}=1.3978\text{\AA}$ is the intralayer carbon-carbon distance obtained with the Tersoff potential, and $r_{0}=0.3187~{}a_{0}$ is the decay length Trambly de Laissardiére et al. (2010); Nam and Koshino (2017).

We compute phonons in tBLG using the force constants of the non-harmonic potentials $U=U_{TERSOFF}+U_{KC}$ that we used to relax the structure:

[TABLE]

where ($i,j$) label the atoms in the unit cell, ($l,s$) the moiré lattice vectors and $\alpha,\beta=x,y,z$. Then, we define the Dynamical Matrix at phonon momentum $\mathbf{q}$ as:

[TABLE]

where $M_{C}$ is the carbon atom mass. Using the above relation we determine the eigenvalue equation for the normal modes of the system and the phonon spectrum:

[TABLE]

where $\omega_{\mathbf{q}}$ denote the energy of the normal mode $\epsilon(\mathbf{q})$. The phonon dispersion of Bernal stacked bilayer graphene obtained with this method is shown Fig. 16. As can be seen, the transverse optical (TO) modes at $K$ are found at $\omega\approx 207$ meV. As a consequence, the Jahn-Teller modes discussed in the main text, which vibrate in the same way on the graphene scale, have similar frequency. However, the frequency of the TO modes in graphene is strongly sensitive to the choice of the intralayer potential used Koukaras et al. (2015), so that these modes can be predicted to have frequencies as low as $\approx 170$ meV Nika and Balandin (2012); Basko and Aleiner (2008). This implies that also the JT modes may be observed at a lower frequencies than ours.

Appendix B Hamiltonian in momentum space

In order to make the invariance under $\text{C}_{3z}$ more explicit, we shall use the transformed spinors

[TABLE]

for the 2d-irreps, which correspond to WOs eigenstates of $\text{C}_{3z}$. Moreover, it is convenient to transform also the spinors $\Psi_{\mathbf{k}\sigma}$ of the 1d-irreps in the same way, i.e.,

[TABLE]

which correspond to WOs still invariant under $\text{C}_{3z}$ but not under $\text{C}_{2x}$, whose representation both in $\Psi_{\mathbf{k}\sigma}$ and $\Phi_{\mathbf{k}\sigma}$ becomes the Pauli matrix $\mu_{1}$. In conclusion the spinor operators defined above satisfy

[TABLE]

We shall, for simplicity, consider only nearest and next nearest neighbor hopping between AB and BA region, which correspond to the following functions in momentum space:

[TABLE]

for first neighbors, and

[TABLE]

for second neighbors, where $\omega=\text{e}^{i2\pi/3}$, $\alpha_{\mathbf{k}}=\text{e}^{i\mathbf{k}\cdot(\mathbf{a}+\mathbf{b})/3}$, and the lattice constants $\mathbf{a}=(\sqrt{3}/2,-1/2)$ and $\mathbf{b}=(\sqrt{3}/2,1/2)$. Since

[TABLE]

then, for $n=1,2$,

[TABLE]

which shows that $\gamma_{n,\pm 1}(\mathbf{k})$ transform like the 2d-irrep $E$.

We assume the following tight-binding Hamiltonian for the 1d-irreps

[TABLE]

where $t^{(1)}_{11}$ and $t^{(2)}_{11}$ are the first and second neighbor hopping amplitudes, respectively, which we assume to be real.

The 2d-irreps have instead the Hamiltonian

[TABLE]

with real hopping amplitudes, where

[TABLE]

Finally, the coupling between 1d and 2d irreps is represented by the Hamiltonian

[TABLE]

with real $t^{(n)}_{12}$.

The Hamiltonian thus reads

[TABLE]

which, through the equations (32) and (36), can be readily shown to be invariant under $\text{C}_{3z}$ and $\text{C}_{2x}$, and is evidently also invariant under the $U_{v}(1)$ generator $\tau_{3}$. In addition, the Hamiltonian must be also invariant under $\text{T}\text{C}_{2z}$, where T is the time reversal operator. Noting that

[TABLE]

one can show that $H$ in (41) is also invariant under that symmetry.

The model Hamiltonian thus depends on eight parameters. The FBs shown in Fig. 13 have been obtained choosing: $\Delta=10$ , $t_{11}^{1}=2$ , $t_{22}^{1}=5$ , $g_{22}^{1}=10$ , $t_{22}^{2}=g_{22}^{2}=-t_{11}^{2}=1.2$ , $t_{12}^{1}=2$ and $t_{12}^{2}=0.5$.