Atomistic Hartree theory of twisted double bilayer graphene near the magic angle

Twisted double bilayer graphene (tDBLG) is a moiré material that has recently generated significant interest because of the observation of correlated phases near the magic angle. We carry out atomistic Hartree theory calculations to study the role of electron–electron interactions in the normal state of tDBLG. In contrast to twisted bilayer graphene, we find that such interactions do not result in significant doping-dependent deformations of the electronic band structure of tDBLG. However, interactions play an important role for the electronic structure in the presence of a perpendicular electric field as they screen the external field. Finally, we analyze the contribution of the Hartree potential to the crystal field, i.e. the on-site energy difference between the inner and outer layers. We find that the on-site energy obtained from Hartree theory has the same sign, but a smaller magnitude compared to previous studies in which the on-site energy was determined by fitting tight-binding results to ab initio density-functional theory (DFT) band structures. To understand this quantitative difference, we analyze the ab initio Kohn–Sham potential obtained from DFT and find that a subtle interplay of electron–electron and electron–ion interactions determines the magnitude of the on-site potential.

Twisted double bilayer graphene (tDBLG) is a moiré material that has recently generated significant interest because of the observation of correlated phases near the magic angle.We carry out atomistic Hartree theory calculations to study the role of electron-electron interactions in the normal state of tDBLG.In contrast to twisted bilayer graphene (tBLG), we find that such interactions do not result in significant doping-dependent deformations of the electronic band structure of tDBLG.However, interactions play an important role for the electronic structure in the presence of a perpendicular electric field as they screen the external field.Finally, we analyze the contribution of the Hartree potential to the crystal field, i.e. the on-site energy difference between the inner and outer layers.We find that the on-site energy obtained from Hartree theory has the same sign, but a smaller magnitude compared to previous studies in which the on-site energy was determined by fitting tight-binding results to ab initio density-functional theory (DFT) band structures.To understand this quantitative difference, we analyze the ab initio Kohn-Sham potential obtained from DFT and find that a subtle interplay of electron-electron and electron-ion interactions determines the magnitude of the on-site potential.
In transport experiments on magic-angle tDBLG, correlated insulators have been observed at a doping of two electrons per moiré unit cell in applied electric fields [23,25,[28][29][30].For the undoped system, a substantial band gap emerges in the presence of a perpendicular electric field; this band gap is not caused by electron interactions [22].Signatures of superconductivity have also been reported, but robust superconductivity has not yet been confirmed in this system [23,25,[28][29][30][31].Scanning tunneling microscopy (STM) experiments have also observed correlated insulating states [32,33] and nematic ordering [34] near the magic angle.
The electronic structure of tDBLG has been studied using a variety of methods including continuum theories [22] and atomistic tight-binding models [24].In Ref. 24, it was found that the band gap of the undoped system obtained from tight binding differs from the result of first-principles density functional theory (DFT).To obtain better agreement, a phenomenological on-site energy of approximately −30 meV was added to the in-ner layers in the tight-binding calculations.This on-site potential has been referred to as the "intrinsic symmetric polarisation" (ISP) or crystal field.Similar results have been found in Refs.35 and 36.These works suggested that the difference in chemical environments between the inner and outer layers of tDBLG can result in charge transfer between layers which in turn gives rise to significant electron-electron interaction effects in the normal state.
In contrast to tBLG, no Fermi level pinning has been observed in tDBLG, in agreement with continuum model calculations of the electronic structure [32][33][34].However, such continuum model calculations only capture the longrange interaction between electrons, but fail to accurately describe the interaction between nearby electrons on different layers.To address this issue, it would be highly desirable to carry out atomistic Hartree theory calculations of tDBLG.
In this paper, we investigate the role of electronelectron interactions in tDBLG within atomistic Hartree theory.We find that the electronic band structure of tD-BLG is not sensitive to electron or hole doping -in stark contrast to tBLG.Next, we consider the effect of a perpendicular electric field and find that electron-electron interactions play an important role in screening the externally applied field.Finally, we evaluate the contribution of the Hartree potential to the crystal field and find that it has the correct sign, but is too small compared to previous findings.To understand this discrepancy, we analyze the Kohn-Sham potential of an ab initio DFT calculation of tDLBG which produces a crystal field in good agreement with previous results.This suggests that ab initio calculations are required for a quantitatively accurate description of the crystal field in tDBLG.

II. METHODS
We study commensurate moiré unit cells of tDBLG consisting of twisted AB stacked bilayers.The bilayers are initially stacked directly on top of each other, similar to the structure of graphite, and the top bilayer is rotated anticlockwise about an axis normal to the bilayers that passes through a carbon atom in each bilayer.The moiré lattice vectors are R 1 = na 1 + ma 2 and R 2 = −ma 1 + (n+m)a 2 [6], where n and m are integers that specify the moiré unit cell in terms of the graphene lattice vectors a 1 and a 2 .
We relaxed these tDBLG structures using classical force fields as implemented in the LAMMPS software package [43].The AIREBO-Morse potential [44] was used for intralayer interactions; while the Kolmogorov-Crespi potential [45] was used for interlayer interactions.More details can be found in Ref. 46.
The electronic structure of tDBLG was investigated with an atomistic Hartree theory.The model outlined here is very closely related to that discussed in Ref. 40.For completeness, we provide all details here.The atomistic Hamiltonian that we solve is given by where ĉ † i and ĉi are, respectively, the electron creation and annihilation operators associated with the p z -orbital on atom i, and ε i is its on-site energy.
The hopping parameters t(τ i − τ j ) between atoms i and j (located at τ i/j ) are determined using the Slater-Koster rules [47,48] where γ 1 = 0.48 eV and γ 0 = 2.81 eV [24] correspond, respectively, to σ-and π-hopping between p zorbitals, with associated decay parameters q σ = 7.43 and q π = 3.14 [6,7].Also, a = 1.397Å is the pristine carbon-carbon bond length, d = 3.35 Å is the pristine interlayer separation parameter, and ϕ is the angle between the z-axis and the vector connecting atoms i and j, and captures the angle-dependence of hoppings.Hoppings between carbon atoms that are separated by more than 10 Å are neglected [49].The Slater-Koster tight-binding parameters are based on a best fit to DFT low-energy bandstructures of graphene and bilayer graphene [6,7,48], with a slightly larger π-hopping parameter as introduced in Ref. 24 to improve agreement with the low-energy DFT bandstructure of tDBLG at large twist angles.
We decompose the on-site energy ε i in Eq. ( 1) into two contributions, ε i = ε αi + ε el i .The first term, ε αi , is constant within each layer (with α i denoting the layer in which atom i resides) and represents the effect of an applied electric field.The second term, ε el i , is the contribution to the on-site energy from electron interactions and is determined self-consistently according to where φ z (r) is the p z orbital of the carbon atoms, and the Hartree potential V H (r) is determined from the electron density n(r) and the screened electron-electron interaction W (r) via Here, n 0 (r) is a reference electron density that ensures overall charge neutrality and whose subtraction in Eq. ( 4) avoids double counting of Hartree interactions in the uniform system [7].The electron density is determined through where denotes the Bloch eigenstates of Eq. ( 1), with subscripts n and k denoting the band index and the crystal momentum, respectively.Also, N k is the number of kpoints in the summation of the electron density, and is the spin-degenerate occupancy of state ψ nk with eigenvalue ε nk (where ε F is the Fermi energy).Inserting the Bloch states into Eq.(5) gives where χ j (r) = R φ 2 z (r − τ j − R) (with R denoting the moiré lattice vectors) and the total number of electrons on the j-th p z -orbital in the unit cell being determined by The reference density is taken to be that of a uniform system, n 0 (r) = n j χ j (r), where n is the average of n j over all atoms in the unit cell, which is related to the filling per moiré unit cell ν through n = 1 + ν/N , where N is the total number of atoms in a moiré unit cell [39].
In transport experiments, there is often a metallic gate above and below the tDBLG, with a hexagonal boron nitride (hBN) substrate separating the gates from tDBLG.These metallic gates add or remove electrons from tD-BLG and can also create electric fields across the system.These gates also screen the electron interactions in tD-BLG, and taking this effect into account has been shown Right panels: results for tBLG, at the same twist angles, generated using the method outlined in Ref. 42.
to be important in tBLG [50][51][52].Therefore, we use a double metallic gate screened interaction where ξ is the thickness of the hBN dielectric substrate, which provides a background dielectric constant bg and separates tDBLG from the metallic gate on each side [50,53,54].We set ξ = 10 nm in all calculations, and bg = 4, which corresponds to experiments in which tDBLG is encapsulated in hBN [55], unless otherwise stated.Note that we do not consider the case of the hBN being closely aligned with the graphene layers which induces significant changes to the electronic structure of the graphene layers [56].
In our atomistic model, we neglect contributions to the electron density from overlapping p z -orbitals that do not belong to the same carbon atom, which is equivalent to treating φ 2 z (r) as a delta-function.Therefore, we calculate the Hartree on-site energies using with W Rij = W (R + τ j − τ i ).If R = 0 and i = j, we set W 0,ii = U/ bg with U = 17 eV [57].
To obtain a self-consistent solution of the Hartree theory, we use a 6×6 k-point grid to sample the first Brillouin zone to converge the density in Eq. ( 5) and we sum over a 21×21 supercell of moiré unit cells to converge the on-site energies of Eq. ( 9).Linear mixing of the electron density is performed with a mixing parameter of 0.1 or less (i.e., 10 percent of the new potential is added to 90 percent of the potential from the previous iteration).Typically, the Hartree potential converges to an accuracy of better than 0.1 meV per atom within 100 iterations.For doping levels where tDBLG is metallic, smaller mixing values and a larger number of iterations are sometimes needed to reach this convergence threshold.
We use the Perdew-Burke-Ernzerhof (PBE) exchange-correlation functional [60].For the electron-ion interaction we employ the projector-augmented-wave (PAW) formalism [61,62] (where PAW pseudopotentials are generated from ultra-soft pseudopotentials [63]).The kinetic energy cutoff is set to 800 eV.A basis consisting of four non-orthogonal generalized Wannier functions (NGWFs) per carbon atom is used, which are optimised in situ, and the ensemble-DFT approach [64,65] is adopted for the minimisation of the ground state energy.

A. Doping dependent band structure
In tDBLG, the low-energy electronic band structure is characterized by a set of four bands that are separated from all other bands.These bands become extremely flat  , the low-energy flat bands near the K and K points exhibit a parabolic dispersion which is inherited from the parent AB-stacked bilayers.However, the two graphene sheets of the AB-stacked bilayers in tDBLG are no longer equivalent and this results in the opening of a band gap at the K and K points -similar to the case when an electric field is applied perpendicular to an AB-stacked bilayer.
In Fig. 1 we show the Hartree theory band structure of doped tDBLG at a twist angle of θ = 1.89 • (top left panel) and θ = 1.41 • (bottom left panel), and compare the results to tBLG (right panels).Results are shown for −3 ≤ ν ≤ 3, where ν represents the number of electrons (ν > 0) or holes (ν < 0) added to each moiré unit cell.In stark contrast with tBLG [37][38][39][40][41][42][66][67][68], it can be seen that the band structure of tDBLG does not significantly change upon doping (all doping levels have been aligned such that the zero energy occurs at the mid-point between the upper and lower flat bands at the K-point).The flat bands distort only by a few meV relative to the charge neutral case (black lines), which is small compared to their bandwidths (129 meV at 1.89 • and 13 meV at 1.41 • ).By comparison, in tBLG the Hartree interactions induce band deformations of ∼25 meV.At the twist angle of 1.89 • (with a bandwidth of 260 meV), these deformations are modest in comparison to the bandwidth, but at the angle of 1.41 • (with a bandwidth of 111 meV), they are starting to become comparable to the bandwidth.Very close to the magic angle of tBLG (approximately 1.1 • ), the bandwidth shrinks to ∼5 meV and the Hartree deformations become the dominant energy scale [37][38][39][40][41][42][66][67][68].
To understand the absence of significant band deformations in tDBLG, we analyze the Hartree potential.Fig. 2 shows the locally-averaged Hartree potential along the diagonal of the moiré unit cell on one of the outer layers (left panels) and one of the inner layers (right panels) for different doping levels ν.We find that the Hartree potential varies by ∼25 meV on the inner layers, but only by ∼10 meV on the outer layers.In contrast, the Hartree potential of tBLG (and also of twisted trilayer graphene [69]) varies by ∼100 meV and is therefore the dominant energy scale for a wide range of twist angles near the magic angle [37].
Upon electron doping tDBLG (ν > 0, top panels of Fig. 2), a positive peak in the Hartree potential emerges in the AA regions.When electrons are removed (ν < 0, bottom panels of Fig. 2) a negative trough is found instead in the AA regions.This is a consequence of the shape of the flat-band wavefunctions near the K and K points.Specifically, these states are localized in the AA regions of the inner layers, similar to tBLG [37][38][39][40][41][42] and tTLG [69].
However, an important difference (besides the strength of the Hartree potential) between tDBLG and tBLG is the degree of real-space localization of flat-band states in different parts of the first Brillouin zone.Fig. 3 shows that in tBLG that flat-band states at K and M are strongly localized in the AA regions, while the states at Γ form rings around the AA regions.As a consequence, the Hartree potential gives rise to large energy shifts of the states at K and M relative to the states at Γ (see Refs. [37][38][39][40][41][42]69 for a more detailed description of these band deformations).
In contrast, Figs. 4 and 5 show that flat-band states in tDBLG are much less localized.For example, the states at Γ form rings around the AA regions in the inner layers (see Fig. 4), but are localized in the AA regions on the outer layers (see Fig. 5).The states at M are localized in the AA regions in the inner layers, but have delocalized stripe-like features on the outer layers.As a consequence, the Hartree potential does not give rise to significant relative shifts of these states.Finally, the valence states at the K-point are localised on the AA regions of the inner layers, but on the AB/BA regions of the outer layers; while the conduction states at K are almost entirely delocalized on the outer layers [36].This explains why significant band deformations are not observed in tDBLG even though the Hartree potential has a similar magnitude as the band width close to the magic angle.A similar explanation for the absence of strong band deformations was also offered in the continuum model of Ref. 70.In the case of tBLG, the interaction-induced dopingdependent band distortions were shown to cause a pinning of the Fermi level at the van Hove singularities [37][38][39][40][41][42], which was observed in tunneling experiments [14][15][16]71].In tDBLG, as no significant band distortions are observed, we do not expect a similar Fermi level pinning as in tBLG.Indeed, recent tunnelling experiments did not observe Fermi level pinning at the van Hove singularity [32][33][34].Note, however, that in these experiments some changes of the electronic structure were observed as function of doping.It was proposed [34] that these changes are a consequence of the perpendicular electric field which accompanies doping in a device with a single metallic gate.We therefore study the effect of such electric fields in the next section.

B. Electric fields
Figure 6 compares the Hartree (red) and tight-binding (black) band structures of 1.89 • tDBLG at charge neutrality for two different electric field strengths.Application of a perpendicular electric field increases the gap between the valence and conduction bands.Without the field, the value of the gap is 9.9 meV, while its value is 19.2 meV for a field strength 30 meV Å−1 .Moreover, the electric field lifts the valley degeneracy of both the valence and conduction bands (except at the Γ and M points).This splitting of the conduction and valence bands increases substantially with the strength of the applied field with the valence bands undergoing more significant distortions than the conduction bands.
In the tight-binding approximation the band distortions are significantly more pronounced than in Hartree theory.The electric field causes the system to polarise such that in one of the bilayers there is an enrichment of electrons and in the other bilayer there is a depletion of electrons, with the outer layers exhibiting larger en- richment/depletion than the inner layers.When Hartree interactions are included, the Hartree potential opposes the external electric field to reduce its effect.The extent to which the Hartree potential screens the electric field can be determined by computing an effective dielectric constant for each layer is the average potential due to the electric field in layer α, with E denoting the electric field strength, z (α) i the z-coordinate of atom i in layer α, and • • • an average over i in layer α.Also, is the averaged Hartree potential in each layer α.For E = 10 meV Å−1 , we find (α) = 2.60 for all layers, while for E = 30 meV Å−1 the effective dielectric constant is reduced to 1.94.This reflects the fact that the electrons cannot screen larger electric fields as effectively.Note that these values were obtained with bg = 4.For freestanding tDBLG, ab initio DFT calculations have found an effective perpendicular dielectric constant of approximately 3 [72,73] for untwisted graphene multilayers in vacuum.These ab initio values are qualitatively similar to our findings.

C. Crystal field
Previous work established that an additional layerdependent on-site potential (often referred to as the crystal field) must be included in tight-binding calculations of tBLG to achieve agreement with ab initio DFT band structures [24,35,36].However, the origin of this on-site potential has remained unclear.In this section we evaluate the effective on-site energy generated by the Hartree potential and compare to previous work.
At charge neutrality, the Hartree potential is more negative on the inner layers than on the outer layers, see Fig. 2, indicating that electrons have transferred from the inner to the outer layers.Adapting the definition of Ref. 35 for the crystal field, we define δ as the difference between the layer-averaged on-site Hartree potential on the outer and inner layers according to where ε el out and ε el in are the on-site energies averaged over atoms in the outer and inner layer of a bilayer, respectively (note that the two bilayers of tDBLG are equivalent by symmetry).
The left panel of Fig. 7 shows the values of the Hartree crystal field δ as function of twist angle for different values of bg .We find that δ does not sensitively depend on the twist angle, but increases as the background dielectric constant is reduced.For freestanding tDBLG (corresponding to bg = 1), we find δ ≈ 8 meV near the magic angle.We have verified that the band structure from a tight-binding calculation with a layer-dependent on-site potential of 8 meV agrees well with the full Hartree theory result indicating that the in-plane variations of the Hartree potential do not play an important role.Note that the value of δ has the same sign, but is somewhat smaller than the value proposed by Haddidi and coworkers [24] who used a value of 30 meV.
To understand the weak twist-angle dependence of δ, we analyze the charge transfer between inner and outer layers.Within an idealized parallel plate capacitor model, δ should be proportional to the charge density per unit area of each graphene layer.The total number of polarized charges per moiré cell in layer α is given by ∆n where j runs through all the atoms in layer α.By symmetry ∆n has the same magnitude but opposite sign for the outer and inner layer of each bilayer.In the right panel of Fig. 7 we show the dependence of ∆n (1/4) (i.e., the charge on the outer layers) on the area of the moiré  7: Left panel: The average Hartree potential difference between the outer and inner layers, δ, for several dielectric constants as a function of twist angle θ.Right panel: The total excess number of electrons ∆n (1/4) on one of the outer layers (which we refer to as layers 1 and 4) within the moiré unit cell at various dielectric constants bg as a function of moiré unit cell area A. unit cell.It can be seen that ∆n (1/4) behaves approximately linearly resulting in a constant charge density per unit area which in turn gives rise to a constant δ.
As the Hartree theory contribution to the crystal field from our atomistic model is significantly smaller than the value determined by Haddidi and coworkers [24], we also analyze the full Kohn-Sham potential obtained from an ab initio DFT calculation of tDBLG with a twist angle of 2.45 • .The Kohn-Sham potential has three contributions: (1) the ion-electron (ion-el) potential which is often approximated with a pseudo-potential; (2) the Hartree contribution from electron-electron (el-el) interactions; and (3) the exchange-correlation contribution.
Figure 8 shows the averaged Kohn-Sham potential of tDBLG as function of z.The potential has been averaged over the x and y directions, and the resulting function of z has been smoothed by taking its convolution with a rectangular function of width 3.20 Å (we have verified that our results do not depend sensitively on the value of this width).In the left panel, we subtract the exchange-correlation contribution from the Kohn-Sham potential which yields an ab initio Hartree theory potential (note, however, that the charge density was obtained from the full potential including exchange-correlation effects).We find that the ab initio Hartree potential is approximately 30 meV smaller on the inner layers than the outer layers, as indicated with the horizontal lines.When the exchange-correlation potential is included (right panel of Fig. 7), the potential difference between inner and outer layers increases to approximately 40 meV, in good agreement with the findings in Refs.24, 35 and 36.These results demonstrate that an ab initio description of the potential is required to obtain a quantitatively accurate description of the crystal field.

IV. CONCLUSION
We have used atomistic Hartree theory calculation to investigate the role of electron-electron interactions in tDBLG and studied the effects of changes in twist angle, dielectric environment, doping and applied electric fields.Our calculations reveal that the band structure of tDBLG is largely insensitive to electron or hole doping in stark contrast to tBLG.Application of a perpendicular electric field changes the band gap and lifts the valley degeneracy of the bands.Electron-electron interactions screen the electric field and we obtain an effective dielectric constant that is quantitatively similar to ab initio results for untwisted graphene multilayers.Finally, we analyze the contribution of Hartree interactions to the crystal field which is defined as the difference between the on-site energies of the inner and outer layers.We find that the difference of the average Hartree potentials in the inner and outer layers has the same sign, but a smaller magnitude compared to previous studies which determined the on-site potential by fitting tight-binding band structures to ab initio DFT results.To understand this difference, we carry out ab initio DFT calculations of tDBLG and analyze the Kohn-Sham potential.We find that the locally averaged Kohn-Sham potential difference between the inner and outer layers agrees well with the previously reported value of the crystal field, indicating that a quantitative description of this effect requires an ab initio description of the subtle interplay between electron-ion and electron-electron interactions.

FIG. 1 :
FIG.1: Left panels: Hartree theory band structure of doped (−3 ≤ ν ≤ 3) tDBLG at a twist angle of θ =1.89 • (top) and θ =1.41 • (bottom) with bg = 4.The horizontal lines denote the Fermi energy at each doping level.The band structure at charge neutrality is shown in black and the band structures at all doping levels are practically identical to it.Note all bands are aligned such that the zero of energy occurs in the middle of the band gap at the K-point.Right panels: results for tBLG, at the same twist angles, generated using the method outlined in Ref.42.

FIG. 2 : 2 FIG. 3 :
FIG. 2: Locally-averaged Hartree potential along the diagonal of the moiré unit cell for electron and hole doped tDBLG with θ =1.89 • and bg = 4.The variable |s| measures the distance from the BA site of the moiré unit cell to a point s along its long diagonal.The vertical solid lines correspond to BA stacking of the inner layers, dotted-dashed lines to the AB stacking, and dotted lines to AA stacking.The left panels show results for the outer layers and the right panels are for the inner layers.Results for electron doped systems (ν > 0) are shown in the upper panels and results for hole doped systems (ν < 0) are shown in the lower panels.The Hartree potential in each layer has been locally averaged: the value at each atomic position was obtained by averaging the Hartree potential at this site with the Hartree potentials of the three nearest neighbours.This removes the atomic-scale oscillations of the Hartree potential.

2 FIG. 4 : 2 FIG. 5 :
FIG. 4: Square modulus of the flat-band wavefunctions on an inner layer of tDBLG at different crystal momenta.Top panels show results for the flat conduction band (whose square moduli have been summed), and the bottom panels show similar results for the two flat valence band states.See the caption of Fig. 3 for further details.

FIG. 8 :
FIG. 8: Left panel: Sum of Hartree potential and ion-electron potential as a function of z (the coordinate perpendicular to the plane of the tDBLG system) obtained from an ab initio DFT calculation of tDBLG at a twist angle of 2.45 • .Right panel: Full Kohn-Sham potential (including the exchange-correlation contribution) as a function of z in 2.45 • tDBLG.The ab initio potentials are first averaged over the x and y coordinates, and the resulting function of z is smoothed by taking its convolution with a rectangular function of width 3.20 Å.The dotted vertical lines correspond to the z-averaged atomic positions of each layer.The horizontal dotted-dashed lines indicate where the potential crosses the z-averaged atomic positions of each layer.