Electronic band structure and pinning of Fermi energy to van Hove singularities in twisted bilayer graphene: a self consistent approach

TL;DR

This paper investigates how long-range interactions influence the band structure and van Hove singularities in twisted bilayer graphene, revealing a Fermi energy pinning effect that could impact electronic properties and manipulation strategies.

Contribution

It introduces a self-consistent approach to study the impact of long-range interactions on the band structure and Fermi energy pinning in twisted bilayer graphene.

Findings

Van Hove singularities remain close to the Fermi energy over a broad filling range.

Long-range interactions significantly alter the band dispersion and shape.

Fermi energy pinning occurs due to the sensitivity of the band structure to interactions.

Abstract

The emergence of flat bands in twisted bilayer graphene leads to an enhancement of interaction effects, and thus to insulating and superconducting phases at low temperatures, even though the exact mechanism is still widely debated. The position and splitting of the flat bands is also very sensitive to the residual interactions. Moreover, the low energy bands of twisted graphene bilayers show a rich structure of singularities in the density of states, van Hove singularities, which can enhance further the role of interactions. We study the effect of the long-range interactions on the band structure and the van Hove singularities of the low energy bands of twisted graphene bilayers. Reasonable values of the long-range electrostatic interaction lead to a band dispersion with a significant dependence on the filling. The change of the shape and position of the bands with electronic filling implies that the van Hove singularities remain close to the Fermi energy for a broad range of fillings. This result can be described as an effective pinning of the Fermi energy at the singularity. The sensitivity of the band structure to screening by the environment may open new ways of manipulating the system.

Equations

$$L=\frac{a}{2\sin \theta /2}$$

$${\mathbf{G}}_1=2\pi (1/\sqrt{3},1)/L\text{ and }{\mathbf{G}}_2=4\pi (-1/\sqrt{3},0)/L,$$

$$\Psi ={\left({\psi }_{A_1},{\psi }_{B_1},{\psi }_{A_2},{\psi }_{B_2}\right)}^T,$$

$$H_{TBG}=\int d^2\mathbf{r}{\Psi }^{†}(\mathbf{r})\left(\begin{array}{cc}{H}_1& U(\mathbf{r})\\ U^{†}(\mathbf{r})& H_2\end{array}\right)\Psi (\mathbf{r}).$$

$$H_l=\hslash v_F{\mathbf{\sigma }}_{{\theta }_l/2}\cdot \left(-i\nabla -K_l\right)$$

$$U(0)$$

$$U\left(-{\mathbf{G}}_1\right)$$

$$U\left(-{\mathbf{G}}_1-{\mathbf{G}}_2\right)$$

$$V_H$$

$$v_H(\mathbf{r})$$

$$\delta \rho (\mathbf{r})=\frac{1}{V_c}\sum _{\mathbf{G}}\delta \rho (\mathbf{G})e^{i\mathbf{G}\cdot \mathbf{r}},$$

$$\delta \rho (\mathbf{G})=4{\int }_{\text{mBZ}}\frac{d^2\mathbf{k}}{V_{\text{mBZ}}}\sum _{{\mathbf{G}}^{\prime }i,l}{\varphi }_{l\mathbf{k}}^{i,\ast }\left({\mathbf{G}}^{\prime }\right){\varphi }_{l\mathbf{k}}^i\left({\mathbf{G}}^{\prime }+\mathbf{G}\right),$$

$$\sum _{i\mathbf{G}}{\varphi }_{l\mathbf{k}}^{i,\ast }\left(\mathbf{G}\right){\varphi }_{l^{\prime }\mathbf{k}}^i\left(\mathbf{G}\right)={\delta }_{l,l^{\prime }}.$$

$$v_H(\mathbf{r})=V_0\delta {\rho }_G\sum _i{e}^{i{\mathbf{G}}_i\cdot \mathbf{r}},$$

$$\delta {\rho }_G\equiv \frac{1}{6}\sum _i{\int }_{V_c}{d}^2\mathbf{r}\delta \rho (\mathbf{r})e^{-i{\mathbf{G}}_i\cdot \mathbf{r}}$$

$$V_F$$

$$⟨{\Psi }^{i,†}(\mathbf{r}){\Psi }^j({\mathbf{r}}^{\prime })⟩={\int }_{\text{mBZ}}\frac{d^2\mathbf{k}}{{\left(2\pi \right)}^2}\sum _{\mathbf{G},{\mathbf{G}}^{\prime }}{e}^{i\mathbf{k}\cdot ({\mathbf{r}}^{\prime }-\mathbf{r})}{e}^{i{\mathbf{G}}^{\prime }\cdot {\mathbf{r}}^{\prime }}{e}^{-i\mathbf{G}\cdot \mathbf{r}}{C}_{\mathbf{k}}^{ji}\left({\mathbf{G}}^{\prime },\mathbf{G}\right).$$

$$C_{\mathbf{k}}^{ji}\left({\mathbf{G}}^{\prime },\mathbf{G}\right)\equiv \sum _{l_{occ}}{\varphi }_{l\mathbf{k}}^j({\mathbf{G}}^{\prime }){\varphi }_{l\mathbf{k}}^{i,\ast }(\mathbf{G}),$$

$$V_F=-{\int }_{\text{mBZ}}{d}^2\mathbf{k}\sum _{ij}\sum _{\mathbf{G},{\mathbf{G}}^{\prime }}{\Psi }_{\mathbf{k}}^{j,†}({\mathbf{G}}^{\prime }){\Psi }_{\mathbf{k}}^i({\mathbf{G}}^{\prime }+\mathbf{G})v_{F,\mathbf{k}}^{ji}\left({\mathbf{G}}^{\prime },{\mathbf{G}}^{\prime }+\mathbf{G}\right).$$

$$v_{F,\mathbf{k}}^{ji}\left({\mathbf{G}}^{\prime },{\mathbf{G}}^{\prime }+\mathbf{G}\right)\equiv -\sum _{{\mathbf{G}}^{\prime \prime }}{\int }_{\text{mBZ}}\frac{d^2\mathbf{p}}{V_{\text{mBZ}}}\frac{v_c\left(\mathbf{p}-\mathbf{k}+{\mathbf{G}}^{\prime \prime }\right)}{V_c}{C}_{\mathbf{p}}^{ji}\left({\mathbf{G}}^{\prime \prime }+{\mathbf{G}}^{\prime },{\mathbf{G}}^{\prime \prime }+{\mathbf{G}}^{\prime }+\mathbf{G}\right).$$

$$v_{\text{scr}}(\mathbf{q})$$

$$v_{\text{scr}}(\mathbf{G})$$

$$\chi (\mathbf{q})$$

$$\delta ϵ$$

$$\chi (\mathbf{G})$$

Full text

I Introduction

The recent observation of insulating behavior at charge neutrality (CN) and unconventional superconductivity close to CN in twisted graphene bilayer (TBG) at the magic angle $\theta\sim 1^{\circ}$ [1, 2], see also [3, 4, 5], has motivated a renewed interest in the physics of twisted graphene superlattices and related structures.

Even more recent experiments have revealed unexpected features of magic-angle TBG close to CN, such as the presence large spectral gaps at fractional filling [6], the appearance of alternating superconducting and insulating phases on an increase of the occupation of the conduction bands [7] as well as evidence of charge ordering and broken rotational symmetry [8, 9].

From a theoretical point of view, the low energy continuum model of TBG originally proposed in [10] and developed in [11, 12, 13, 14] leads to isolated flat bands at CN, with vanishing Fermi velocity at the certain angles, the “magic angles”. In contrast to tight-binding models [15, 16, 17], which require a substantial numerical effort because of the large size of the moiré superlattices at small twist angles (with more than $10^{4}$ atoms per supercell), such continuum models are a much more efficient alternative. They have recently been generalized to describe the electronic spectrum of related twisted materials, such as twisted graphene tri-layer [18, 19], quadruple-layer [20, 21, 22] and multi-layer graphitic structures [23, 24, 25].

Near the magic angle, the quenching of the kinetic energy induced by non dispersive (flat) bands can strongly enhance the role of the electronic interactions, especially if the Fermi energy is close to the van Hove singularity. A number of recent theoretical studies have addressed this topic [26, 27, 28, 29, 30, 31, 32, 33, 34, 35], suggesting this can lead o the formation of novel strongly-correlated phases. The proliferation of these singularities, due to the increasing number of band-crossings in the flat bands as we decrease the twist angle in TBG [34, 9, 36, 37], can further enhance van Hove peaks in the DOS. In addition, the strong filling dependence of the electronic compressibility observed in Ref. [6] and the experimental evidence of filling-induced symmetry-broken phases [7] emphasize the role of the carrier density in driving the electronic interactions. Note, finally, that the peaks due to van Hove singularities in the density of states can be observed directly by scanning tunneling spectroscopy[38, 9, 8, 39]

As shown in Ref. [40], the long-range Coulomb interactions in TBG away from CN may be comparable or larger than the bandwidth, giving rise to an inhomogeneous Hartree potential which substantially modifies the shape of the electronic bands. The dominant role of the Coulomb interaction in the low energy physics of the TBG near the magic angles has also been discussed out in Refs. [9, 37, 8].

In this work we report on a fully self-consistent Hartree calculation of the electronic bands of TBG near the magic angle, upon varying the filling of the conduction bands away from CN. In agreement with the findings of the Ref. [40], we find that the shape, width, and position of the bands near the neutrality point depend significantly on the filling of the same bands. Remarkably, we find an approximate pinning of the DOS at the Fermi energy, where the self-consistent bands display the tendency to flatten. As mentioned above, the proximity of the van Hove singularity to the Fermi level can lead to enhanced electron-electron interactions. The results reported here suggest that a fine tuning of the carrier density is not required for achieving an enhanced density of states at the Fermi level. Lattice relaxation leads to the reduction in size of the $AA$ regions [41, 42]. It is taken into account approximately in our calculations by assuming the interlayer hopping between $AA$ and $BB$ sublattices is lower than the hopping between $AB$ and $BA$ sublattices [14].

Our findings are valid in a wide range of values for the carrier density and the dielectric properties of the substrate. We intended to study the broken symmetry phases[43] in a future work. In this paper we make a comprehensive study of the effects of long-range electrostatic interactions, which is a crucial prerequisite to the understanding of more complex phases.

The work is organized as follows: in Sec. II we provide a technical description of the model used, introducing an inhomogeneous Hartree potential within the continuum model of TBG. In Sec. III we show and comment the numerical results concerning the self-consistent electronic bands and the corresponding DOS for different fillings and twist angles. We also provide an exhaustive analysis of the shape of the Fermi surface and of the positions of the van Hove singularities in the momenta space. In this section we study situations where the electrostatic potential is strongly screened, in order to show better the continuity between non interacting and interacting results. In Sec. IV we consider the role of the exchange potential. We show that its contribution is sub-leading when compared to the Hartree term and, more importantly, has little effect on the shape of the bands and on the position of the van Hove singularity. In Sec. V we extend the results to cases with more realistic and smaller screening. We find a substantial enhancement of the bandwith. The total bandwidth is further studied in Sec. VI were we complete the analysis by adding the effect of possible residual strains, see Ref. [44]. These combined effects might account for a more precise quantitative description of the spectra of TBG. The appendices complete the information in the main text by analyzing further: i) screening by metallic gates, ii) screening by high-energy bands in the bilayer iii) extension of the calculations to a simplified model with infinitely narrow bands at neutrality, where the dispersion is solely determined by electrostatic effects, and iv) a description of the optical conductivity and the contribution from the van Hove singularities, which could be used to study this problem experimentally.

II The model

We model the TBG at small twist angles within the low energy continuum model considered in Refs. [10, 11, 12, 13, 14]. At angles $\theta$ where the two layers are commensurate, the twist between neighboring layers leads to the formation of a moiré superlattice. The size of the corresponding unit cell

[TABLE]

dramatically increases as the angle is reduced, where $a=2.46$Å is the lattice constant of monolayer graphene. The continuum models are meaningful for small enough angles, where any twist forms approximately commensurate structures. The moiré mini Brillouin zone (mBZ), resulting from the folding of the two BZs of each monolayer (see Fig. 1(a)), has the two reciprocal lattice vectors

[TABLE]

shown in green in Fig. 1(b).

For small angles of rotation the coupling between different valleys of the two monolayers can be safely neglected, as the interlayer hopping has a long wavelength modulation. In what follows we describe the model for the behavior of the two $K$-valleys of the twisted monolayers, where $K=4\pi(1,0)/3a$ is the Dirac point of the unrotated monolayer graphene. The case corresponding to the opposite valleys at $K^{\prime}=-K$ directly follows from time reversal symmetry $\boldsymbol{k}\rightarrow-\boldsymbol{k}$.

The fermionic field operators of the TBG are 4-component Nambu spinors:

[TABLE]

where $A,B$ and $1,2$ denote the sub-lattice and layer indices, respectively. We introduce a relative twist $\theta$ between the two monolayers by rotating the layer $1$ by $-\theta/2$ and the layer $2$ by $\theta/2$. Without loss of generality, we assume that in the aligned configuration, at $\theta=0$, the two layers are $AA$-stacked. In the continuum limit, the effective Hamiltonian of the TBG can be generally written as [10, 11, 14]:

[TABLE]

Here

[TABLE]

is the Dirac Hamiltonian for the layer $l$, where $v_{F}=\sqrt{3}ta/(2\hbar)$ is the Fermi velocity, $t$ is the hopping amplitude between $p_{z}$ orbitals located at nearest neighbors carbon atoms, $\theta_{1,2}=\mp\theta$, $\boldsymbol{\sigma}_{\pm\theta/2}=e^{i\sigma_{z}\theta/4}\boldsymbol{\sigma}e^{-i\sigma_{z}\theta/4}$, $\sigma_{i}$ are the Pauli matrices, and $K_{1,2}=4\pi\left(\cos\theta/2,\mp\sin\theta/2\right)/(3a)$ are the Dirac points of the two twisted monolayers, which identify the corner of the mBZ shown in Fig. 1(a). The inter layer potential, $U(\boldsymbol{r})$, is a periodic function in the moiré supercell. In the limit of small angles, its relevant harmonics are given by the three reciprocal lattice vectors $\boldsymbol{G}=0,-\boldsymbol{G}_{1},-\boldsymbol{G}_{1}-\boldsymbol{G}_{2}$ [10]. The corresponding amplitudes $U\left(\boldsymbol{G}\right)$ are given by:

[TABLE]

In the following we adopt the parametrization of the TBG given in the Ref. [14]: $\hbar v_{F}/a=2.1354$ eV, $g_{1}=0.0797$ eV and $g_{2}=0.0975$ eV. The difference between $g_{1}$ and $g_{2}$, as described in Ref. [14], accounts for the corrugation effects where the interlayer distance is minimum at the $AB/BA$ spots and maximum at $AA$ ones, or can be seen as a model of a more complete treatment of lattice relaxation [42]. The Hamiltonian Eq. (4) then hybridizes states of layer $1$ with momentum $\boldsymbol{k}$ close to the Dirac point with the states of layer $2$ with momenta $\boldsymbol{k},\boldsymbol{k+G_{1}},\boldsymbol{k+G_{1}+G_{2}}$. For the parametrization used here we find a pair of flat bands, with vanishing Fermi velocity and no dispersion near the Fermi level for a magic angle $\theta=1.08^{\circ}$.

Next, we take into account the long-range Coulomb interaction by adding an additional Hartree term in the Hamiltonian of the Eq. (4). As emphasized in Ref. [40] by means of a Wannier-based analysis (see also [45, 14, 46]), the inhomogeneous charge distribution in TBG leads to an electrostatic potential of the order of $V_{\text{Coulomb}}\sim e^{2}/\epsilon L$, where $e$ is the electron charge and $\epsilon$ the dielectric constant of the substrate. At angles $\theta\sim 1^{\circ}$, $V_{\text{Coulomb}}$ turns out to be generally much bigger than both the mini-band width, which is a few meV, and the intra-atomic Hubbard interaction, and consequently strongly modifies the band structure close to the Fermi level. In general the Hartree contribution to the Hamiltonian is given by

[TABLE]

Here $v_{C}(\boldsymbol{r})=\frac{e^{2}}{\epsilon\left|\boldsymbol{r}\right|}$ is the Coulomb potential and $\delta\rho(\boldsymbol{r})$ is the fluctuation of the charge density with respect to CN. It is worth noting that, in the continuum model, the charge density does depend on the cutoff and, consequently, a convention has to be applied. In what follows we impose the condition that the Hartree potential vanishes at CN, where $\delta\rho=0$ by definition. The density of carriers in TBG is periodic in the moiré unit cell, with maxima in the $AA$ spots and minima in the $AB/BA$ ones [12, 13, 40, 47]. Thus it can be generally written as the superposition of plane waves:

[TABLE]

where $V_{c}=L^{2}\sqrt{3}/2$ is the area of the unit cell and the $\boldsymbol{G}$’s are reciprocal lattice vectors. The amplitudes $\delta\rho(\boldsymbol{G})$ are given by:

[TABLE]

where $V_{\text{mBZ}}=(2\pi)^{2}/V_{c}$ is the area of the mBZ, $l$ is the band index resulting from the Bloch diagonalization of the full Hamiltonian, $H_{TBG}+V_{H}$. The parameter $i=A_{1},B_{1},A_{2},B_{2}$ labels the sub lattice and layer and $\phi^{i}_{l\boldsymbol{k}}\left(\boldsymbol{G}\right)$ is the amplitude for an electron to occupy a state with Bloch momentum $\boldsymbol{k+G}$ in the $l$-th band. The $\phi$’s are normalized as

[TABLE]

In Eq. (10), the sum over the band index $l$ extends only over energy levels with $E_{l}(\boldsymbol{k})$ between the CN point and the Fermi level. The factor of $4$ accounts for the spin/valley degeneracy of the TBG Hamiltonian. The Hartree potential of the Eq. (8) implicitly depends on the fractional filling of the conduction band: $n\in[-4,4]$, where $n=4$ corresponds to fully filled valence and conduction bands, while they are both empty for $n=-4$.

As shown in Ref. [40], the dominant contribution to $\delta\rho$ comes from the first star of reciprocal lattice vectors, $\boldsymbol{G}=\boldsymbol{G}_{i}=\pm\boldsymbol{G}_{1},\pm\boldsymbol{G}_{2},\pm\left(\boldsymbol{G}_{1}+\boldsymbol{G}_{2}\right)$. As a consequence, we will neglect higher order harmonics in the following. In addition, exploiting $C_{6}$ symmetry we can assume that the Fourier components of $\delta\rho$ are equally weighted in the first star of vectors $\boldsymbol{G}_{i}$. Under these assumptions, the Hartree potential of the Eq. (8) takes the form:

[TABLE]

where $V_{0}=v_{C}\left(\boldsymbol{G}_{i}\right)/V_{c}$, with $v_{C}\left(\boldsymbol{G}_{i}\right)=2\pi e^{2}/\left(\epsilon\left|\boldsymbol{G}_{i}\right|\right)$ the Fourier transform of $v_{C}$ evaluated at $\boldsymbol{G}_{i}$, and thus $V_{0}=e^{2}/\left(\epsilon L\right)$. The real quantity

[TABLE]

is the single parameter which defines the Hartree potential, and thus needs to be determined self-consistently. Note that the Hartree potential (12) contributes to the Hamiltonian of the TBG as a diagonal term in the sub-lattice/layer subspace, coupling the generic state of momentum $\boldsymbol{k}$ to the six states of momentum $\boldsymbol{k}-\boldsymbol{G}_{i}$, with an amplitude equal to $V_{0}\delta\rho_{G}$.

For a given value of the filling of the conduction band, $n$, we compute the valence and conduction bands by diagonalizing the Hamiltonian for the TBG in the presence of the Hartree potential of Eq. (12), $v_{H}$, in which $\delta\rho_{G}$ is obtained as the self-consistent solution of the Eq. (13). This procedure is iterated until convergence. Numerical details concerning the iterative method are provided in Appendix A.

III Bands distortion induced by the inhomogeneous Hartree potential

In this section we show the numerical results for the self-consistent diagonalization described in the previous section, as a function of the twist angle and the filling $n$.

The self-consistent order parameter, $\delta\rho_{G}$, is a dimensionless number varying between $-4$ and $4$. As shown in Ref. [40] and in Fig. 2 below, the value of $\delta\rho_{G}$ depends approximately linearly on $n$ in a wide range of fillings and parameters. For the results in this section, we use a very high value of the dielectric constant $\epsilon\simeq 66$. to reduce the strength of the Coulomb force. For this choice we find that the amplitude of the fluctuations of the Hartree potential between the $AA$ and $AB/BA$ spots of the moiré unit cell, $V_{H}(AA)-V_{H}(AB)\simeq 5\,\mathrm{meV}$, for a band filling $n=1$. This energy scale is comparable to the bandwidth of the TBG in the absence of interactions, so that a simple one-to-one comparison between the two cases is possible.

A study of the effects of more realistic electrostatic interactions that include a more sensible screening parameter, $\epsilon\sim 4-8$, is carried out in Sec. V. The details of the screening due to external metallic gates, and to the high-energy subbands in the bilayer is discussed in Appendix A.

Figure 3 shows the bands obtained along the high-symmetry circuit of the mBZ, as shown in Fig. 1(b), for various fillings and angles near the magical one. As we anticipated, the bands do not rigidly shift when varying the filling, but they rearrange their shape in a quite non trivial way. Most interestingly they tend to flatten in the vicinity of the Fermi level. We note that the bands shown in Fig. 3b at CN (red line) behave rather differently than those in the other panels. This is because $\theta=1.08^{\circ}$ is a magic angle for our choice of parameters, so that the corresponding low energy bands minimize the bandwidth as compared to the other angles considered in Fig. 3. Also the appearance of new Dirac points is a typical feature of the bands of TBG for angles near the magic one.

This is confirmed by the DOS, shown in Fig. 4, which is clearly pinned at the Fermi energy for $n\neq 0$.

Moreover, the reshaping of the bands is much more evident around the points $K$ and $M$ than at $\Gamma$, where the bands are not sensitive to a change in the density of carriers.

In Fig. 5 we can see the details of the relation of the position of the Fermi surface and the van Hove singularities in the mBZ, for the angles $\theta=1.08^{\circ},1.05^{\circ}$ and fillings $n=1,2$. As already suggested by the results of Figs. 3,4, the van Hove singularities generally match the Fermi surface closely for many parameters. This effect is more evident for $\theta=1.08^{\circ}$ and $n=1$, and for $\theta=1.05^{\circ}$ and $n=2$, than for the other two cases.

In Fig. 6 we show the value of the DOS at the Fermi level as a function of the filling of the conduction band, for the angles $\theta=1.08^{\circ}$ and $\theta=1.05^{\circ}$. In both cases the DOS displays intense isolated peaks. This feature suggests the possibility of enhancing or suppressing the electron-electron interactions in TBG by fine tuning the number of charge carriers, as recently pointed out in experiments [6, 7, 9, 8]. Furthermore, the features of the DOS are generally not particle-hole symmetric, which agrees with the experimental findings of Ref. [6].

To gain a more detailed understanding of the evolution of the band structure with filling, and to stress the main differences between the interacting and non-interacting case, we show the bands (Fig. 7(a)) and the DOS (Fig. 7(b)) obtained at the magic angle $\theta=1.08^{\circ}$, for the fillings $n=-1,0,1$, where $n=1$ represents the non-interacting case.

IV Role of the exchange potential

If we wish to perform a fully consistent mean-field calculation, we need to include the Fock term. The Fock or exchange potential can be written as

[TABLE]

where the brackets, $\left\langle\dots\right\rangle$, denote the expectation value in the filled states of the full Hamiltonian $H_{TBG}+V_{H}+V_{F}$. Using an expansion in Bloch waves, we can write this as

[TABLE]

Here

[TABLE]

with the sum over all the occupied levels, $l_{occ}$, in the two bands in the middle of the spectrum. The Hamiltonian (IV) is more conveniently expressed in momentum space. Substituting Eq. (15) into (IV) and folding the BZ of monolayer graphene gives

[TABLE]

Here $\Psi_{\boldsymbol{k}}(\boldsymbol{G})\equiv\int\frac{\,d^{2}\boldsymbol{r}}{2\pi}\Psi(\boldsymbol{r})e^{i\left(\boldsymbol{k+G}\right)\cdot{r}}$ is the Fourier component of the field operator $\Psi(\boldsymbol{r})$, corresponding to the Bloch momentum $\boldsymbol{k}+\boldsymbol{G}$, and

[TABLE]

The off-diagonal elements, $v_{F}^{j\neq i}$, can be safely neglected, as the interactions coupling different layers and/or sub lattices are supposed to be local in space within the continuum model, while the exchange potential describes the non-local effects of the Coulomb interaction. Assuming that the Coulomb potential is negligible at momenta outside the mBZ, we retain only the term corresponding to $\boldsymbol{G^{\prime\prime}}=0$ in the sum of the Eq. (18). In addition, to link with the truncation of the Hartree potential, we only consider matrix elements of the exchange potential (17) coupling states of momentum $\boldsymbol{k+G}$ with states of momentum $\boldsymbol{k+G}$ and $\boldsymbol{k+G+G_{i}}$, where $\boldsymbol{G_{i}}$ belongs to the first star of reciprocal lattice vectors.

In Fig. 8 we compare the DOS obtained using the Hartree-Fock approximation (continuum lines) to those with only Hartree potential (dotted lines), for a twist angle $\theta=1.05^{\circ}$ and different fillings of the conduction band. The main effect of the exchange potential is visible at CN, $n=0$, where the Hartree potential is absent, and at full filling, $n=4$, where the exchange potential is enhanced by the contribution of all the electronic states. However, in both these cases the difference between the results corresponding to the full Hartree-Fock calculations and those obtained with only the Hartree potential is not qualitatively significant, even if the DOS differs quantitatively in the two cases. In panel (b), for two electrons per unit cell in the conduction bands (including valley and spin), the Fock contribution is negligible when compared to the Hartree one. For empty bands, panel (c), the exchange potential is zero. We can conclude that adding the exchange potential only makes minor modifications to the DOS or the pinning at the Fermi level, but does not give a qualitative change to the main results of this work. We conclude that the main effect causing distortion of the bands originates from the Hartree potential.

V Effect of screening

As we mentioned above, the results shown in the previous sections correspond to a strong screening, in order to facilitate the comparison between the non-interacting and the interacting bands. The screening from metallic gates, and from the higher energy bands of the bilayer, at the relevant wavevector for the electrostatic potential, $|\boldsymbol{G}|=(4\pi)/(\sqrt{3}L_{M})$, is not expected to be very large. Typical values for the dielectric function of the hBN slabs that encapsulate the bilayer are $\epsilon\approx 4-8$ (note that, in addition, a very thin hBN slab does not screen electrostatic potentials at low wavevectors, see Appendix B. In this section, we analyze the role of the Hartree potential for screening constants of $\epsilon=4$ and $\epsilon=7.5$.

Fig. 9 shows the bands (left panels) and the DOS (right panels) obtained for $\epsilon=7.5$ (a) and $\epsilon=4$ (b), for $n=-1,0,1$ and twist angle $\theta=1.08^{\circ}$. Analogous results corresponding to $\theta=1.05^{\circ}$ are shown in Fig. 10. The distortion of the bands induced by the Hartree potential is much stronger than for the results shown in Sec. III. Though, the data still display a clear pinning of the DOS at the Fermi level and the main qualitative features pointed out in the previous sections remain nearly unchanged. In addition, the overall bandwidth sensitively grows as $\epsilon$ is lowered. For the case corresponding to $\epsilon=4$, which is a quite realistic value, we obtain a bandwidth of approximately $30\,\mathrm{meV}$. This value is much larger than what predicted by the continuum models [10, 11, 12, 13, 14] ($\lesssim 5\,\mathrm{meV}$). Consequently, we argue that the Coulomb interaction will make an important contribution in broadening the mini-band width of the TBG near the magic angle.

VI Effect of residual strains

For completeness, we consider the additional effects of residual strains, which may play a role in realistic situations. The presence of strain is revealed by the elliptical shape of the moiré dots observed in STM experiments [38]. As shown in [41, 44], the combined effects of relaxation and strain might explain the large distance between the van Hove singularities, up to $55\,\mathrm{meV}$, observed in [38, 6]. Here we consider a uniaxial strain of $\varepsilon=0.6\%$, applied along the direction forming an angle of $30^{\circ}$ with respect to the horizontal axis. This leads to a distortion of the mBZ as shown in Fig. 11(b). For the technical details concerning the parametrization of the strain in the continuum model of TBG, we refer to Ref. [44]. Figure 11(a) shows the DOS obtained including strain and a dielectric constant $\epsilon=4$, for the twist angle $\theta=1.05^{\circ}$ and different fillings of the conduction bands, as labeled by the vertical right axis. The overall bandwidth is delimited by the two peaks associated to the van Hove singularities, with a distance between them of approximatively $\sim 35\,\mathrm{meV}$, in agreement with the experimental observations [38, 6]. In general, we find a less pronounced pinning of the DOS at the Fermi energy as compared to the strain-less case of the Sec. III. Interestingly, moving away from CN, the peak closest to the Fermi energy broadens while the farthest one gains intensity.

VII Conclusions

We have analyzed the effect of the electrostatic interaction on the band structure of a TBG near the magic angle $\theta\sim 1^{\circ}$, with emphasis on the width of the bands, and the change in the role of the van Hove singularities in the density of states. The scale of the electrostatic potential, $\sim e^{2}/\epsilon L_{M})$, is comparable to or larger than the bandwidth $W$ of the non-interacting system. Even for an unrealistically large screening, $\epsilon\sim 50-80$, the scale of the Coulomb force equals the band splitting near the first magic angle, suggesting that such effect may drive the electronic properties of such systems.

An analysis of the screening by the hBN substrate, metallic gates, and high energy bands of the twisted bilayer allows us to estimate that realistic values for the screening of the electrostatic potential are in the range $\epsilon\sim 6-10$. For these values of the screening, the bandwidth is mostly determined by the Hartree potential–the exchange terms have little effect, apart for CN where the Hartree terms vanish. The strong sensitivity of the band structure to the screening by the environment opens new ways of manipulating the system.

The effect of the Hartree potential varies significantly as function of filling. This dependence of the shape of the bands on electronic filling shifts states from above (below) the van Hove singularities to below (above) them, depending on whether the carriers are electrons or holes. This leads to an approximate pinning of the van Hove singularities at the Fermi energy at certain fillings, in agreement with experimental results.

This suggest that interactions stabilize such pinning, making it more stable over a wider range of fillings, and thus much more amenable to experimental measurement and exploitation in future devices.

After the completion of this work[48], a related manuscript has been posted[49]. This work analyzes the role of the Hartree potential as function of band filling, using a real space tight binding description of a twisted graphene bilayer at a magic angle. As far as the two papers overlap, the results are consistent.

Acknowledgments. We would like to thank José González and Tobias Stauber for useful conversations. This work was supported by funding from the European Commission under the Graphene Flagship, contract CNECTICT-604391. NRW acknowledges support by the UK STFC under grant ST/P004423/1.

Appendix A Numerical calculation of the order parameter

As we already mentioned in the main text, the order parameter $\delta\rho_{G}$, which defines the Hartree potential, Eq.  (12), is determined as the self-consistent solution of Eq. (13). We compute the value of $\delta\rho_{G}$ numerically by means of an iterative procedure starting with an initial value, which in general we choose as $1$. At the $m$-th step of this procedure, we compute the value $\delta\rho^{(m+1)}_{G}$ by diagonalizing the Hamiltonian of the TBG in the presence of an Hartree potential given by $\delta\rho^{(m)}_{G}$. We set the threshold for convergence by imposing $\left|\delta\rho^{(m+1)}_{G}-\delta\rho^{(m)}_{G}\right|<10^{-12}$, where $\delta\rho^{(m)}_{G}$ is dimensionless and generally of $O(1)$. In the majority of the cases, the algorithm converges in less than 10 steps. As an example, in Fig. 12 we show the values of $\delta\rho^{(m)}_{G}$ versus the number of steps, $m$, for the angle $\theta=1.05^{\circ}$ and filling $n=-1$ (cyan), $n=1$ (blue) and $n=2$ (green). As can been seen, there are no significant fluctuations of the order parameter while reaching the solution, which ensures the stability of our method.

Appendix B screening

B.1 Screening due to metallic gates.

The screening induced by a metallic gate at distance $d$ from the twisted graphene bilayer can be analyzed using the method of image charges.A single metallic gate changes the Fourier transform of the effective interaction between two point charges in the bilayer. For a gate at distance $D$ from the bilayer, the potential becomes

[TABLE]

where $\boldsymbol{q}$ is a two dimensional vector in the plane of the bilayer. The relevant Fourier components of the Hartree potential are at reciprocal lattice vectors with $|\boldsymbol{G}|=(4\pi\sqrt{3})/(3L_{M})$ where $L_{M}$ is the moiré lattice length. Hence, the effective interaction is reduced by a factor

[TABLE]

For bilayer gate distances $D\sim L_{M}$ we obtain that the reduction factor is $1-e^{(4\pi)/\sqrt{3}}\approx 0.9993$.

A second gate on the other side of the bilayer leads to an infinite series of image charges, with alternating signs. The effect of this series can be summed. For simplicity, we assume that the top and bottom gates are at the same distance from the bilayer, $D$, and obtain

[TABLE]

For $e^{-|\boldsymbol{G}|D}\ll 1$, the effect of two gates is, approximately, twice the effect of a single gate.

B.2 Screening due to the upper and lower bands.

The polarization of the bands not considered in the main text does not contribute to the static dielectric constant, $\boldsymbol{q}\rightarrow 0$. The system formed by these bands describes a two dimensional insulator. The dependence of the susceptibility of these bands, $\chi(\boldsymbol{q},\omega)$, at $\omega=0$ and $|\boldsymbol{q}|\rightarrow 0$ is $\chi(\boldsymbol{q},\omega=0)\propto|\boldsymbol{q}|^{2}/\Delta$, where $\Delta$ is an energy scale of order of the gap of the insulator. The contribution to the static susceptibility goes as $\delta\epsilon\propto v_{\boldsymbol{q}}\chi(\boldsymbol{q},\omega)=0)\propto|\boldsymbol{q}|\rightarrow 0$. As function of momentum, $|\boldsymbol{q}|$, the screening of a slab with $N$ layers will be significant for $|\boldsymbol{q}|\gtrsim\Delta/(Ne^{2})$, provided that the total thickness of the slab, $w$, is such that $w\lesssim|\boldsymbol{q}|^{-1}$. Note that this analysis is also valid for thin slabs of hBN layers.

The relevant components of the Hartree potential in twisted bilayer graphene have Fourier components at the lowest reciprocal lattice vectors, $|\boldsymbol{G}|=(4\pi\sqrt{3})/(3L)$. We first approximate the screening at this wavevector by considering the screening due to two decoupled graphene layers. The susceptibility of undoped graphene is[50]

[TABLE]

where ${\cal N}_{f}=8$ is the number of Dirac cones. The contribution to the screening function is

[TABLE]

This result is independent of the wavevector, in our case $\boldsymbol{G}$. For isolated graphene layers, we can take $e^{2}/(\hbar v_{F})\approx 2$, and we obtain $\delta\epsilon\sim 6$.

A more precise estimate of the internal screening can be obtained from directly studying the continuum model, The screening from the upper bands at low wavevectors, $\boldsymbol{q}\rightarrow 0$, was considered in Ref. [51]. The calculation is simplest in reciprocal space. The density operator $\tilde{\rho}_{\boldsymbol{G}}$ only has finite matrix elements between states with the same momentum in different bands. The static susceptibility can thus be written as

[TABLE]

where labels $l$ and $m$ denote band indices, the factor 4 stands for the spin and valley degeneracy, and the factor 2 follows from the double counting of the occupied states in the Lindhard function. Numerical calculations give $\chi(\boldsymbol{G})\approx\alpha\,V_{c}^{-1}$ with $\alpha=0.07\,\mathrm{meV}^{-1}$ where we use the parameters described in the main text. The correction to the screening is $\delta\epsilon_{\boldsymbol{G}}\approx 8$, in reasonable agreement with the previous analytical estimates.

Appendix C Band structure and pinning in twisted graphene layers with infinitely narrow bands.

The effects discussed in this work, the change of the band structure of the low energy bands, and the approximate pinning of the van Hove singularities at the Fermi level for certain band fillings, is clearly manifested in a model where the interlayer hopping between the $AA$ lattices is set to zero[12, 47]. This model has electron-hole symmetry, a number of almost equally spaced magic angles, where the bandwidth of the lowest bands is strictly zero[12].The model gives a simpler description of the emergence of magic angles in twisted bilayer graphene, while keeping the most relevant physical features intact.

If we assume that at half filling the Hartree potential is zero throughout the moiré unit cell, we can expect that a self consistent Hartree potential will arise away from half filling. This potential induces a finite bandwidth, which is solely due to the electrostatic interaction. As mentioned in the main text, the fact that he number of electronic flavors is 4 makes the exchange term significantly smaller than the Hartree contribution.

Results for the bands and for the density of states are shown in Fig. 13. The plots in Fig. 13 indicate that the pinning of the van Hove singularities at integer fillings is enhanced when the bands are infinitely narrow at neutrality. The two bands remain degenerate away from neutrality, although they become dispersive. This result can be understood if we assume that the Wannier functions which describe the system at neutrality are the same for the bands away from neutrality. At neutrality, the flatness of the bands implies that all hoppings between Wannier functions are zero. The dispersion shown in Fig. 13 is due solely to new hoppings induced by the Hartree potential. As the model used assumes that the Hartree potential changes in the same way both layers, the two bands acquire a dispersion, but remain degenerate.

Appendix D Optical absorption.

A good alternative to STM experiments for studying the van Hove singularities in the density of states experimentally is by observing the optical absorption. In twisted bilayer graphene at the magic angles intra- and interband transitions between the two lowest bands occur at very low energies, and they may thus be difficult to observe, however.

We have performed a calculation of the optical absorption following the framework set out in [52, 53]) for the four bands and fillings shown in Fig. 3c. The results are shown in Fig. 14. It is interesting to note that that the energy of the main peak, related to transitions between van Hove singularities, hardly changes with filling. The van Hove singularities in the density of states associated to the peaks in the optical absorption are in different bands for $n=0$, and in the same bands for $n\neq 0$. Notice that the main peak occurs at $0.5\,\text{meV}$, which is of course a hard wavelength to observe.