Correlated Insulators in Twisted Bilayer Graphene

Experiments on graphene bilayers, where the top layer is rotated with respect to the one below, have displayed insulating behavior when the moir\'e bands are partially filled. We calculate the charge distributions in these phases, and estimate the excitation gaps.


I. Introduction
Graphene bilayers, where the top layer is rotated with respect to the bottom, show remarkable properties 1-10 . These arise from the presence of a long wavelength moiré superlattice. For twists near certain "magic angles", the low energy bands become very flat, and interactions dominate 11 . Moreover, the bands have non-trivial topological indices. Experimentally one observes insulating phases at certain rational fillings of the bands. Electrostatically doping away from these rational fillings leads to superconducting phases, whose transition temperatures are large compared to the bandwidth 1,2 . An important part of understanding the physics of these systems is to identify the structure of the correlated insulating states. Here we conduct a variational study of the various possible charge-density, spin-density, and valley-density waves, which are the most natural candidates. We find that there are "stripe" ordered spin and valley ferromagnets at fillings of {1/8, 3/8, 5/8, 7/8}. At fillings of {1/4, 1/2, 3/4}, the ground state is a spatially homogeneous spin and valley ferromagnet.
There is significant prior work on this problem. Two separate groups, Seo, Kotov, and Uchoa 12 , and, Kang and Vafek 13 recently argued for a ferromagnetic state at 1/8 filling. Kang and Vafek's results are similar to ours, in that they find charge ordering in addition to the spin ordering. Our approach is complementary in that we explore different models for the electron-electron interactions. Kang and Vafek use a short-range interaction, while we consider two models: a long-range Coulomb interaction, and a dipolar interaction which accounts for image charges in the back-gate. All three models give similar results, pointing to the robustness of the phenomena.
Related analysis is found in the supplementary material of Choi et al. 3 , which considers the half-filled case. Lu et al. 14 found signatures of an orbital ferromagnetic state at 3/8 filling. Also notable is the slave-spin treatment of Pizarro, Calderon, and Bascones 15 , the Hartree-Fock analysis of Xie and MacDonald 16 , and the projection-based technique of Repellin, Dong, Zhang, and Senthil 17 . All of these reveal different aspects of (Color Online) Schematic of twisted bilayer graphene: Green/orange discs represent carbon atoms in each layer. In the lower left, lines connect nearest-neighbor atoms in each layer. AA (pink), AB (red), and BA (blue) regions are labeled. Two schematic trefoil-shaped Wannier states are shown, which centered on BA and AB regions respectively. the correlated insulators.
In addition to those already described 1,2 , experimental studies of the correlated insulators include measurements of the compressibility 3,4 , magnetic response 5 , tunneling spectroscopy 6-8 and further transport properties 9,10 .
Beyond finding the lowest energy charge configurations, we investigate the energy cost of adding a particle, adding a hole, or adding a particle-hole pair. We find that the interaction energy from these defects is large compared to the width of the moiré mini-bands. Consequently, one expects large interaction-driven mixing in of higher bands, which would need to be included in quantitative models. arXiv:2008.13752v1 [cond-mat.str-el] 31 Aug 2020

II. Model
Graphene forms a honeycomb lattice, with a unit cell containing two sites -denoted A and B. For small twistangles (near 1 • ), the bilayer system displays a large moiré unit cell. There are three notable regions, denoted by AA, AB, and BA, which each form triangular lattices (see Fig. 1). In the AA regions, the two lattices align. In the AB regions, the A sites of the lower layer line up with the B sites of the upper layer. The BA region is the opposite. The AB and BA sites together form a honeycomb lattice.
As described in Ref. 18, the electronic structure of the low energy bands is built from Wannier states which are centered at the AB and BA sites. Each Wannier state has a three-lobed spatial structure, with each lobe centered on an AA region. Thus the Wannier state centers are in the AB-BA regions, but the charge density is peaked in the AA regions. This structure is schematically shown in Fig. 1. In addition to a sublattice index (AB or BA), the Wannier states in this model are labeled by spin and valley indices, leading to an 8-band effective Hamiltonian.
There is a number of subtleties about the construction of the Wannier states. For example, Refs. [19][20][21] argued that there is a topological obstruction that prevents the existence of exponentially localized Wannier states obeying all the emergent symmetries. Indeed, while the Wannier states from Ref. 18 that we work with are exponentially localized, one of the protecting symmetries is broken 20,21 . These considerations only affect physics on an exponentially small energy scale, where symmetry breaking effects matter: They are not relevant to the questions we are asking, and our results would be unchanged if we worked with a larger basis set.
The effective electronic hopping matrix element between two nearest Wannier orbitals is about 0.3 meV, while the effective direct Coulomb interaction between them is of the order of 10 meV 12,18 . This clear separation of scales motivates us to consider stationary charge distributions, where we break translational symmetry, and distribute the electrons among various Wannier orbitals. We neglect kinetic energy and find the configurations which minimize the Coulomb interaction energy.
We systematically consider periodic arrangements of charges defined by the generators of the supercell v 1 = m 1 a 1 + m 2 a 2 and v 2 = n 1 a 1 + n 2 a 2 , where a 1 and a 2 are the generators of the moiré lattice. We include all unit cells containing fewer than 13 sublattice sites, with each of m 1 , m 2 , n 1 , and n 2 restricted to be smaller than or equal to 4. We also include a small number of larger unit cells. In all, we consider well over 3000 charge configurations -though many of them are equivalent, or symmetry related.
We use the interaction model derived by Koshino et al. 18 . They found that the Coulomb interaction energy from lobes in the AA regions r i and r j is well approximated by where −e/3 is the charge in each lobe, is the dielectric constant, which we take equal to 10 0 . Here L M is the distance between neighboring AB and BA regions. We use L M = 13.4 nm, corresponding to a twist angle of θ = 1.05 • . We also include an exchange interaction, which lowers the energy by J when two overlapping occupied orbitals have the same spin and valley quantum numbers. We use the values for the exchange energies from Ref. 18. These are dominated by the nearest neighbor and next-nearest neighbor terms: J = 0.376, 0.0645 e 2 /( L M ). All other exchange energies are smaller than 0.0014 e 2 /( L M ). Terms beyond the fifth-nearest neighbour are ignored. The Coulomb interaction energy from a 2D array of charges is linearly divergent with system size. The relative energies of the configurations with the same average density, however, are well defined. Some care must be taken in calculating these energies, since even after subtracting off the leading divergence, the remaining sums are conditionally convergent. The process is regularized by using 2D Ewald sums 22,23 : Here η is an arbitrary parameter which separates the interactions into long-range and short-range parts. The sum is over a Bravais lattice L, which in our case is generated by the supecell vectors v 1 , v 2 . The reciprocal lattice isL, andL\0 represents the reciprocal lattice with the origin removed. Here C is an irrelevant infinite constant, which physically represents the Coulomb energy from a uniform charge distribution. The sums in Eqs. (3) and (4) are absolutely convergent. We chose η to give a reasonably fast convergence rate -typically taking it to be 1/3 of the lattice constant.
In the experimental set-up of Ref. 2, a metallic backgate sits roughly 30 nm behind the sample. We supplement our long-range Coulomb calculations by also calculating energies where we include a set of image charges in this layer. Appendix A gives the resulting expressions for the energies. To leading order, these images uniformly shift the energies of all configurations by the same amount (corresponding to the classical capacitance). Configuration-dependent corrections are exponentially small in the ratio of the supercell periodicity to the distance from the back-gate, and we find that for the experimental geometry the results for ground state energies are unchanged when we add the images.

III. Ground States
As previously explained there are 8 states per unit cell -labeled by sublattice (AB or BA), spin, and valley. Due to the relatively strong exchange interaction, physically corresponding to large spatial overlap between neighboring Wannier states, it is always favorable to maximize the number of electrons which have the same spin and valley quantum numbers. Thus we find that at fillings ν = 1/8, 1/4, all electrons have the same spin and valley quantum numbers. At ν = 3/8, there are not enough states for all electrons to carry the same spin and valley quantum numbers. Instead, two-third of the electrons carry one set, while the other one-third carry another set. At ν = 1/2, things are similar, but an equal number of electrons carry each of the two sets of quantum numbers. Larger ν's are simply particlehole mirrors of these configurations. At the level of our model, the exact quantum numbers do not matterwhat matters is just that they follow this pattern.
We find that for ν = 1/8, the occupied modes (all of which have the same spin and valley quantum numbers) form a striped configuration, as shown in Fig. 2. For ν = 1/4, the occupations are uniform, with every Wannier state occupied by a single electron -all of which share the same spin and valley quantum numbers. For ν = 3/8, the electrons with the dominant quantum numbers form a uniform pattern, while the ones with the other Despite the symmetry breaking at ν = 1/8, the charge density is nearly uniform: Each AA site has a charge of −e. The spatial symmetry breaking, however, should be apparent in optical measurements through a birefringence. One may also find that transport measurements are anisotropic.
The energy per particle in each of these four cases is E = −1.72, −2.39, 35.5, 49.4 e 2 / ( L M ) . In physical units, this corresponds to E = −18.4, −25.6, 381, 530 meV per particle. The relatively small difference between the energy per particle for ν = 1/8 and ν = 1/4 occurs because in these spin/valley polarized states, the direct Coulomb interaction is largely cancelled by the exchange interactions. At higher fillings, we are adding electrons with different spin/valley quantum numbers, for which there is no such cancellation. These energies should all be understood as relative to a classical uniform charge distribution.
One way to quantify the stability of the striped state at ν = 1/8, is to look at the next-lowest energy state that we found (Fig. 3). This excited configuration has energy −17.7 meV per particle, which is only modestly higher than the ground state. The charge density for this pattern is again uniform, even though the pattern of occupied Wannier states is non-uniform. Note, one can create ν = 1/8 configurations whose energies lie between these by alternating these two patterns.
Given the small energy difference, lattice defects, impurities, or other irregularities could play a role in the experimentally observed pattern. Furthermore, kinetic effects could mix them.

IV. Excitations
We calculate the energy of four types of excitations: Adding a single electron to one of our ground states, removing a single electron (i.e. adding a hole), adding an electron-hole pair, or flipping the spin/valley of a single electron: E e , E h , E eh , E s . Physically, E e and E h correspond to the chemical potentials required to add or remove a charge and are experimentally determined by measuring the charge density as a function of gate voltage: the charge density is constant over a voltage range of width ∆V = (E e + E h )/e. The electron-hole excitation energy can experimentally be measured from transport: one expects the low temperature resistivity to obey an Arrhenius law with energy scale E eh . The spin excitation energy can similarly be inferred from the temperature dependence of the spin susceptibility. While we find that the presence of a back-gate plays no role in the ground state energies, it can shift the excitations energies by ∼ 1 meV. We just quote numbers for the case without the back-gate.
Experiments find the energy scales for excitations are fractions of meV, while our model yields tens of meV. The discrepancy has a number of sources. First, since the Coulomb interaction scale is large compared to the width of the moiré mini-bands, the low energy excitations undoubtedly involves electronic states which are outside of the minibands (and hence not included in our model). Second, it is likely that impurities and other defects are playing a role in the experiment. Third, collective effects can lower the excitation energy. Nonetheless, our estimates are quite valuable as a starting point for understanding the phenomena.
For ν = 1/8, the lowest energy electron excitation consists of adding a particle with the same quantum numbers as the other atoms to one of the empty sites: all possibilities are equivalent. This has energy E e = −1.33 × e 2 L M = −14.3 meV. Similarly, all ways of removing an electron are equivalent, with E h = 36.9 meV. As is required for stability, −E h < E e , which is a discrete analogy to requiring that the compressibility is positive, ∂ 2 E/∂N 2 > 0. The fact that E e < 0 should not cause alarm -this is a common physical occurrence. The physical picture is that our energies are measured relative to having a uniform charge distribution, and we effectively have contributions from a uniform positive background.
The lowest energy particle-hole excitation at ν = 1/8 occurs when the particle is moved to a neighboring unoccupied site. This has the energy E eh = 10.2 meV, which is smaller than E e + E h because of the electronhole interactions. It is also larger than the energy-perparticle difference between the configurations in Fig. 2 and Fig. 3. One can generate Fig. 3 from Fig. 2 by moving one-fourth of the particles. The interactions between these particle-hole pairs reduce the cost of this rearrangement. A spin or valley flip excitation costs an energy of E s = 10.1 meV, which is solely due to the exchange interaction.
For ν = 1/4, any extra electron must have different quantum numbers than the others. All sites are equivalent, and we find E e = −13.8 meV. A naive expectation might be that this excitation energy would be twice what we found at ν = 1/8, since the charge is twice as large. In fact adding an electron at ν = 1/4 costs much more than twice the ν = 1/8 value. Given that the charge density is very same in both cases, this difference can be attributed to the exchange interaction. The energy to remove an electron is E h = 51.2 meV, which is somewhat less than twice the energy of the equivalent excitation at ν = 1/8. The difference can be attributed to the detailed structure of the charge distribution in the Wannier states. The particle-hole excitation energy is E eh = 21.1 meV. A spin or valley flip excitation costs an energy of E s = 17.5 meV.
At ν = 3/8, the ground state contains electrons with two different quantum numbers. Let us denote these by α and β. The α electrons are spread uniformly, while the β electrons form the zigzag pattern in Fig. 2. The lowest energy electronic excitation corresponds to adding another electron in the β state, with an energy of E e = −28.1 meV. Clearly, this is just the sum of the values of E e for the ν = 1/8 and ν = 1/4 configurations. Note that if we add an electron with a different spin/valley, the value would have been −20.7 meV, which is energetically costlier. The lowest energy hole excitation comes from removing an electron from the β state, which yields E h = 33.2 meV. The particle-hole excitation energy is E eh = 10.2 meV. A spin or valley flip excitation costs an energy of E s = 10.1 meV.
At ν = 1/2, any added electron must have different quantum numbers than the rest of electrons. The required energy is E e = −27.6 meV. As expected, this comes out to be twice the value of E e for ν = 1/4. The energy required to remove an electron is E h = 65.0 meV. The particle-hole excitation energy is E eh = 21.1 meV. A spin or valley flip excitation costs an energy of E s = 17.5 meV. We note that the value of E s is the same for ν = 1/8 and ν = 3/8 (equal to 10.1 meV), and for ν = 1/4 and ν = 1/2 (equal to 17.5 meV).
Much of the structure in the excitation energies is related to the exchange interaction: The direct Coulomb energy (in units of e 2 L M ) for adding a particle to the grounds states corresponding to 1/8, 1/4, 3/8, 1/2 fillings are: 0.642, 1.284, 1.926, 2.568, which are in the ratio 1 : 2 : 3 : 4. This scaling arises because the charge distribution is uniform in all cases.

V. Summary and Outlook
To summarize, we have calculated the charge distributions of the correlated insulators at all integer fillings of magic angle twisted bilayer graphene. We have also estimated all possible excitation gaps over the insulating ground states. In all these computations, we have neglected the kinetic terms, as the effective tightbinding models show that the strength of the hopping terms is much much less than the Coulomb interactions, and will only have a perturbative effect. One essential caveat is that because the Coulomb scale is so large, it likely leads to mixing between the flat bands and the spectator bands.
It would be rewarding to apply similar techniques to characterize the correlated insulating phases of transition metal dichalcogenide (TMD) homobilayers and heterobilayers 24-28 , as they have much simpler moiré band structures. Moreover, in certain TMDs like WSe 2 , the flat bands and the resulting correlated states are found to exist over a continuum of angles 29 rather than a narrow range around some "magic angles".

VI. Acknowledgements
This work was supported by the NSF Grant PHY-1806357.

A. Interactions with Image Charges
We model screening by a back-gate sitting a distance d/2 from the sample by placing equal and opposite image charges at a distance d. We treat the charge distribution of each Wannier state as three delta functions, each of charge q = e/3. The Coulomb interaction between a single charge in one plane, and a Bravais lattice of charges in another is where λ = q 2 /(4π ), w is a 2D lattice vector, and r is the 2D location of the charge in the device. By Fourier transforming the potential, one can write: where the first term contains the same infinite constant as Eq.
(2). The second term is the potential from a uniformly charged plane, with Ω being the area of the unit cell. The third term falls off exponentially with the separation. ThereL is the 2D reciprocal lattice. Due to this exponential cutoff, the image charges play a very minor role in the energetics of different periodic charge configurations. We follow the procedure in Lee and Cai 23 to write