Hartree theory calculations of quasiparticle properties in twisted bilayer graphene

A detailed understanding of interacting electrons in twisted bilayer graphene (tBLG) near the magic angle is required to gain insights into the physical origin of the observed broken symmetry phases including correlated insulator states and superconductivity. Here, we present extensive atomistic Hartree theory calculations of the electronic properties of tBLG in the (semi-)metallic phase as function of doping and twist angle. Specifically, we calculate quasiparticle properties, such as the band structure, density of states (DOS) and local density of states (LDOS), which are directly accessible in photoemission and tunnelling spectroscopy experiments. We find that quasiparticle properties change significantly upon doping - an effect which is not captured by tight-binding theory. In particular, we observe that the partially occupied bands flatten significantly which enhances the density of states at the Fermi level and explains the experimentally observed Fermi level pinning. We predict a clear signature of this band flattening in the LDOS in the AB/BA regions of tBLG which can be tested in scanning tunneling experiments. We also study the dependence of quasiparticle properties on the dielectric environment of tBLG and discover that these properties are surprisingly robust as a consequence of the strong internal screening. Finally, we present a simple analytical expression for the Hartree potential which enables the determination of quasiparticle properties without the need for self-consistent calculations.

To understand the properties of interacting electrons in tBLG, many different theoretical approaches have been used. In strong correlation techniques, such as dynamical mean-field theory [25], Quantum Monte Carlo [26,27] or exact diagonalization methods [28], an effective Hamiltonian for the flat-band electrons is often constructed by adding Hubbard-like interaction terms to a kinetic energy that describes the hopping between flat band Wannier functions [29][30][31][32]. However, the resulting Hamiltonian is relatively complicated and contains long-ranged hoppings [29]. Moreover, the accurate determination of the Hubbard parameters is difficult and the construction of flat-band Wannier functions can be hindered by obstructions [33,34].
In contrast, mean-field treatments of electron-electron interactions are conceptually more straightforward and do not require the construction of flat band Wannier functions. Cea, Walet and Guinea [35] used Hartree theory within a continuum model of tBLG to calculate band structures and densities of states (DOS) as function of doping and twist angle. They found that the band structure of doped tBLG changes significantly when electronelectron interactions are included, with results that are in qualitative agreement with recent scanning tunnelling spectroscopy (STS) studies [7][8][9][10] which showed that the Fermi level of the doped system is pinned at the van Hove singularity (VHS). Several groups have also carried out Hartree-Fock calculations of tBLG and studied brokensymmetry phases [36][37][38][39][40][41][42]. These calculations are also all based on a continuum theory for the electronic structure of tBLG. While continuum model calculations are numerically very efficient, they typically employ a shortwavelength cutoff for the plane-wave expansion [35] of the charge density and do not capture the effect of atomic scale Hubbard interactions. Klebl and Honerkamp [43] carried out atomistic calculations of the spin susceptibility of tBLG with shortranged atomic Hubbard interactions using the randomphase approximation and found that tBLG inherits magnetic properties from the untwisted bilayer [44]. To date, the first and only atomistic study of the effect of long-ranged electron-electron interactions on the electronic structure of tBLG was carried out by Rademaker, Abanin and Mellado [45] who used Hartree theory within an atomistic tight-binding model to calculate the charge density, band structure and local density of states (LDOS) in the AA-stacked region of both undoped and hole-doped tBLG at a single twist angle (θ = 1.05 • ). They find that electron-electron interactions smoothen the charge density and observe significant changes in the band structure upon doping in qualitative agreement with the continuum Hartree theory calculations of Cea and coworkers [35].
In this article, we present a systematic study of the effect of long-ranged Coulomb interactions on the band structure, DOS and LDOS as a function of twist angle and doping in tBLG near the magic angle. Specifically, we carry out self-consistent atomistic Hartree calculations. For tBLG suspended in air, we find that electronelectron interactions induce significant changes to the band structure of doped tBLG. In particular, for twist angles near (but not directly at) the magic angle, the partially occupied bands flatten while the unoccupied or fully occupied bands become more dispersive. This explains both the Fermi level pinning and the different shapes of the VHSs observed in recent STS experiments. While most STS experiments have focused on the AA regions of tBLG, we predict a significant enhancement of the peak in the AB regions and hypothesize that this peak is responsible for instabilities to broken-symmetry states even when the system is not at the magic angle. We also study the dependence of the band structure on the dielectric environment and find that the environmentinduced changes are relatively small. This is a consequence of the large internal dielectric constant of tBLG. We do not explicitly investigate broken-symmetry states in this work.

II. METHODS
We study commensurate moiré unit cells of tBLG, defined using the convention of Ref. [20]. At small twist angles, tBLG undergoes significant in-plane and out-ofplane atomic relaxations [46][47][48][49][50][51]. We calculate these relaxations using classical force fields: interactions between atoms in the same graphene layer are modelled using the AIREBO-Morse potential [52], whilst the Kolmogorov-Crespi potential [53] is used for interactions between atoms in different layers. All relaxations are carried out with the LAMMPS software package [54].
To calculate electronic properties of tBLG, we use atomistic Hartree theory and diagonalize the following Hamiltonian where ε i andĉ † i (ĉ i ) denote the on-site energy of a carbon atom and the electron creation (annihilation) operator associated with the p z -orbital on atom i, respectively. The hopping parameters t(r i −r j ) between atoms i and j are obtained using the standard Slater-Koster rules [55] with parameters derived from DFT [20,21,56]. A spherical cut-off of 20Å is employed for retaining hopping terms [57], symbolically denoted by ij in Eq. (1). The on-site energy is determined by the Hartree potential energy V H (r) according to where φ z denotes the carbon p z -orbital at position t i in the first unit cell (note that V H (r) is periodic in the moiré unit cell). The Hartree potential is obtained from the electron density n(r) and the screened electron-electron interaction W (r) via where n 0 (r) is a reference electron density that ensures overall charge neutrality. We consider two cases for the screened interaction. One is for tBLG encapsulated by a dielectric substrate with background dielectric constant bg and, hence, a screened interaction given by W (r) = e 2 /(4π 0 bg |r|). The other is for the case when there is the additional presence of metallic gates on both sides of the dielectric substrate. Assuming that the tBLG lies in the x-y plane, the screened interaction in this case is given by where ξ is the thickness of the dielectric substrate separating tBLG from the metallic gate on each side [58,59]. The charge density can be expressed in terms of the Bloch eigenstates ψ nk (r) (with subscripts n and k denoting a band index and the crystal momentum, respectively) of the Hamiltonian in Eq. (1) according to where f nk = 2Θ(ε F − ε nk ) is the occupancy of state ψ nk with eigenvalue ε nk (where ε F is the Fermi energy), χ j (r) = R φ 2 z (r − t j − R) (with R denoting the moiré lattice vectors) and n j is the total number of electrons in the j-th orbital. Note that we neglect contributions to the density that result from the overlap of p z -orbitals on different atoms.
To construct the reference electron density n 0 (r), we note that the hopping parameters of Eq. (1) were obtained from fits to band structures of graphene and untwisted graphene bilayers calculated using densityfunctional theory (DFT) and, therefore, include the Hartree potential energy of the uniform system (when the occupancy of all carbon atoms is equal) [20,21,56]. To exclude this contribution to V H in our tBLG calculations, we use the reference density wheren is the average of n j . The average filling can be expressed asn = 1 + ν/N , where ν denotes the number of electrons that have been added to or removed from the moiré unit cell and N is the total number of states (also atoms in the moiré unit cell).
To obtain a self-consistent solution of the atomistic Hartree equations, we proceed as follows. We first set ε i = 0 and diagonalize the Hamiltonian, Eq. (1), using an 8 × 8 k-point grid to sample the first Brillouin zone. From the eigenstates, we construct an initial guess for n j andn. Next, we calculate the updated on-site energies via where . This is equivalent to treating φ 2 z as a delta-function when considering interactions between different atoms. We carry out calculations for both bg = 1 (tBLG suspended in air) and bg = 3.9 (tBLG sandwiched between hexagonal boron nitride). To converge the sum over moiré lattice vectors, we use a 21 × 21 supercell.
To speed up convergence, we divide ε i by a factor of 8 in the first iteration as the converged Hartree charge density is much more uniform compared to the noninteracting tight-binding result. In each subsequent iteration of the self-consistent cycle, we mix a fraction of the new Hartree potential with the Hartree potential from the previous iteration. A mixing fraction of 0.1, i.e., the addition of 10 percent of the new potential to 90 percent of the potential from the previous iteration, was found to give satisfactory results in most cases. In a few cases, however, smaller values for the mixing fraction were used to improve convergence. Typically, the Hartree potential converges within 60 iterations to an accuracy of better than 0.1 meV for all doping levels and twist angles considered.
In order to calculate the density of states (DOS) per moiré cell, we sample the first Brillouin zone using approximately 6,000 k-points and represent the contribution from each energy level as a gaussian. A similar procedure is used for the local density of states (LDOS). Note that we average the LDOS over atoms within a radius of 15Å [we found that the results do not depend qualitatively on the radius chosen, provided it is larger than the length scale of the carbon-carbon bond length O(1Å) but smaller than the moiré length scale O(10 nm)].

III. RESULTS AND DISCUSSION
A. Band structure Figure 1 shows the band structures of tBLG suspended in air ( bg = 1) from Hartree theory at six twist angles between θ = 1.54 • and θ = 1.05 • (below, we show that a dielectric substrate only leads to small changes in the Hartree band structures). Only the four flat bands closest to the Fermi energy are shown. We refer to the lower two of the flat bands as the flat valence bands and the upper two as the flat conduction bands. These four bands are separated from all other bands by energy gaps that result from the atomic corrugation of tBLG [23,29,31,47,48,57]. The width of the flat band manifold is smallest at θ = 1.16 • and we refer to this twist angle as the magic angle.
We first discuss the band structures of undoped tBLG, corresponding to ν = 0 (black curves in Fig. 1). The band structures at all twist angles except the magic angle are semi-metallic and feature linear bands at the K and K points. As the magic angle is approached, the total width of the flat band manifold decreases rapidly, see Fig. 2(a). Interestingly, at charge neutrality, the valence band widths are always smaller than the conduction band widths, see Figs. 2(c) and (d).
At the magic angle, the band structure of undoped tBLG is qualitatively different as compared to the other twist angles [23,61]. In particular, the two valence bands at Γ are pushed up and are now higher in energy than the states at K and K . As a consequence, at this level of theory, tBLG is metallic at the magic angle even without doping with additional carriers.
The Hartree band structures of undoped tBLG are very similar to the non-interacting tight-binding band structures (see Supplementary Materials for a comparison). This can be understood by analyzing the charge density and the corresponding Hartree potential. Without doping the charge density oscillates on the atomic scale, but each region of the moiré unit cell is approximately charge neutral (when the atomic oscillations are averaged over) [45] resulting in a small Hartree potential [Eq. (3)], as we shall discuss further later. Figure 1 also shows Hartree band structures for electron-doped (upper two rows) and hole-doped (lower two rows) tBLG. In agreement with previous Hartree calculations [24,35,45], we observe that doping results in significant changes in the band structures. In contrast, the tight-binding band structures that are widely used to understand the electronic properties of doped tBLG do not change upon doping. Focusing first on the largest twist angle considered, θ = 1.54 • , electron doping (corresponding to ν = 1, 2 or 3) flattens the conduction bands, while the valence bands become more dispersive. Figure 2(c) shows that the conduction band width decreases by approximately 5 meV for each added electron. However, the valence band width increases by the same amount [ Fig. 2(d)] and the total band width of the flat band manifold remains constant at this twist angle, as shown in Fig. 2(a). For hole doping (ν = −1, −2 and −3), the situation is similar but the valence bands flatten and the conduction bands become more dispersive.
To understand why electron-electron interactions are more relevant for the doped system, we analyze again the charge density and the corresponding Hartree potential (the explanation here follows that outlined in Ref. [45]). As the local density of states is larger in the AA regions than in the AB/BA regions, additional carriers (both electrons and holes) preferentially localize in the AA regions [45]. This creates a highly non-uniform charge distribution, which gives rise a strong Hartree potential [45]. gions) increases by approximately 30 meV for each added electron. States near the K and K points are localized in the AA regions and are pushed up in energy relative to the states at Γ (which have a ring-like shape surrounding the AA regions) for electron-doped systems [45]. In contrast, the K/K states are pushed down in energy relative to the Γ-states for hole-doped systems [45]. For smaller twist angles, doping induces even more significant changes in the band structure. At θ = 1.41 • , the valence bands are almost completely flat between Γ and M for ν = −3. In contrast, the flattening of the conduction bands upon electron doping is not as pronounced at this twist angle. For θ = 1.2 • , the Γ-states have moved past the K/K -states so that the curvature of the conduction band at Γ changes sign at all doping levels (both electron and hole doping) except ν = 1. For this doping level, the conduction band is very flat in the vicinity of the Γ-point. Interestingly, for ν = 2 the width of the conduction bands exhibits a local minimum at θ = 1.2 • , see Fig. 2(c), and is even smaller than at the magic angle (defined as the twist angle that exhibits the smallest total band width of the flat band manifold, θ = 1.16 • ). Similarly, for ν = 3 the width of the conduction bands exhibits a local minimum at θ = 1.3 • . This suggests that long-ranged Coulomb interactions between electrons can modify the twist angle at which electron correlation phenomena are strongest and that this may not necessarily be at the magic angle.
These qualitative changes in the band structures of doped tBLG close to the magic angle can be understood by analyzing the twist angle dependence of the Hartree potential. Fig. 6(a) shows that ∆V H only depends weakly on the twist angle. In contrast, the band widths decrease rapidly as the magic angle is approached and therefore the importance of long-ranged electron-electron increases strongly.
At the magic angle (θ = 1.16 • ), the band structures of hole-doped tBLG (ν = −1, −2 and −3) look qualitatively similar to the undoped band structure, but with a significantly larger band width. For example, for ν = −2 we find a band width of 13 meV (compared to 1 meV for the undoped system). For electron-doped systems, the conduction bands 'invert' such that both the valence and conduction bands at Γ are at lower energies than the states at K and K .
For twist angles smaller than the magic angle, the band structures of doped tBLG are quite similar to those of twist angles larger than the magic angle. In particular, the band structures at θ = 1.12 • correspond closely to those of θ = 1.2 • (both differ from the magic angle by 0.04 • ) and the band structures of θ = 1.05 • are similar to those of θ = 1.3 • (which differ from the magic angle by 0.11 • and 0.14 • , respectively). . The LDOS is shown both for the AA (solid lines) and AB (dashdotted lines) regions, averaged over a region around the centre of the respective region (as discussed in the Methods section). When the tunnelling matrix elements are constant (which is likely a good approximation for the flat bands of tBLG), the LDOS is proportional to the measured tunnelling spectrum and thus directly accessible in experiments. Several STS studies of tBLG have been reported recently [8][9][10] and we will discuss the similarities and differences of our calculations with these experimental measurements. Below, we analyse each of the three twist angles in turn; the results for the other twist angles can be found in the Supplementary Materials.

B. DOS and LDOS
For θ = 1.54 • , the DOS exhibits two pronounced peaks at all doping levels. At zero doping, these van Hove singularities (VHSs) are located at ±20 meV on both sides of the Fermi energy (so their energy separation is approximately 40 meV). Comparing the DOS to the LDOS, we find that the dominant contribution to the DOS derives from the AA regions [20]. The LDOS in the AB regions also exhibits small peaks in the vicinity of the VHSs [9]. The valence band VHS is somewhat larger than the conduction band one because the valence bands are flatter than the conduction bands, see Fig. 1. These findings are in agreement with experimental STS measurements [8][9][10]. Note, however, that our values for the energy difference between valence and conduction VHSs are smaller (for the same twist angle) than the experimental results. In Ref. [9], this was attributed to the use of DFT-derived tight-binding parameters for the interlayer hopping which are about 20 percent smaller than more accurate GW values.
Upon electron doping, the conduction VHS becomes larger while the valence VHSs becomes smaller. This is a consequence of the doping-induced band flattening of the conduction bands, while the valence bands become more dispersive, see Fig. 1. In contrast, hole doping increases the valence VHS while the conduction VHS becomes smaller. Again, these findings are in agreement with experimental measurements [9] and cannot be ex-plained by tight-binding theory. Note that at this twist angle the Fermi level of the doped system is not pinned at the VHSs.
At θ = 1.41 • , the separation between the VHSs is reduced to 30 meV. Upon hole doping, the difference between valence and conduction band VHSs is much clearer than at 1.54 • . This is caused by the strong distortion of the doped valence bands resulting in extremely flat valence bands throughout large regions of the Brillouin zone, see Fig. 1 (recall that the distortion of the valence bands is always more pronounced that that of the conduction bands). For ν = −2 and ν = −3, we observe that the Fermi level is pinned at the valence VHS. This Fermi level pinning has also been reported in experimental STS studies and is a consequence of electron-electron interaction induced changes in the band structure [9]. The LDOS in the AA region is again very similar to the DOS. However, we find that the valence peak of the LDOS in the AB regions grows significantly upon hole doping (see SM for further details). This is because the wave functions of the flat valence bands are partly localized in the AB regions (in particular, the valence states near Γ). This prediction can be tested by STM measurements and would provide direct evidence of the doping-induced band flattening in Hartree theory. Figure 2(b) shows that the separation of the VHSs is reduced by hole doping for twist angles larger than the magic angle and increased by electron doping. The opposite trend is observed for twist angles smaller than the magic angle. While this is in qualitative agreement with experimental measurements [9], the absolute magnitude of the change in VHS separation is typically smaller than in experiments [7][8][9][10][62][63][64][65].
Besides Fermi level pinning, the enhancement of the DOS at the Fermi level due to the doping-induced flattening of the partially occupied bands is also relevant for understanding broken-symmetry phases, such as correlated insulator or superconducting states. In particular, the values of the transition temperatures to these states are usually very sensitive to the DOS at the Fermi energy, DOS(E F ). For example, the superconducting critical transition temperature is given by with V describing the coupling strength of the electrons to the superconducting glue (e.g., phonons or spin waves). The doping-induced increase of the DOS at the Fermi level should therefore result in a dramatic increase of the critical temperature. Again, this effect is not captured by tight-binding theory.
At θ = 1.2 • , very close to the magic angle, the VHS separation is only 5 meV in the undoped system and the valence VHS is much larger than the conduction VHS. Fermi level pinning is observed both for electron and hole doping. In the DOS, the shape of the VHS of the partially filled band is highly asymmetric. In particular, the leading edge of the peak (i.e., the side of the peak facing towards the other VHS) rises more sharply than the trailing edge (i.e., the side facing away from the other VHS). Interestingly, we also observe a double peak in the conduction VHS at ν = 1. The second peak is caused by a peak of the LDOS in the AB regions which does not coincide with the main peak of the LDOS in the AA regions. Again, this double peak structure is caused by the electron-electron interaction induced distortion of the conduction band near Γ. Fig. 1 shows that the conduction bands are extremely flat near Γ, but have a slightly higher energy than the states at M which give rise to the main peak of the VHS.

C. Environmental screening and comparison to the continuum model
So far, we have presented results for tBLG suspended in air ( bg = 1). In most experiments, however, the tBLG is placed on or sandwiched by a dielectric substrate (typically, hBN) and the presence of this dielectric environment screens the interaction between electrons in the tBLG [1-4, 8-10, 59]. In transport experiments, the dielectric substrate separates the tBLG from a metallic gate which is used to control the charge density in the tBLG and the presence of gates further modifies the  hBN). Surprisingly, the difference between the two band structures is small on the scale of the band width of the flat bands (similar band widths to those in experiments too). To understand this finding, we analyze the Hartree potentials of the two systems. Fig. 5(b) shows ∆V H (the difference between the Hartree potential in the centers of the AA and AB regions) as a function of doping for the two cases. While one might naively expect that the slope of ∆V H should be reduced by a factor of bg = 3.9 when the dielectric environment is included, we find that the reduction is much smaller  (∆V H is only reduced by 30% when the dielectric environment is included). The inclusion of metallic gates on both sides of hBN-encapsulated tBLG at a distance of 10 nm also has little effect on the band structure [ Fig. 5(a), purple dash-dotted line] because the Hartree potential does not change significantly, as shown in Fig. 5(b). It is worth noting that most experiments use larger gate distances than this, which would result in an even smaller effect. This surprising robustness of the Hartree band structure of tBLG towards changes in the dielectric environment has two reasons. First, the weakening of the Coulomb repulsion by the dielectric substrate allows for a greater inhomogeneity of the charge density. This results in a larger Hartree potential energy than the one that would have been obtained if the charge density had been frozen in its unscreened configuration. Second, the change in the dielectric environment only leads to small changes in the total screening response because the internal screening of the tBLG is already quite strong [32].
Except at the magic angle, the doping and twist-angle dependent atomistic Hartree potential energy is accurately described by where ν 0 (θ) is the doping level where the Hartree potential vanishes, V (θ) is a twist angle dependent energy parameter and G j denote the three reciprocal lattice vectors that are used to describe the out-of-plane corrugation of tBLG in Ref. [29]. Table I shows the optimal values of these parameters for the twist angles that we have studied and Fig. 6(b) compares the fit to the calculated Hartree potential as function of doping for different twist angles. Using Eq. (8) as an onsite energy in a tight-binding calculation allows the determination of Hartree-theory band structures without the need for selfconsistent calculations. We believe that this approach is a useful starting point for understanding broken symmetry phases in doped tBLG. IG. 6: (a) Hartree potential difference ∆V H between the AA and AB region as a function of twist angle for undoped (black) and electron-doped twisted bilayer graphene for ν = 1 (cyan), ν = 2 (blue) and ν = 3 (purple). (b) ∆V H as function of doping for three twist angles near the magic angle and linear fits obtained from Eq. (8).

IV. CONCLUSION
We have calculated quasiparticle properties, such as band structures and (local) densities of states, of interacting electrons in twisted bilayer graphene as function of doping and twist angle using atomistic Hartree theory. We find that doping results in significant changes to quasiparticle properties which are not captured by tightbinding approaches. In particular, we find that the partially occupied bands flatten between Γ and M in the Brillouin zone and even invert upon doping. The resulting local densities of states are in good agreement with recent scanning tunneling spectroscopy experiments. In particular, we capture the Fermi level pinning and the shapes of the van Hove singularities in the AA regions of tBLG. We predict that the band flattening gives rise to a strong enhancement of the peak in the AB regions. We also study the dependence of quasiparticle properties on the dielectric environment and find that they are not highly sensitive. This is a consequence of the strong internal screening of tBLG which suggests that the en-ergy scale of long-ranged Coulomb interactions is much smaller than previously estimated. As a consequence, the properties of broken symmetry phases of tBLG could result from a delicate interplay of long-ranged Coulomb interactions arising from the emergent moiré lattice and short-ranged atomic Hubbard interactions inherited from the untwisted bilayer. This will be the subject of future work.

V. ACKNOWLEDGEMENTS
We wish to thank K. Atalar, P. Guinea, N. Walet, D. Kennes and F. Corsetti for helpful discussions. We also wish to thank A. Kerelsky and L. Xian for sharing their data and for helpful discussions. ZG was supported through a studentship in the Centre for Doctoral Training on Theory and Simulation of Materials at Imperial College London funded by the EPSRC (EP/L015579/1). We acknowledge funding from EPSRC grant EP/S025324/1 and the Thomas Young Centre under grant number TYC-101. Hartree FIG. S1: Flat band structure along high symmetry path for various twist angles with the Hartree theory and tight-binding. All plots have been aligned at the K-point for clarity such that the energy, E, is relative to the K-point, and the energy scale of each plot is different to clearly show the flat band structures. Dotted-dashed horizontal lines denote the Fermi energy. The band distortions from the Hartree interaction are of 5-10 meV at charge neutrality, which can be significant right at the magic angle.