# Charge-transfer insulation in twisted bilayer graphene

###### Abstract

We studied the real space structure of states in twisted bilayer graphene at the ‘magic angle’ . We found that the flat bands close to charge neutrality are composed of a mix of ‘ring’ and ‘center’ orbitals around the AA stacking region. Long-range Coulomb interaction causes a charge-transfer at half-filling of the flat bands from the ‘center’ to the ‘ring’ orbitals. Consequently, the Mott phase is a featureless spin-singlet paramagnet. The superconducting state depends on the nature of the dopants: hole-doping yields -wave whereas electron-doping yields -wave pairing symmetry.

###### pacs:

Upon the discovery of superconductivity in twisted bilayer graphene (TBG)Cao et al. (2018a, b), the theoretical community jumped on the ‘flat band’-wagon and started the hunt for a simple yet sufficient model that might shed light on the nature and symmetries of this systemEfimkin and MacDonald (2018); Pal et al. (2018); Xu and Balents (2018); Yuan and Fu (2018); Po et al. (2018); Roy and Juricic (2018); Guo et al. (2018); Padhi et al. (2018); Irkhin and Skryabin (2018); Dodaro et al. (2018); Zhang (2018). The fact that the phase diagram, when squinting, looks similar to that of cuprates or pnictides made this discovery even more exhilarating. Namely, at half-filling of the double degenerate flat band just below the Dirac cones a Mott insulator was found, and superconductivity appeared upon doping of this Mott state.

The nature of a Mott insulating phase, especially when dealing with degenerate bands, can only be fully understood by studying the real-space structure of the orbitals that are involved. However, in the case of TBG at the magic angle, the unit cell consists of 11,164 atoms. This is not very suitable for a simple model, so we need to reduce the number of orbitals involved. To this purpose, we constructed an effective model consisting of eight orbitals taken from the full tight binding problem, see Fig. 1. At the -point, the low energy states are localized at the center of the AA stacking regions. At the -point, however, the wavefunctions form ‘rings’ around the AA centers. We show that the interaction energy is minimized in the insulating phase by a charge-transfer from the ‘center’ to the ‘ring’ orbital.

In the following we first discuss the symmetries and structure of the twisted bilayer graphene lattice. We then construct the real-space wavefunctions and a corresponding effective tight binding model. Using the hopping model and general energy considerations, we show that the Mott phase is a spin singlet paramagnet. Doping this paramagnet can lead to superconductivity, and the different pairing symmetries for different doping are discussed. We end this Letter with a brief discussion of experimental signatures of the suggested charge-transfer.

The model - Our aim is to derive a simple model that captures the essential properties of twisted bilayer graphene. A single layer of graphene has a honeycomb structure, with lattice unit vectors where nmReich et al. (2002); Castro Neto et al. (2009). The unit cell contains two inequivalent sites, labeled A and B. The corresponding band structure is known for having Dirac cones at and . The symmetry group is , which means the unit cell contains one 6-fold rotation center, two 3-fold rotation centers and six distinct reflection axes.

The band structure of a bilayer depends on how the two layers are stacked. A literal stacking of two layers on top of each other, where the A sites in one layer are on top of the A sites of the other, is known as AA stacking. A more natural stacking occurs when the A sites of one layer are above the B sites of the other layer, this is known as AB stacking.

‘Twisted’ bilayer graphene can be made by starting with either stacking and then rotate one of the layers around its origin. To achieve a commensurate rotation, we can choose two integers and define the rotation angle as the angle between and Suárez Morell et al. (2010). The new unit cell has unit vectors and and contains atoms. This new unit cell, shown in Fig. 2, possesses a region where the layers are effectively AA stacked, as well as regions with AB stacking and BA stacking - the latter just being the same as AB stacking but with the two layers exchanged.

Interestingly, the space group of the twisted bilayer is different from the single layer one. It does contain 6-fold rotations around the AA centers, and 3-fold rotations around the AB and BA centers. However, there are no mirror symmetries, as can be inferred from Fig. 2. Reflection along one of the symmetry axes interchanges the two layers. Therefore, the space group is the double cover of , rather than .Dresselhaus et al. (2007); Bradlyn et al. (2017); Mele (2011)

The magic angle of Refs. [Cao et al., 2018a, b] was , which can be generated by choosing and . In this case, the new unit cell contains 11,164 atoms. Needless to say, an exact description of the band structure from first principles is a very difficult task. However, we can write down a tight-binding model including all the 11,164 bands using parameters from literature. We choose the in-plane nearest neighbor hopping eV and interlayer hopping described by , where is the total distance between two atoms including the interlayer distance nm, and is chosen such that eV for the AA stacked atomsSuárez Morell et al. (2010); Fang and Kaxiras (2016); Jung and MacDonald (2014); Malard et al. (2007). The resulting band structure for the magic angle system is shown in Fig. 3.

As predicted by continuum models,Bistritzer and MacDonald (2011); Lopes dos Santos et al. (2012) there are double degenerate Dirac cones at the and points, but with drastically reduced velocity. The bandwidth of the bands around the Dirac points is also severely reduced to about meV. Any effective model should include these degenerate Dirac cones and a corresponding narrow band.

To understand the Mott physics, however, a precise knowledge of the real space structure of the bands is necessary. We display the spatial structure of the low energy wavefunctions at and in Fig. 1. Consistent with tunnelling experiments,Wong et al. (2015) the electrons close to charge neutrality are strongly localized at the AA regions. However, at the point a clear transfer of charge is seen to a ring around the AA center. Any effective model of the whole band should capture this charge-transfer. We will therefore include all of the 8 orbitals shown in Fig. 1.

Let us focus first on the -orbitals or ‘center’ orbitals, whose weight at the AA centers form a triangular lattice. Yet the symmetry of the corresponding band structure including Dirac cones suggests a honeycomb-like structure! However, in TBG the existence of Dirac cones is not tied to a honeycomb structure, but to the nature of the little co-group at the point. The and are exchanged under six- and two-fold rotations with respect to the axis perpendicular to the TBG plane. As a direct consequence, the little co-group at the -point is , a subgroup of , the point group of TBG. Orbitals at this point could realize irrep and irrep of . Consider the equivalent representation at which involves the transformation of under the symmetry operations of . Rotation of the radius vector by the angle anticlockwise is equivalent to rotation of the vector in the opposite direction, that is to substitution of the three equivalent corners of the small Brillouin zone. When this reducible representation is decomposed into irreps of the group , we find that it is exactly contained in the two-dimensional irrep of . Therefore orbitals at the point should be four-fold degenerate, taking into account the double cover. In practice, the debate between ‘team triangular’Dodaro et al. (2018); Xu and Balents (2018); Efimkin and MacDonald (2018); Guo et al. (2018) and ‘team honeycomb’Zhang (2018); Yuan and Fu (2018) was elegantly circumvented by Po et al.Po et al. (2018), who showed that it is possible to have a honeycomb band-structure while having spectral weight centered at the AA regions.

The -orbitals or ‘ring’ orbitals, on the other hand, have a simple interpretation in terms of a gapped triangular-lattice orbitals. In the absence of hybridization between ‘ring’ and ‘center’ orbitals, we would have four gapped triangular-lattice orbitals and four-fold degenerate Dirac cones. Hybridization between and -orbitals causes the lowest energy band to be a mix of both orbitals. To construct an effective model we split the orbitals in two degenerate ‘valleys’. Each ‘valley’ consists of two ring and two center-orbitals, and there are four parameters that determine the band structure: hopping and , hybridization between different type orbitals , and the gap of the -orbitals . The resulting band-structure is shown in Fig. 4, and this model will serve as the starting point for our analysis of the insulating phase.

Interactions and the Mott phase - The emergence of insulating behavior when non-interacting theories predict conducting can be due to either Wigner or Mott localization. A Wigner crystalPadhi et al. (2018) as the origin of the insulating behavior can be dismissed due to commensurability. Wigner crystals can exist in lattice systems at sufficiently low densities, but in these cases there is a devil’s staircase of many Wigner crystals at different wavelengths,Rademaker et al. (2013) which is not the case in TBG. Because the insulating phase only appears at half-filling of the flat bands, the insulating phase is of a Mott character.

A Mott insulator can be realized by adding to the tight-binding model an onsite Hubbard repulsion for each orbital, and it is known that single-layer graphene has a relatively strong onsite interaction eVSchuler et al. (2013); Wehling et al. (2011). However, for orbitals that span thousands of different atoms, the long-range nature of the Coulomb interaction plays a central role. The Coulomb energy is nonzero whenever there are macroscopic charge inhomogeneities ,

(1) |

where measures the deviation from the average electron density at position .

At charge neutrality, the electron charge density is evenly distributed over all carbon atoms. Charge fluctuations at the scale of the unit cell are weak: as the central limit theorem predicts they scale as , where is the area of the unit cell. However, doping away from charge neutrality causes the electronic charge density to cluster around the AA centers. Based on the full band-structure of Fig. 3, adding two electrons or holes per unit cell causes a double occupancy of the ‘center’ -orbitals. Only when the flat bands are fully filled, meaning four electrons or holes, the charge will be smoothed out by also occupying the ring-like -orbitals. Consequently, the interaction energy will be largest in the middle of the flat bands.

We thus propose that the Mott insulating state is caused by a charge-transfer from the ‘center’ () to the ‘ring’ () orbitals, to minimize the macroscopic charge inhomogeneity. Notice that this is very similar to the physics of cuprates, where a charge-transfer from the copper to the oxygen orbitals causes insulating behaviorZaanen et al. (1985).

Having established that there is one localized charge in the -orbital and one in the -orbital, there are 3 degrees of freedom remaining: spin, -valley and -valley. In multi-orbital atomic Mott insulators the Hund’s coupling favors a large total spin in the unit cellImada et al. (1998); Georges et al. (2013). However, in TBG there is a nonzero overlap between ‘ring’ and ‘center’ orbitals, which causes in 2nd order perturbation theory an effective antiferromagnetic coupling. Contrary to Hund’s expectations, the two spins in each unit cell will therefore form a singlet. We predict that the Mott phase in TBG is a non-entangled featureless spin-singlet paramagnet, consistent with the experimentally measured magnetic responseCao et al. (2018a). Also note that this is the natural ground state for a model with an even number of spins per unit cellHastings (2004).

Similarly, second order perturbation theory suggests a ‘valley antiferromagnetic’ coupling between orbitals at neighboring sites. For the charge in the -orbitals, whose hopping pattern follows honeycomb symmetry, this implies a two-sublattice alternating valley order. For the -orbitals, the triangular symmetry implies a three-sublattice valley ordering.

Superconductivity - Consider the insulating state that has two electrons per unit cell removed relative to charge neutrality. When doping away from this Mott state, the dynamics of the dopants is described by a model. A key feature of such models is the effective nearest-neighbor attraction between dopants. This attraction can lead to superconductivity of the dopants. The symmetry of the pairing state depends crucially on the nature of the dopants themselves.

A prominent consequence of the proposed charge-transfer is the difference between electron and hole-doping relative to the insulating phase. Adding electrons to the insulating state - that is, moving closer to charge neutrality - will add dopants on the -orbitals. The effective model has nearest-neighbor attraction on a triangular lattice. This has been studied before and the most likely superconducting state would be spin-singlet -waveChen et al. (2013); Guo et al. (2018).

On the other hand, hole-doping adds carriers to the -orbitals, which realize an effective honeycomb lattice. Nearest-neighbor attraction on a honeycomb lattice leads to exotic spin singlet -wave superconductivity, as was proposed for single layer graphene away from charge neutralityUchoa and Castro Neto (2007). A symmetry difference between the electron and hole-doped superconducting phases relative to the Mott state would be a clear proof of the charge-transfer occurring in TBG.

Outlook - We showed that the observed insulating state in TBG can be described in terms of charge-transfer from ‘center’ to ‘ring’ orbitals around the region of AA stacking. An experimental signature of this transfer can be found in tunnelling experiments: the electron density at the AA region center should be less than expected based on non-interacting theories. Note that the lattice relaxation of the Moiré patterns due to electron-phonon coupling can influence the expected charge densitySlotman et al. (2015); van Wijk et al. (2015); Nam and Koshino (2017). Therefore a full first-principles computation of the interactions in magic angle twisted bilayers is necessary to quantify the suggested charge-transfer.

A direct consequence of the charge-transfer is that the Mott phase is a featureless spin singlet paramagnet. The lowest energy spin excitations will be propagating triplets, who can be observed using thin film resonant inelastic X-ray scatteringAment et al. (2011). Also the symmetry difference between electron- and hole-doped superconductors is a result of the charge-transfer, and is observable in experiments.

Finally, we want to emphasize that we based our predictions on a simple analysis of the real space wavefunctions of an 11,164-bands model. The effective model is a hybrid mixture of triangular and honeycomb symmetries, and it is not a trivial task to construct a low energy effective model out of those ingredients. However, we think that developing such a model and studying it using both analytical and numerical methods might provide key insights towards the understanding of twisted bilayer graphene.

###### Acknowledgements.

Acknowledgments - We are thankful to T. Hsieh, J. Venderbos and M. I. Katsnelson for discussions. P. M. acknowledge Fondecyt Grant No. 1160239. This research was supported in part by Perimeter Institute for Theoretical Physics. Research at Perimeter Institute is supported by the Government of Canada through the Department of Innovation, Science, and Economic Development, and by the Province of Ontario through the Ministry of Research and Innovation.## References

- Cao et al. (2018a) Y. Cao, V. Fatemi, A. Demir, S. Fang, S. L. Tomarken, J. Y. Luo, J. D. Sanchez-Yamagishi, K. Watanabe, T. Taniguchi, E. Kaxiras, et al., Nature 556, 80 (2018a).
- Cao et al. (2018b) Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxiras, and P. Jarillo-Herrero, Nature 556, 43 (2018b).
- Efimkin and MacDonald (2018) D. K. Efimkin and A. H. MacDonald, arXiv p. 1803.06404 (2018).
- Pal et al. (2018) H. K. Pal, S. Spitz, and M. Kindermann, arXiv p. 1803.07060 (2018).
- Xu and Balents (2018) C. Xu and L. Balents, arXiv p. 1803.08057 (2018).
- Yuan and Fu (2018) N. F. Q. Yuan and L. Fu, arXiv p. 1803.09699 (2018).
- Po et al. (2018) H. C. Po, L. Zou, A. Vishwanath, and T. Senthil, arXiv p. 1803.09742 (2018).
- Roy and Juricic (2018) B. Roy and V. Juricic, arXiv p. 1803.11190 (2018).
- Guo et al. (2018) H. Guo, X. Zhu, S. Feng, and R. T. Scalettar, arXiv p. 1804.00159 (2018).
- Padhi et al. (2018) B. Padhi, C. Setty, and P. W. Phillips, arXiv p. 1804.01101 (2018).
- Irkhin and Skryabin (2018) V. Y. Irkhin and Y. N. Skryabin, arXiv p. 1804.02236 (2018).
- Dodaro et al. (2018) J. F. Dodaro, S. A. Kivelson, Y. Schattner, X.-Q. Sun, and C. Wang, arXiv p. 1804.03162 (2018).
- Zhang (2018) L. Zhang, arXiv p. 1804.09047 (2018).
- Reich et al. (2002) S. Reich, J. Maultzsch, C. Thomsen, and P. Ordejón, Phys. Rev. B 66, 2794 (2002).
- Castro Neto et al. (2009) A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, and A. K. Geim, Rev. Mod. Phys. 81, 109 (2009).
- Suárez Morell et al. (2010) E. Suárez Morell, J. D. Correa, P. Vargas, M. Pacheco, and Z. Barticevic, Phys. Rev. B 82, 121407 (2010).
- Dresselhaus et al. (2007) M. S. Dresselhaus, G. Dresselhaus, and A. Jorio, Group theory: application to the physics of condensed matter (Springer Berlin Heidelberg, 2007).
- Bradlyn et al. (2017) B. Bradlyn, L. Elcoro, J. Cano, M. G. Vergniory, Z. Wang, C. Felser, M. I. Aroyo, and B. A. Bernevig, Nature 547, 298 (2017).
- Mele (2011) E. J. Mele, Phys. Rev. B 84, 9242 (2011).
- Fang and Kaxiras (2016) S. Fang and E. Kaxiras, Phys. Rev. B 93, 235153 (2016).
- Jung and MacDonald (2014) J. Jung and A. H. MacDonald, Phys. Rev. B 89, 035405 (2014).
- Malard et al. (2007) L. M. Malard, J. Nilsson, D. C. Elias, J. C. Brant, F. Plentz, E. S. Alves, A. H. Castro Neto, and M. A. Pimenta, Phys. Rev. B 76, 201401 (2007).
- Bistritzer and MacDonald (2011) R. Bistritzer and A. H. MacDonald, Proc. Natl. Acad. Sci. 108, 12233 (2011).
- Lopes dos Santos et al. (2012) J. M. B. Lopes dos Santos, N. Peres, and A. H. Castro Neto, Phys. Rev. B 86, 155449 (2012).
- Wong et al. (2015) D. Wong, Y. Wang, J. Jung, S. Pezzini, A. M. DaSilva, H.-Z. Tsai, H. S. Jung, R. Khajeh, Y. Kim, J. Lee, et al., Phys. Rev. B 92, 155409 (2015).
- Rademaker et al. (2013) L. Rademaker, Y. Pramudya, J. Zaanen, and V. Dobrosavljević, Phys. Rev. E 88, 032121 (2013).
- Schuler et al. (2013) M. Schuler, M. Rosner, T. O. Wehling, A. I. Lichtenstein, and M. I. Katsnelson, Phys. Rev. Lett. 111, 036601 (2013).
- Wehling et al. (2011) T. O. Wehling, E. Sasioglu, C. Friedrich, A. I. Lichtenstein, M. I. Katsnelson, and S. Blugel, Phys. Rev. Lett. 106, 236805 (2011).
- Zaanen et al. (1985) J. Zaanen, G. A. Sawatzky, and J. W. Allen, Phys. Rev. Lett. 55, 418 (1985).
- Imada et al. (1998) M. Imada, A. Fujimori, and Y. Tokura, Rev. Mod. Phys. 70, 1039 (1998).
- Georges et al. (2013) A. Georges, L. d. Medici, and J. Mravlje, Annu. Rev. Condens. Matter Phys. 4, 137 (2013).
- Hastings (2004) M. B. Hastings, Phys. Rev. B 69, 104431 (2004).
- Chen et al. (2013) K. S. Chen, Z. Y. Meng, U. Yu, S. Yang, M. Jarrell, and J. Moreno, Phys. Rev. B 88, 041103 (2013).
- Uchoa and Castro Neto (2007) B. Uchoa and A. H. Castro Neto, Phys. Rev. Lett. 98, 496 (2007).
- Slotman et al. (2015) G. J. Slotman, M. M. van Wijk, P.-L. Zhao, A. Fasolino, M. I. Katsnelson, and S. Yuan, Phys. Rev. Lett. 115, 186801 (2015).
- van Wijk et al. (2015) M. M. van Wijk, A. Schuring, M. I. Katsnelson, and A. Fasolino, arXiv (2015).
- Nam and Koshino (2017) N. N. T. Nam and M. Koshino, Phys. Rev. B 96, 075311 (2017).
- Ament et al. (2011) L. J. P. Ament, M. van Veenendaal, T. P. Devereaux, J. P. Hill, and J. van den Brink, Rev. Mod. Phys. 83, 705 (2011).

## Appendix A Appendix: Hamiltonian of effective model

Our effective model has 8 orbitals centered around the AA regions of the large unit cell. The Hamiltonian can be expressed as a block diagonal matrix as a function of momentum , where each block is large.

The honeycomb structure of the -orbitals is reflected in the hopping factors

(2) |

whereas the triangular nature of the -orbitals is shown in the factors

(3) |

The coupling between and -orbitals can have two shapes,

(4) | |||||

(5) |

reflecting the two possible ways to combine honeycomb and triangular symmetries.

Each -block of the Hamiltonian now reads

(6) |

with the following parameters:

(7) | |||||

(8) | |||||

(9) | |||||

(10) |

## Appendix B Orbital nature throughout the Brillouin zone

In the main text we have discussed the orbital character of the low-energy states at and . At momenta in between, the orbital character smoothly changes from ring-like to center-like, as is shown in the figure below.

The overlap is defined as the norm of the wavefunction at momentum projected onto the subspace of the four lowest energy wavefunctions at either (‘ring’) or (‘center’). Notice that the overlap between ring and center orbitals is .