The nature of correlations in the insulating states of twisted bilayer graphene

The recently observed superconductivity in twisted bilayer graphene emerges from insulating states believed to arise from electronic correlations. While there have been many proposals to explain the insulating behaviour, the commensurability at which these states appear suggests that they are Mott insulators. Here we focus on the insulating states with ±2 electrons or holes with respect to the charge neutrality point. We show that the theoretical expectations for the Mott insulating states are not compatible with the experimentally observed dependence on temperature and magnetic field if, as frequently assumed, only the correlations between electrons on the same site are included. We argue that the inclusion of non-local (inter-site) correlations in the treatment of the Hubbard model can bring the predictions for the magnetic and temperature dependencies of the Mott transition to an agreement with experiments and have consequences for the critical interactions, the size of the gap, and possible pseudogap physics. The importance of the inter-site correlations to explain the experimental observations indicates that the observed insulating gap is not the one between the Hubbard bands and that antiferromagnetic-like correlations play a key role in the Mott transition.


Introduction
Unexpected insulating states have been recently observed upon doping a graphene bilayer with a small twist angle ∼1.08°−1.16° [1]. The stacking misorientation creates a moiré pattern with a superlattice modulation corresponding to thousands of atoms per unit cell. The insulating states, believed to be due to electronic correlations [1], were originally observed when the charge per moiré cell is ±2 with respect to the charge neutrality point (CNP). More recently, insulating states when the twisted bilayer graphene (TBG) is doped with 1 or 3 carriers have also been observed [2,3].
The experimental characterization of the insulating states with ±2 electrons or holes, in which we focus from now on, reveals that the system becomes metallic with increasing temperature and in a magnetic field [1]. In addition, there is no sign of symmetry breaking due to magnetic order. Superconductivity emerges from the insulating statel [2][3][4]. Understanding the nature of the insulating states is key to uncover the origin of the superconductivity.
The band structure of TBG at the charge neutrality point features two doubly-degenerate Dirac cones [5]. The nearby bands at both positive and negative energies become very narrow at the experimental twist angle [6]. The electronic charge in these narrow bands is believed to be concentrated at the center of each moiré cell, in the so-called AA regions which form a triangular superlattice [7].
Several authors have proposed different origins for the insulating states, many of them relying on specific properties of the density of states at the Fermi level [8][9][10][11][12][13][14][15]. However, an important clue comes from the fact that the fillings at which these states are experimentally found correspond to an integer number of electrons/ holes per moiré unit cell, as for Mott states in Hubbard models. In fact, the exact commensurability has been now observed for a larger range of twist angles [1][2][3], for which considerably different density of states are Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. expected. The proximity of a metallic gate, at a distance of only 1−2 moiré lattice constants from the TBG [1], is expected to screen the long-range character of the Coulomb interaction. For all these, an effective Hubbard model for the moiré superlattice seems a plausible starting point to describe the insulating and superconducting states in a TBG. Due to the valley degree of freedom the effective Hubbard model should include at least two orbitals [1].
In contrast to single orbital systems, for which the insulating Mott state appears only at half-filling (one electron per site in the single orbital case), in multi-orbital Hubbard models Mott states are also found for other fillings with an integer number of electrons per site [16]. At large interactions the energy of the system is reduced if an integer number of electrons is localized at each atomic site and the system goes through a metal-insulator Mott transition at a critical on-site interaction U c which depends on the filling and on other interaction parameters such as Hund's coupling J H . Though localized electrons have a tendency to order magnetically, the Mott transition does not require any symmetry breaking. For interactions smaller than the critical interaction for the Mott transition, a metallic correlated state is found with the spectral weight partially redistributed to many-body Hubbard bands.
Local correlations, namely correlations between electrons on the same site, induce a Mott transition at integer fillings via the opening of a charge gap. Approaches usually used to address the Mott transition, such as single-site Dynamical Mean Field Theory (DMFT) [17] or slave particle techniques [18][19][20][21], only include these on-site correlations. However, even if in the Hubbard model interactions are restricted to electrons on the same site, inter-site (non-local) correlations are established. The interplay between these magnetic and orbital nonlocal correlations and the charge degree of freedom affects significantly the Mott transition. Here we show that experimental observations on the insulating states in TBG are not consistent with expectations including only local correlations. We support our claim by studying their effect on the two-orbital Hubbard model on the honeycomb lattice for the moiré superstructure with a U(1) single-site slave-spin mean-field technique [21]. We argue that the inclusion of non-local correlations can heal this disagreement, indicating that antiferromagneticlike correlations play a key role in the localization of the carriers.

Model
The hamiltonian for the moiré superlattice comprises the kinetic energy H t and the on-site interaction term H U . The tight-binding models for the moiré superlattice are not defined in the basis of the atomic p z orbitals of the carbon atoms, but on effective orbitals from Wannier projections for the low energy bands. These effective orbitals are not defined in a particular graphene layer but in the TBG. However, the most suitable model is still under discussion [8-10, 12-15, 22-26]. Although the AA regions form a triangular lattice, tight binding models with such symmetry spanning only the flat bands do not preserve the Dirac cones [10]. Tight-binding models for different lattices and number of bands have been proposed.
Here we focus on the simplest proposed model which reproduces the Dirac cones: a two-orbital tightbinding hamiltonian on the honeycomb lattice [8][9][10][11]22] with orbitals presumably centered on the AB and BA regions featuring lobes directed toward the AA stacking centers.
The energy bands of the two-orbital honeycomb model [8][9][10][11]22] mimic the flat bands around the charge neutrality point with bandwidth W∼10 meV [6]. There are two bands with energy E>0 and two bands with energy E<0 originating in the valley degree of freedom of each of the graphene layers. For simplicity, we only include the intraorbital hopping t to first nearest neighbors (t∼2 meV). Although this assumption makes the model particle-hole symmetric with respect to half-filling, see the inset of figure 1 for the density of states, it does not affect our conclusions. Below, all energies are in units of t.
The filling per site in the honeycomb lattice is defined as x=n/2N with n (N) the number of electrons (orbitals) per site. Note that we define the filling from the bottom of the bands, not from the charge neutrality point as in the original experimental papers [1,4]. In our notation, the filling at the Dirac points (the charge neutrality point of the TBG) corresponds to half-filling x=1/2 (on average 2 electrons per site, one per orbital, in the honeycomb lattice, 4 per moiré cell). Experimentally, at this filling the system has the semi-metallic character expected from the Dirac cone dispersion and the density of states vanishes, see the inset of figure 1. At x=1/4 (n=1) and x=3/4 (n=3), the density corresponding to the m2 electrons with respect to the CNP, the Fermi surface is centered around Γ and the density of states does not show any special features.
For multi-orbital systems, the interaction term H U includes the intra-and inter-orbital interactions U and U′, the Hund's coupling J H , and pair-hopping ¢ J terms [8,27]: with γ and β labelling the orbitals, j the sites and s s , the spin. Here ¢ = J J H and ¢ = - 27]. U has been estimated to be ∼20 meV [1]. The value of the Hund's coupling is not known. Yuan and Fu [8] assume J H =0 and = ¢ U U . On the other hand, Dodaro et al [14] argue that once the effect of the phonons of the TBG with frequencies ω∼200 meV is included the effective J H becomes negative, as for alkali-doped fullerides [28]. We do not make any a priori assumption on the value or the sign of J H . We assume the equality ¢ = -U U J 2 H remains valid for J H <0 [28].
We analyze this model with a single site U(1) slave-spin mean-field [21] technique. In the slave-spin approach the electron operator is written in terms of slave-spin and auxiliary fermion operators. This procedure enlarges the Hilbert space and the unphysical states are eliminated requiring the system to satisfy a certain constraint. In the single-site approximation the operators in different lattice sites are assumed to be uncorrelated, namely only local correlations are included. The pair-hopping and the spin-flip terms are dropped in the calculation. We refer the reader to [20,21], where the technique is explained with great detail.
The slave-spin approach has been extensively used to study electron correlations and Mott transitions in multi-orbital systems, such as iron superconductors, Hund metals and non-degenerate systems with orbital selective Mott physics [21,[29][30][31][32][33][34]. In particular, in multiorbital systems, slave-spin approaches [20,21,[29][30][31][32][33][34], in Z 2 [20] or U(1) [21] versions, have proven to be very useful to calculate the critical value of the interaction U c x . The predictions from Dynamical Mean Field Theory (DMFT) and those from slave-spin are compared in [33] for the quasiparticle weight and the critical interaction for several electronic fillings and values of the Hund's coupling. Very good agreement between both techniques is found [33]. 1 The slave-spin technique does not give information on the spectral weight redistribution which takes place with increasing interactions: Spectral weight is transferred from the quasiparticle peak to the incoherent/many body Hubbard bands, producing a three peak spectral function. We obtain this information from previous calculations performed with DMFT [17] in the literature and cited accordingly in the text.

Doping dependent critical interaction
In a two-orbital model the insulating states can be found at fillings x=1/4, 1/2 and 3/4. The insulating states with ±2 electrons with respect to the CNP are found at x=1/4 and x=3/4, with 1 and 3 electrons per site . The dotted line corresponds to U 1 2 1orb , the critical interaction for the Mott transition in a single-orbital model in the honeycomb lattice. The inset shows the non-interacting density of states as a function of electronic filling. The lower horizontal axis shows the number of electrons n per lattice site in the two orbitals and the upper axis shows the electronic filling with respect to its maximum value. Half-filling x=1/2 corresponds to the charge neutrality point. Dashed lines mark the fillings at which the Mott insulating states, corresponding to ±2 electrons from the CNP, are found. Note that the van Hove singularities are found at different fillings.
(2 and 6 per moiré unit cell), respectively, with respect to the bottom of the band. Due to the particle-hole symmetry of the model, below we restrict the discussion to  x 1 2. Agreement with the experiment requires that < U U c c 1 4 1 2 , and that the interaction U is intermediate between these two values, such that the system is an insulator at quarter-filling and a metal at half-filling. U c x depends on the charging energy cost D = . On spite of having the same charging energy cost, as seen in figure 1 and previously discussed [8,16,19,35] for small values of J H , orbital fluctuations stabilise the metallic state up to larger values of U at half-filling than at quarter-filling, such that > U U c c 1 2 1 4 , consistent with experiments [8]. From now on, we focus on the small range of Hund's coupling 1 2 is satisfied. We approximate such small values by J H =0. We note that due to the larger degeneracy of the multi-orbital insulating states, these orbital fluctuations make both U c 1 4 and U c 1 2 at small J H larger than U 1 2 1orb , the critical interaction for the Mott transition in the single-orbital model in the honeycomb lattice.

Magnetic field dependence
Experimentally, the insulating state at x=1/4 becomes metallic in a magnetic field~-H 3 6T (m~-H 0.18 0.36 B meV) [1]. This behaviour does not depend on whether the field is applied perpendicular or parallel to the TBG sample, suggesting that the suppression of the insulating state is due to the Zeeman effect. This is consistent with the fact that for localized electrons the orbital effect is not expected to play a role.
Agreement with experiments requires that the critical interaction for the Mott transition U c 1 4 increases in a magnetic field. However, contrary to the experimental observations, we find that in the local approximation a magnetic field decreases U c 1 4 promoting the insulating tendencies, see figure 2. The effect of the magnetic field H is introduced via the Zeeman term å - does not depend on H. When a magnetic field is applied, the Zeeman term favors spin polarization. When the bands are completely spin polarized (for instance, the spin up band is emptied) the behaviour of the spin down band becomes equivalent to that of a single orbital model at half-filling. As discussed above, the critical interaction for a single orbital system U 1 2 1orb is smaller than U c 1 4 at H=0. The magnetic field lifts the spin degeneracy suppressing the orbital fluctuations. As a consequence, U c 1 4 decreases with the magnetic field towards U 1 2 1orb . This is shown in figure 2. More details are given in the appendix.

Temperature dependence and gap size
The insulator to metal transition observed experimentally with increasing temperature T is also at odds with the behaviour expected from a charge gap emerging from local correlations. In single-orbital Hubbard models, the T dependence of the Mott transition in the paramagnetic phase has been extensively studied using single-site DMFT [17,36]. The results are summarized in figure 3 (see the blue line on the right) where U c,local refers to U 1 2 1orb . If at low T the system is insulating (for > U U c,local ) the resistivity decreases with increasing T, but the metallic behaviour is not restored [36]. Below U c,local the system is metallic and at a critical T there is a first order transition, blue solid line, to an insulating state. This happens because in the local approximation the entropy in the Mott state is larger than in the metallic one [17,37]. This phenomenology has been experimentally observed in oxides [38] and remains unchanged in the two-orbital case [16]. Hence, if only local correlations are included, the insulating behaviour is not suppressed with T, contrary to the experimental observations. Previous works have pointed out that the observed small gap ∼0.3 meV is also in contradiction with the expected gap in a Mott insulator [11,14]. In this statement it is implicitly assumed that the insulating gap equals the gap between the Hubbard bands~-~( ) U W U W , , with W the bandwidth, as obtained by the single-site DMFT treatment of the Hubbard model [17].
The lack of agreement between the experimental observations (dependence on magnetic field and T, and the gap size) and theoretical predictions including only local correlations is restricted neither to the two-orbital character of the Hubbard model nor to particular lattices 2 or tight-binding models. These disagreements indicate that the experimentally observed insulating state cannot be described as a Mott state considering only local correlations. Below we argue that the inclusion of inter-site correlations in the Mott transitions could be the clue to reconcile theory and experiment.

Effect of inter-site correlations
Even within the Hubbard model (with only on-site interactions) inter-site correlations are generated. The best known example is the single orbital case at half-filling for which antiferromagnetic (Heisenberg-like) correlations between the local electrons in neighboring sites emerge. In multi-orbital systems the inter-site correlations involve both magnetic and orbital degrees of freedom, and depend on J H [40,41] and the filling. For a two-orbital system with x=1/4 positive (negative) Hund's coupling promotes ferromagnetic and antiferroorbital (AFM and ferroorbital) correlations [42,43] while for J H =0 correlations are AFM and antiferroorbital [44].
Previous calculations in several lattices, mostly restricted to single-orbital models, have found that including the effect of short-range inter-site correlations in the Mott transition decreases the critical interactions and can considerably change the expected dependences [39,45]. Their effect in non-ordered states have been studied with cluster approaches [45] in the context of organic superconductors and cuprates where they play an important role. Unlike local correlations, these effects do depend on the lattice and become more important when the effective dimension is reduced, as in the honeycomb lattice, with small coordination number z=3.
Cluster treatments of multi-orbital models are computationally challenging and very scarce [40,41,44,46,47]. To our knowledge there are no studies which include these non-local correlations in the twoorbital honeycomb lattice. However, the expected qualitative behaviour in the non-ordered Mott state of the present model can be inferred from previous results in the two-orbital square lattice and in single-orbital models with different lattice geometries.
As shown in figure 3, non-local correlations shift the Mott transition to a smaller critical interaction -U c,non loc [39,44,45,[48][49][50] even in the absence of magnetic order. The correlated insulator decreases its energy Figure 3. Sketch of the temperature-interaction phase diagram for the Mott transition in the paramagnetic single-orbital case, previously obtained in DMFT (only on-site (local) correlations are included, blue lines at the right) and CDMFT (inter-site (non-local) correlations are included, green lines at the left) [17,39]. Above a given temperature the phase transition turns into a crossover (dashed lines). With only on-site correlations, the insulating behaviour is enhanced when the temperature is increased. For U below the blue line, the metallic state turns into an insulator with increasing temperature. Above the blue line, the system is always an insulator and its resistivity decreases with increasing temperature. When short-range inter-site correlations are included, the critical interaction decreases (green line) and the low T slope reverses due to the short-range singlet formation. Non-monotonic behaviour is observed and at sufficiently high-temperature both local and non-local crossover lines converge into a single one. In the region of parameters between the local and non-local lines, the system is metallic if only local correlations are considered, and insulating with inter-site correlations. In the later case, the system is insulating below the green line and metallic above it. The red arrows indicate the effect of a magnetic field on the critical interactions.
with respect to the correlated metal. While for large T and U the local approximation gives a good description of the electronic properties, at lower temperatures and close to -U c,non loc (green line on the left in figure 3), the behaviour is controlled by the inter-site correlations. In the single-orbital case at half-filling and with AFM correlations, the entropy of the insulator decreases due to short-range singlet formation [39]. As a consequence, the critical interaction -U c,non loc reverses its slope as a function of T with respect to the local-correlations prediction [39,48,52] (compare blue and green lines in figure 3). Hence, the low T phase is the insulator, as observed experimentally. Similar behaviour has been found in a two-orbital model at zero and finite J H [40]. Cluster DMFT (CDMFT) calculations in the paramagnetic two-orbital square lattice at x=1/4 and J H =0 [44] have also found that the AFM correlations are suppressed by temperature. With increasing T the crossover between the metal and the insulator is non-monotonic, see figure 3, and it matches the crossover line found in the single-site approximation at high temperatures (both dashed lines merge). The qualitative behaviour when non-local correlations are included is also found for other lattices. In particular, it is expected to be found in the honeycomb lattice, bipartite as the square one. This has been already studied in the single-orbital case at halffilling [39,49]. A larger degree of magnetic frustration decreases the magnitude of the effect, i.e.
-U c,non loc becomes closer to the local predictions and the range of temperature in which non-local correlations control the experimental behaviour gets narrower [53]. We emphasize that the temperature dependence compatible with experimental observations is found only very close to -U c,non loc , at interactions for which the system is metallic in the absence of inter-site correlations (namely, below the green line in figure 3). This result suggests that intersite correlations play a key role in the localization and insulating states observed experimentally in TBG.
The gap size controversy can also be accounted for by the inclusion of non-local correlations. While for local correlations the quasiparticle peak disappears at the Mott transition at U c,local and the gap is the one between the Hubbard bands, for < < -U U U c c ,non loc ,local (namely, between the green and blue lines in figure 3) non-local correlations open a small gap in the quasiparticle peak for smaller interactions [39,54]. Close to the transition this gap is much smaller than the one between the Hubbard bands and its size would be consistent with experimental observations. Close to -U c,non loc , in the region of parameters where the temperature dependence and gap size place the experimental system, the insulating state is controlled by the inter-site magnetic and orbital correlations. If these correlations are AFM, as it happens for zero or negative Hund [42,44], a sufficiently large magnetic field could suppress them and induce a transition to a metallic state, as observed experimentally. A magnetic field does not suppress ferromagnetic correlations. As a consequence if the inter-site correlations were ferromagnetic a magnetic field could not produce an insulator to metal transition, in disagreement with experiments. This indicates that the inter-site correlations which assist the Mott transition have antiferromagnetic character.
Finally, we note that non-local correlations can alter the ratio, discussed above, between the critical interactions at quarter-and half-fillings. The critical interaction is controlled by the non-local correlations which have a different effect at each filling, as it becomes evident by their different ordering tendencies. For example, at J H =0 the ground state of the honeycomb lattice is a valence bond solid [55] at x=1/2 but a quantum spin-orbital liquid [56,57] at x=1/4 in the strong coupling limit. Specific calculations should address whether the inequality ,non loc can be achieved for zero or negative Hund's coupling in the honeycomb lattice.

Conclusion and discussion
We have explored the phenomenology of the insulating states of twisted bilayer graphene with ±2 electrons with respect to the charge neutrality point. We have shown that the experimental observations dependences are not in agreement with the predictions for a non-ordered Mott state if only local correlations are included in the description but could be consistent once non-local correlations are taken into account. In particular, the insulator to metal transition with a magnetic field or with temperature and the small gap size are in contradiction with expectations for Mott insulating states driven by local correlations. Even in the absence of magnetic or orbital order, the inclusion of inter-site correlations reduces the critical interaction for the transition, reverses the dependence of the metal-insulator transition with temperature and strongly reduces the gap size close to the transition. A magnetic field can suppress the inter-site correlations, and consequently the insulating behaviour, if the former have antiferromagnetic character, as for  ,non loc , necessary to reproduce the metallic character of the charge neutrality point while the system is insulating at quarter-filling, is fulfilled in the local approximation, but could be affected by non-local effects. Whether it is satisfied in the presence of inter-site correlations should be studied in future work. An interesting aspect of non-local correlations is the possibility of finding pseudogap physics. Local correlations result in a momentum-independent self-energy [39]. Non-local correlations make the selfenergy momentum dependent. These momentum dependent correlations have been studied in the singleorbital square lattice when the system is doped away from the insulating state [53,[58][59][60][61][62][63]. The predicted momentum dependence has been observed in cuprate superconductors and it is a key signature of the pseudogap phenomenology [64].
In this manuscript we have described the TBG with a two-orbital Hubbard model on the honeycomb lattice. The predicted behaviour for the metal insulator transition with temperature and magnetic field, as well as the small gap close to the transition, are not restricted to this model. We expect that our qualitative conclusions would remain valid in other models proposed for the TBG which feature Mott transitions at the experimental fillings. On the other hand, the range of interaction parameters, temperature and magnetic field in which the inter-site correlations control the experimental behaviour, the ratio --U U c c 1 4 ,non loc 1 2 ,non loc , or the possible pseudogap physics will depend on the model. The description of the insulating states with ±1 and ±3 electrons from the CNP, recently observed [2] and not addressed in this work, could require different explanations depending on the lattice model. For instance, in a two-orbital honeycomb model these states could correspond to Wigner-Mott insulating states in an extended Hubbard model [65].
In summary, our work reconciles the experimental observations for the insulating states with ±2 electrons with the physics of Mott insulators, in agreement with expectations from the exact commensurability. The necessity to go beyond the local approximation to reproduce the experiments suggests that the inclusion of inter-site correlations is not just a more precise way to describe the electronic properties but a key ingredient in the localization of the carriers. The behaviour in a magnetic field indicates that these short-range non-local correlations have antiferromagnetic character.
to zero. If the orbitals are completely spin polarized the system becomes equivalent to a single orbital model at x=1/2 with interactions = U U 1orb for J H =0. Except for very small H the spin up band is emptied at interactions < U U pol 1 2 1orb , as it can be seen in figure A1(a). Beyond U pol the evolution of  Z follows the one expected in a single orbital model and ( ) U H c 1 4 tends to U 1 2 1orb , see figure A1(b) and figure 2  ORCID iDs M J Calderón https:/ /orcid.org/0000-0001-5827-4731 E Bascones https://orcid.org/0000-0003-2170-4472