Jahn–Teller coupling to moiré phonons in the continuum model formalism for small-angle twisted bilayer graphene

We show how to include the Jahn–Teller coupling of moiré phonons to the electrons in the continuum model formalism which describes small-angle twisted bilayer graphene. These phonons, which strongly couple to the valley degree of freedom, are able to open gaps at most integer fillings of the four flat bands around the charge neutrality point. Moreover, we derive the full quantum mechanical expression of the electron–phonon Hamiltonian, which may allow accessing phenomena such as the phonon-mediated superconductivity and the dynamical Jahn–Teller effect.


Introduction
The discovery of superconductivity first in small-angle twisted bilayer graphene (tBLG) [1][2][3][4], and later in trilayer [5] and double bilayer graphene [6], has stimulated an intense theoretical and experimental research activity. In these systems, the twist angle tunes a peculiar interference within a large set of energy bands, compressing energy levels to form a set of extremely narrow bands around charge neutrality [7][8][9][10]. In twisted bilayer graphene, these flat bands (FBs) have a bandwidth of the order of ≈ 10 − 20 meV and are isolated in energy by single-particle band gaps of the order of 30-50 meV. Superconductivity is observed upon doping such narrow bands, often surrounding insulating states at fractional fillings that contradict the metallic behaviour predicted by band structure calculations. The observed phenomenology of these insulating states, which turn metallic above a threshold Zeeman splitting or above a critical temperature, suggests that they might arise from a weak-coupling Stoner or CDW band instability driven by electron-electron and/or electron-phonon interactions, rather than from the Mott's localisation phenomenon in the presence of strong correlations. This is further supported by noting that the effective on-site Coulomb repulsion U must be identified with the charging energy of the supercell, which can be as large as tens of nanometers at small angles, projected onto the flat bands. If screening effects due to the gates and to the other bands are taken into account, the actual value of U is estimated of the order of few meVs, suggesting that tBLG might not be more correlated than a single graphene sheet [11]. On the contrary, there are evidences that the coupling to the lattice is instead anomalously large if compared with the FBs bandwidth. For instance, ab initio DFT-based calculations fail to predict well-defined FBs separated from other bands, unless atomic positions are allowed to relax [7,9,[12][13][14][15], in which case gaps open that are larger than the FBs bandwidth. Further evidences supporting a sizeable electron-phonon coupling come from transport properties [16][17][18], but also from direct electronic structure calculations.
Specifically, in Ref. [8] it was shown that a pair of optical phonon modes are rather strongly coupled to the FBs and thus might play an important role in the physics of tBLG. These modes, which have a long wavelength modulation on the same moiré scale, have been dubbed as "moiré phonons" and recently observed experimentally [19]. However, the large number of atoms contained in the small-angle unit cell of twisted bilayer graphene (more than 11000) makes any calculation more involved than a simple tight-binding one, rather tough, if not computationally impossible. In this paper, we try to cope with such problem by implementing the effect of these phonon modes on the band structure in the less computationally demanding continuum model of Ref. [20]. This method can serve as a suitable starting point for BCS [21,22], Hartree-Fock [23][24][25][26] and many other calculations, which may involve both phonons and correlations. The work is organised as follows. In Sect. 2, we derive the Bistritzer and MacDonald continuum model for twisted bilayer graphene. In Sect. 3, we implement in the continuum model formalism the effect of a static atomic displacement. By using lattice deformation fields which are similar to the Jahn-Teller moiré phonon modes of Ref. [8], we show how the band structure and the density of states of the system evolve as a function of the lattice distortion intensity. Finally, 4 is devoted to concluding remarks.

Continuum model Hamiltonian
We start by introducing the Bistritzer-MacDonald continuum model for twisted bilayer graphene [20] and recall that the single-layer Dirac Hamiltonian is, around the K and K = −K valleys, respectively, with ζ = ±1 the valley index, and where the Pauli matrices σ a act on the two component wavefunctions, each component referring to one of the two sites per unit cell that form the honeycomb lattice, which we shall hereafter denote as sublattices A and B. Our analysis must start by defining the specific twisted bilayer graphene we shall investigate, and by setting some conventional notations. We assume that the twisted bilayer is obtained by rotating two AA-stacked layers (i = 1, 2) by an opposite angle ±φ with respect to the centre of two overlapping hexagons, see Fig. 1a, where tan φ = 1/ (2k + 1) √ 3 , with large positive integer k. With this choice, the moiré pattern forms a superlattice, which is still honeycomb and endowed by a D 6 space group symmetry [7] that is generated by the threefold rotation C 3z around the z-axis perpendicular to the bilayer and the twofold rotations around the in-plane x-and y-axes, C 2x and C 2y , respectively.
The Dirac nodes of each monolayer are, correspondingly, R ±φ (K) ≡ K ±φ for the valley we shall denote as ζ = +1, and R ±φ (K ) ≡ K ±φ for the other valley, ζ = −1. With our choices, K +φ and K −φ fold into the same point K 2 of the MBZ, as well as K −φ and K +φ into the point K 1 , see Fig. 1.
We introduce the (real) Wannier functions derived by the p z orbital of each carbon atom: where d ⊥ = (0, 0, d), with d the interlayer distance, R (i) label the positions of the unit cells in layer i = 1, 2, while r (i) α the coordinates with respect to R (i) of the two sites within each unit cell, α = A, B denoting the two sublattices. From the Wannier functions, we build the Bloch functions Conventionally, one assumes the two-centre approximation [20], so that, if V ⊥ (r) is the interlayer potential, then the interlayer hopping dr φ 1 αR (2) depends only on the distance between the centres of the two Wannier orbitals. We define T ⊥ (q), the Fourier transform of T ⊥ (r): where r and q are vectors in the x-y plane. Hereafter, all momenta are assumed also to lie in the x-y plane. The interlayer hopping between an electron in layer 1 with momentum p and one in layer 2 with momentum k is in general a matrixT kp , with elements T αβ kp , α, β = A, B, which, through equations (5), (6) and (7), reads explicitly (2) , p+G (1) Since we are interested in the low energy physics, k and p must be close to the corresponding Dirac points, namely K φ and K φ for p in layer 1, while K −φ and K −φ for k in layer 2.
Therefore,T kp can in principle couple to each other states of different layers within the same valley or between opposite valleys. Since T ⊥ (q) decays exponentially with q = |q| [20], the leading terms are those with the least possible k + G (2) compatible with momentum conservation k + G (2) = p + G (1) . At small twist angle φ, only the intra-valley matrix elements, p ∼ k, are sizeable, while the inter-valley ones are negligibly small, despite opposite valleys of different layers fold into the same point of the MBZ. For instance, if p K +φ and k K −φ , momentum conservation requires very large 1 , thus an exponentially small T ⊥ k + G (2) . The effective decoupling between the two valleys implies that the number of electrons within each valley is to high accuracy a conserved quantity, thus an emergent valley U v (1) symmetry [20,33] that causes accidental band degeneracies along high-symmetry lines in the MBZ [8,33].
We can therefore just consider the intra-valley inter-layer scattering processes. We start with valley ζ = +1 and thus require that k is close to K −φ = K 1 and p close to K +φ = K 2 , see Fig. 1d. Since the modulus of k ∼ K 1 is invariant under C 3z rotations, where C 3z 2 , maximisation of T ⊥ k + G (2) compatibly with momentum conservation leads to the following conditions, see Eq. (3), Upon defining T (k) ≡ t ⊥ , and using Eq. (9) to evaluate the phase factors in (8), we finally obtainT where we explicitly indicate the valley index ζ , and with ω = e 2πi/3 . We now focus on the other valley, ζ = −1, and take k close to K −φ = −K 2 , and p to K φ = −K 1 , see Fig. 1d. In this case, Eq. (9) is replaced by Let us briefly discuss how one can take into account lattice relaxation, which is known to shrink the energetically unfavourable AA regions enlarging the Bernal-stacked triangular domains in the moiré pattern [7,12,14,15,27,28,30,31]. As a consequence, the inter-and intra-sublattice hopping processes acquire different amplitudes, which is taken into account by modifying the operatorsT i in Eq. (11) according tô with u generally smaller than u . We conclude by showing how this formalism allows recovering the untwisted case, where G (1) which is what one would expect from an AA-stacked bilayer.

A more convenient representation
For our purposes, it is actually more convenient to use the alternative representation of the Hamiltonian derived in Ref. [32]. We translate K φ = −K 1 so that it falls on K −φ = K 1 , and similarly, K −φ = −K 2 on K φ = K 2 , see Fig. 1d. This implies that the diagonal parts of the Following Ref. [32], we define a set of vectors which span the vertices of the MBZs, where Q A , black circles in Fig. 2, correspond to valley ζ = +1 in layer 2 and valley ζ = −1 in layer 1, while Q B , red circles in Fig. 2, correspond to valley ζ = +1 in layer 1 and valley ζ = −1 in layer 2. In addition, we define Next, we redefine the momenta for layers 1 and 2 as, respectively, thus so that the selection rules transform into With those definitions, and denoting the conserved momentum k as k, the Hamiltonian of valley ζ now readsĤ whereT In particular,T One can further simplify the notation introducing the Pauli matrices τ a , a = 0, x, y, z, with τ 0 the identity, that act in the valley subspace, and thus writê whereT i (u, u ) ≡T +1 i (u, u ), and Ω is the real unitary operator being P ζ the projector onto valley ζ , which actually interchanges sublattice A with B in the valley ζ = −1. Applying the unitary operator Ω, we thus obtain which has the advantage of having a very compact form. For convenience, we list the action of Ω applied to σ and τ operators, In this representation, any symmetry operation G ∈ D 6 corresponds to a transformation whose explicit expressions are given in Ref. [32].

Perturbation induced by a static atomic displacement
We now move to derive in the continuum model the expression of the perturbation induced by a collective atomic displacement. Under a generic lattice deformation, the in-plane atomic positions x iα change according to where i is now labelling a generic unit cell position. Since the phonon modes we are going to study involve only in-plane atomic displacements, we assume that z-coordinate of each carbon atom does not vary. It follows that a generic potential in the two-centre approximation and at linear order in the displacement reads We further neglect the dependence on z, which we will take into account by distinguishing at the end between different scattering channels, intra-and inter-layers, so that: namely assuming, as before, that T (q) depends only on q = |q|.

Jahn-Teller moiré phonon modes
Reference [8] pointed out the existence of a pair of high-frequency optical modes at the Γ point of the MBZ, which are extremely efficient in lifting the valley degeneracies observed in the band structure. These phonon modes are schematically drawn in Fig.3a, and they both share the same modulation on the moiré length scale. However, microscopically, they both look as the well-known in-plane optical phonon modes of graphene at K, which transform as the A 1 and B 1 irreducible representations, see Fig.3b. These two irreducible representations differ by the fact that B 1 is odd with respect to C 2z and C 2y , while A 1 is even with respect to all symmetries of the D 6 space group. Although the complexity of these modes is hard to capture by a simple analytical expression, their effect on the band structure can be well approximated introducing the following deformation fields (a) (b) (c) Fig. 3 a Real space phonon amplitude of the pair of Jahn-Teller phonon modes found in Ref. [8]. These modes, often referred to as "moiré modes", vibrate mostly within the AA regions involving no atomic movement at all within the Bernal (AB/BA)-stacked regions. b Sketch of the microscopic vibration of these modes. The atomic movement is in-plane and behaves as two irreducible representations (A 1 and B 1 ) of the group D 6 . c Q-vectors which connect inequivalent valleys in different layers. These vectors has been used to approximate the vibration of the modes in a and b, see (34) where Q i j are the k-vectors connecting different valleys and depicted in Fig. 3c, while a = 1, 2 is the layer index. Since the transformation C g (g = 3z, 2x) is a symmetry operation even in the distorted lattice, we have that By noting that the set of momenta {Q i j } is invariant under C 3z and C 2x , it immediately follows that, for g = 3z, 2x On the contrary, {Q i j } is not invariant under C g (g = 2y, 2z), and the phonon modes are either even (A 1 ) or odd (B 1 ) under these symmetries. Therefore, recalling that C 2z exchanges the two sublattices, If we choose u (a) , as A 1 . They both can be shortly written as where u (a) Q i j * ≡ u (a) − Q i j and u (a) Q i j is real for the B 1 distortion and imaginary for A 1 .
We end by pointing out that W (q) satisfies for all symmetry operations of the lattice, in particular

Phonon-induced Hamiltonian matrix elements
A lattice distortion involving the A 1 or B 1 phonons generates a matrix element between layer a momentum k ∼ K (a) and layer b momentum p ∼ −K (b) , where we recall that K (1) = K 2 and K (2) = K 1 in Fig. 1d: We can readily follow the same steps outlined in Sect. 2 to identify the G (a) and G (b) reciprocal lattice vectors that enforce momentum conservation and maximise the matrix element W (−q) = W (q). Therefore, we shall not repeat that calculation and jump directly to the results.
The lattice distortion introduces a perturbation both intra-layer and inter-layer. The former, in the representation introduced in Sect. 2.1, has the extremely simple expression: where τ x refers to the A 1 mode, τ y to the B 1 one, and the matricesT i (g, g ) have the same expression as those in Eq. (14), with u and u replaced, respectively, by g and g .
The inter-layer coupling has a simpler expression, since, as we many times mentioned, the opposite valleys in different layers fold on the same momentum in the MBZ, and thus, the coupling is diagonal in Q and Q and reads As before, τ x and τ y refers to the A 1 and B 1 modes, respectively. It is worth remarking that, because of the transformation (26), which exchanges the sublattices in the valley ζ = −1, the diagonal elements of the matricesT i (g, g ) in (43) and σ 0 in (44) refer to the opposite sublattices, while the diagonal elements to the same sublattice, right the opposite of the unperturbed Hamiltonian (27).
Let us rephrase the above results in the second quantisation and introducing the quantum mechanical character of the phonon mode. In the continuum model, a plane wave with momentum k + G, where G = n A + m B is a reciprocal lattice vector of the MBZ, in layer i = 1, 2, valley ζ = +1 and with sublattice components described by a two-component spinor χ k+G , can be associated with a two-component spinor operator according to For any G, we can write, see Eq. (17), and thus define We note that K +φ + K −φ ≡ G k = (2k + 1) A + B , which allows us defining the operators in valley ζ = −1 as where, in accordance with our transformation in Eq. (26), we interchange the two sublattices in valley ζ = −1 through σ x . We note that the mismatch momentum G k is just what is provided by the phonon modes. Absorbing the valley index into two additional components of the spinors, and introducing back the spin label, the second quantised Hamiltonian can be written in terms of four-component spinor operators Ψ kσ,Q , where, Q = Q A refer to layer 2 if the valley index ζ = +1 and layer 1 if ζ = −1, while Q = Q B to layer 1 if ζ = +1, and layer 2 if ζ = −1.
Next, we introduce a two-component dimensionless variable q 0 = (q 1 , q 2 ), and its conjugate one, p 0 = ( p 1 , p 2 ), where q 1 and q 2 are the phonon coordinates of the A 1 and B 1 modes at Γ = 0, respectively. Using the above-defined operators, the full quantum mechanical Hamiltonian reads where ω 0 is the phonon frequency, equal for both A 1 and B 1 modes, τ = (τ x , τ y ), and We observe that the Hamiltonian (49) still possesses a valley U v (1) symmetry, with generator where T z is half the difference between the number of electrons in valley ζ = +1 and the one in valley ζ = −1, while L z is the angular momentum of the phonon mode. The Hamiltonian (49) actually realises a e ⊗ E Jahn-Teller model. It is straightforward to generalise the above result to an atomic displacement modulated with the wave vectors Q i j + P, where P ∈ MBZ. Since Q i j are multiples of the MBZ reciprocal lattice vector, such displacement is at momentum P and can be considered as the previous one at Γ , as shown in Fig. 3, on top of which we add an additional incommensurate long wavelength component. Since P is tiny as compared to the vectors Q i j , we shall assume that the displacement has the same expression of Eq. (39), with the only difference that The full quantum mechanical Hamiltonian becomes where δĤ QQ ( P) is the same as δĤ QQ in Eq. (50) with P-dependent constants g P , g P and γ P , invariant under the little group at P. In this general case, the generator of U v (1) reads

Frozen phonon band structure
We can perform a frozen phonon calculation neglecting the phonon energy, last term in Eq. (49), and fixing q = (q 1 , q 2 ) to some constant value. Because of the U v (1) symmetry, what matters is just the modulus q of q. In practice, we have taken q = (1, 0) and studied (a) (b) (c) (d) Fig. 4 Band structure of 2φ = 1.08 • twisted bilayer graphene with increasing frozen phonon deformation intensity g. The other parameters used in this calculation are: u = 0.0761, u = 0.1031, g = g/10 and γ = g/2.5. a Undistorted case. b Slightly distorted lattice. Two small avoided crossings occur along the K 1 → Γ and M → Γ lines. c By further increasing the distortion intensity, the four bands further separate. The avoided crossing along K 1 → Γ persists, while a genuine crossing occurs along Γ → M. d A gap at charge neutrality has finally opened the band structure varying the coupling constants g, setting g = g/10 and γ = g/2.5, and assuming the following parameters:hv/a 0 = 2.1354 eV [29]; u = 0.0761 eV and u = 0.1031 eV [14]. This choice fits well the microscopic tight-binding calculations in Ref. [8]. As shown in Fig. 4, as soon as the frozen phonon terms are turned on, all the degeneracies in the band structure arising due to the valley symmetry are lifted. This occurs with a set of avoided crossings which move from K → M and from K → Γ (Fig. 4b).
In particular, the crossings that move from K → M eventually meet at M, forming (six) Dirac nodes, which then move towards Γ (Fig. 4c). Finally, at a threshold value of g, a gap opens at the charge neutrality point (Fig. 4d). Such gap keeps increasing as the deformation amplitude increases.

Moiré phonons at M
The phonon modes considered in the previous section were at position Γ of the MBZ, thus preserving the periodicity of the moiré superlattice. As pointed out in Ref. [8] and shown before, these modes are able to open a gap in the band structure only at charge neutrality. Gap opening at different commensurate fillings requires freezing finite momentum phonons [8]. Here, we consider a multicomponent distortion which involves the modes at the three inequivalent M points in the MBZ: Freezing a multiple distortion at all these points reduces by a quarter the Brillouin zone, see Fig. 5, which has now the reciprocal lattice vectors Since M i , i = 1, 2, 3, are tiny as compared to the vectors Q i j introduced in the previous section, we can make the same assumption (52) that leads to the Hamiltonian (53); namely, assume that the displacement induced by the multiple distortion has the same expression of Eq. (39), with the only difference that Both p i and q i , i = 1, 2, 3, are shown in Fig. 5. Considering all momenta k within the new BZ, the light blue hexagon in Fig. 5, and assuming that, besides the multicomponent distortion at M, there is still a distortion at Γ , the Hamiltonian can be written again as a matrixĤ QQ (k), which now readŝ where the matricesT i (x, x ) are those in Eq. (14), though they depend on different sets of parameters, (u, u ), (g, g ) and (a, a ). The crucial difference with respect to the Hamiltonian (50) with only the Γ -distortion is that the Q vectors span now the sites of the honeycomb lattice generated by the new fourfold-smaller Brillouin zone; hence, they are defined through (a) (b) (c) Fig. 6 Density of states around charge neutrality obtained with the Hamiltonian (60) with g = g/10, γ = g/2.5, a = a/10 and α = a/2.5. Gaps are highlighted in red, and the corresponding filling factor is ν. a) The undistorted lattice density of states. b) The density of states obtained by deforming the lattice with only the distortion at Γ , which opens a gap at charge neutrality. c) Density of states obtained using both the Γ and the M multicomponent distortions. Here, a gap opens at filling of one electron (three holes with respect to neutrality) per unit cell and shown in Fig. 5 as black and red circles, respectively, and must not be confused with those in Eq. (17). In Fig. 6, we show the density of states around neutrality of the Hamiltonian (60). The first two cases correspond to undistorted and Γ -only distorted structures, while the third panel involves also the M multicomponent distortion. As can be seen, a gap now opens at the partial filling of one electron per unit cell. As it was shown in Ref. [8], other phonons or combinations of them can open gaps at any integer filling of the four electronic flat bands.

Conclusions
We have shown that the moiré phonons of Ref. [8], which are coupled to the valley degrees of freedom of the electrons so to realise an E ⊗ e Jahn-Teller model, can be successfully implemented in the continuum model formalism of small-angle twisted bilayer graphene. This method is more manageable than the realistic tight-binding modelling of Ref. [8], whose results have been here reproduced with much less effort. In addition, the continuum model formalism has the great advantage of providing a full quantum mechanical expression of the electron-phonon Hamiltonian, which may allow going beyond the simple frozen phonon calculation of [8], and thus describing phenomena like a dynamical Jahn-Teller effect and the phonon-mediated superconductivity.