Layer-Dependent Interaction Effects in the Electronic Structure of Twisted Bilayer Graphene Devices

Near the magic angle, strong correlations drive many intriguing phases in twisted bilayer graphene (tBG) including unconventional superconductivity and chern insulation. Whether correlations can tune symmetry breaking phases in tBG at intermediate (≳ 2°) twist angles remains an open fundamental question. Here, using ARPES, we study the effects of many-body interactions and displacement field on the band structure of tBG devices at an intermediate (3°) twist angle. We observe a layer- and doping-dependent renormalization of bands at the K points that is qualitatively consistent with moiré models of the Hartree–Fock interaction. We provide evidence of correlation-enhanced inversion symmetry-breaking, manifested by gaps at the Dirac points that are tunable with doping. These results suggest that electronic interactions play a significant role in the physics of tBG even at intermediate twist angles and present a new pathway toward engineering band structure and symmetry-breaking phases in moiré heterostructures.


Sample Preparation:
Flakes of single-layer Graphene and hexagonal Boron Nitride were exfoliated onto Silicon Wafers with 90nm-thick oxide. The sample was constructed using a method similar to that used in [1]. A stamp comprised of Polypropylene carbonate (PPC), and Polydimethylsiloxane (PDMS), and transparent tape was used to pick up Graphite, hBN, and Graphene in sequential order. The PPC stamp holding the stack was flipped onto a 90nm oxidized Si wafer with the Graphene facing up, and the polymer was subsequently removed by annealing in a vacuum furnace at 350C for 10 hours. Contacts were patterned onto each sample surface using electron-beam lithography followed by evaporation of 5nm Cr and 50nm Au.

ARPES Measurements and Analysis:
Samples were measured using a Scienta R4000 Hemispherical Analyzer at the nanoARPES branch of beamline 7.0.2 (MAESTRO) at the Advanced Light Source using a photon energy of 74 eV, a temperature of 300K,and a pressure better than 1e-10 Torr. The beam was capillary refocused [2] to a spot size of ∼ 1µm x 1µm. The overall energy and momentum resolution was 30meV and 0.014Å −1 , respectively. The sample was doped electrostatically using a Keithley 2450 Source Meter.
All ARPES data in this paper were analyzed using pyARPES, an open-source pythonbased analysis framework [3]. Spectra presented in the figures have been deconvolved by the experimental energy and momentum resolution using the Lucy-Richardson method as described in ref [4]. Spectra presented in Figs 2-4 have additionally been divided by the Fermi-Dirac distribution following the methods in ref [4]. Second derivative spectra have been smoothed by 10 meV and 0.007Å −1 which is smaller than the energy and momentum resolution of the experiment. layers, calculated using Luttinger's theorem: n e = k F 2 /π, and normalized by the superlattice filling n s , whose formula is given in the text below. Error bars represent errors propagated from estimates of k F .

Carrier Density Measurements
The charge-carrier density can be obtained from the size of the Fermi surface using Luttinger's theorem. For each graphene layer, the Fermi surface (Fig 1a) is a circle with a radius of the Fermi wavevector k F , i.e. n e = k 2 F /π. We extract k F as half the distance between spectral peaks in the MDCs at the Fermi level, which are displayed in Fig 1b,c for the upper and lower layers of the 3 • twisted graphene, respectively. The summary of extracted k F values, and calculated n e values are presented in Fig 1d and Fig 1e, respectively.
Here we normalize the doping to the total filling of the moiré unit cell, which at small twist angles is approximated as n s = 4ν ≈ 10 4 8 * θ 2 a 2 √ 3 , where θ is the twist angle and a is the graphene lattice constant of 2.46Å [5]. Error bars for k F are estimated based on the broadness of the bands and the momentum resolution, which was ∼ 0.014Å −1 , while error bars for n e are obtained by propagating errors in k F .
As mentioned in the main text, the bottom layers of each sample consistently receives larger doping than does the upper layers due to the inability to screen the field produced by a back gate voltage [6,7]. Away from neutrality, both samples appear to have a roughly linear dependence between doping and gate voltage i.e. n e ≃ CV g , which is expected when treating the system as a parallel-plate capacitor with geometric capacitance C. At the neutrality point, a dip in geometric capacitance indicates the presence of a small gap, perhaps from the inversion symmetry broken by the hBN substrate [8] (see Supplementary Note 6 for more details).
Due to the relatively large twist angle of the samples measured here, the overall doping range is not very large with respect to the filling required to occupy 4 electrons per moiré unit cell. The presence of trace amounts of PPC, slightly hole dopes the sample, leaving the neutrality point around 1V applied back gate voltage [9].

Hartree Interaction Effects in 3 • tBG
The doping-induced renormalization effects presented in the main manuscript are an order of magnitude larger than what is predicted from Hartree and Hartree Fock interaction models in the literature for small twist (1.4 • ) angle graphene [10][11][12]. Here we present a Hartree model for 3 • twisted graphene in Fig2 which incorporates the band velocity enhancement from the long range electron-electron interaction in single layer graphene [9], and exhibits band renormalization of a qualitatively comparable to that of the experiment. We then  compare this model to a bare Hartree model which exhibits minimal band renormalization.
Below we breifly summarize the methodology behind these calculations.
The mean-field Hartree Hamiltonian we used is given in the following expression: where the first term is tight-binding Hamiltonian described in Supplementary Note 3, the second term is the on-site Hubbard Hamiltonian and the last is Hartree Hamiltonian with ϵ|r i −r j | in a unit of Hartree. The local density fluctuation δρ i,σ ≡ ⟨n i,σ ⟩ − n 0 i,σ is considered to avoid double-counting of coulomb interaction where n 0 i is local neutral density. In order to check the trend of doping-induced effects, the typical values of U = 5 and ϵ = 4 were used. The self-consistent calculation was performed on a 5 × 5 Monkhorst-Pack mesh.
To check the effect on the long-range exchange interaction, we defined an effective exchange model, which incorporates the long-range coulomb interaction effects present in single layer graphene on dielectric substrates near the charge neutrality point [9,[13][14][15]. In this model, the first order parameter t 0 of intra-layer hopping is redefined as a function of k point as ). suggesting they are indeed replicas of these original cones. Indeed, within each valley the distinct cones are separated by exactly 100meV, reflecting the energy separation between the upper and lower layers. As these replica cones are faint, they may be present but unreported in previous zone unfolded calculations [19].
Below we describe the methodology of the band structure calculations performed here, which are described extensively in ref [20]. The calculations use KLIFF [21]-fitted DRIP potential [22] parameters reproducing EXX-RPA level long-range interactions in bilayer graphene [23]. The energy mimization is performed using the LAMMPS package [24] using the cg algorithm with a timestep of 0.001 ps. The REBO2 Brenner potential [25] is used for We define the tight-binding (TB) Hamiltonian aŝ where |i⟩ is a basis of localized states at site i and the eigenfunctions are given as with n at the number of atoms and where k = (k x , k y ). The onsite energies ϵ i on the one hand are defined using the F2G2 model of graphene [26] and can be modulated due to an electric field by potential energy shifts given by where L 1 or L 2 indicates if an atom belongs to layer 1 or 2. On the other hand, the hopping terms t ij can be separated into an intralayer t intra ij contribution and an interlayer t inter ij contribution where the former follows the previously mentioned F2G2 model of graphene [26] and the latter are using a parametrizatioin that puts the twisted bilayer graphene magic angle at 1.08 • due to the specific choice of S = 0.895 when using the DRIP potential for interlayer interactions and the REBO2 potential for intralayer interactions during the LAMMPS minimization procedure in the following expression [27] t inter where p = 3.25Å and q = 1.34Å control the interlayer distance-dependent fitting of the tunneling at the K-point and and with the interlayer distance c 0 = 3.35Å, the rigid interatomic carbon distance in graphene a 0 = 1.42Å, the transfer integral between nearest-neighbor atoms V 0 ppπ = −2.7 eV, the transfer integral between two vertically aligned atoms V 0 ppσ = 0.48 eV, the decay length of the transfer integral set to r 0 = 0.184a such that the next-nearest intralayer coupling becomes 0.1 V 0 ppσ and the magnitude of the interatomic distance r ij = |r ij |. The cutoff for this distance-dependent model is set to 4.9Å beyond which additional contributions do not affect the observables anymore [28].
For the spectral function [29][30][31][32][33][34][35] calculations, we use the implementation outlined in Ref. [20] based on Ref. [29] for which we remind the expressions here. In such calculations, the zone-folded large supercells that capture the moire physics can be represented in the Brillouin zone of a smaller periodic unit cell through with |KJ⟩ eigenbands of the supercell that are labeled with capital letters and smaller letter n labels the Bloch function basis |kn⟩ with the localized orbital n in the reference small unit cell. The latter are used to to distinguish the layer and sublattice and in our case, we chose to represent the results for n = 1. In the above expression A KJ,KJ (E) reduces to a δ(E − ϵ KJ ) function at the eigenvalue of the superlattice system and ⟨kn|KJ⟩ is a structure factor that is modulated by a position-dependent phase term as (BZ) respectively. w N ≤ 1 is a coefficient that allows to tune the relative contribution of certain atoms to capture the top layer contribution that is usually stronger in experiment.
Here we set this value to 1 as we plot the top and bottom layer contributions separately.
Photon polarization effects are ignored although they could alter the momentum distribution anisotropy [36,37].

Gap Analysis at the K Point
The analysis here is very sensitive to the exact locations of the K upper and K lower points.
Indeed, given the steep linear dispersion of graphene, a small misalignment can result in an overestimate of the band gap size [8,38,39]. Figs 4 and 5 address the gap location with high momentum resolution ARPES spectra surrounding K upper and K lower points, respectively.

Twist Angle Inhomogeneity Effects on Electronic Structure
We can estimate the twist angle homogeneity of the sample measured in this report using spatially resolved ARPES. Clearly, the sample has a spatially inhomogeneous twist angle. Since distinct electronic structure can be seen for the two sample regions, the electronic structure appears to be locally homogeneous over the length scale of the 1 µm beam spot, similar to findings in the literature [40].
Twist angle inhomogeneity often corresponds to doping inhomogeneity, both of which are present in tear-and-stack tw-BLG samples in the literature [10,40,41] as well as our sample measured in this report. The careful reader will notice that the valence band Dirac point is This is confirmed by the carrier density as a function of doping in panel g, extracted using Luttinger's Theoreom, i.e. n e = k 2 F /π, which has a minimum at 1.5V for the 3 • region, and 0.5V for the 2.3 • region. Such doping inhomogeneity present in our sample may provide an additional explanation for why macroscopic measurements such as transport do not find very large gaps at the neutrality point [42,43]: since different regions have different dopings, one sample region that is in the gap and insulating will be adjacent to another sample region that has carriers populated and thus is conducting.
To our knowledge the local gap size is not strongly affected by the twist angle within this particular twist angle range measured in the experiment . Fig 6h and i present equilibrium ARPES second derivative spectra along Γ − K upper for the 3 • and 2.3 • regions, respectively.
Both twist angle regions exhibit ≈ 100 meV separation between conduction and valence bands, which is confirmed by the region in energy where there is deviation from linearity in the dispersions (right sides of panels h and i). Indeed, the strain that produces the twist angle inhomogeneity on our sample clearly does not have such a strong affect on the gap size at the neutralty point in our experimental conditions, as suggested by early theoretical works on tw-BLG [44]. This is clearly not the case for the magic angle twist regime at low temperature, where there is a competition between Kekule spiral order and K-IVC states [45,46].

Spatial Inhomogeneity Effects on Carrier Density Measurements
The spatial inhomogeneity ubiquitous in materials often precludes bulk probes such as transport from producing accurate readings of sample carrier concentration. Typically, this shows up when the surface carrier density is different from that of the bulk, e.g. from surface dosing [47,48], or the presence of surface and interface states [49,50] present on bulk conductive materials. A local and surface sensitive probe such as nanoARPES is therefore much better suited to access the carrier density because, as mentioned in the previous note, it can locally measure Luttinger volume of the Fermi surface [47,49,50].
The spatial inhomogeneity of our sample comes into account when estimating the carrier density using a simple capacitance model, i.e. n = C(n)V g where the capacitance C(n) is determined by 1 C = 1 Cg + 1 Cq , and C g and C q are the geometric and quantum contributions to the capacitance, respectively [51]. Let's focus first on the geometric capacitance, which dominates away from charge neutrality when our sample is more or less metallic: the hBN contribution is given by C g,hBN = ϵ hBN ϵ 0 /d hBN ≈ 0.19 µF cm 2 for an hBN thickness of 15nm and dielectric strength ϵ hBN = 3.4. Interestingly, this value vastly overestimates the carrier concentration obtained via ARPES (see Fig 7d away from neutrality).
We attribute this discrepancy to the presence of an air gap between the sample and the graphite back gate, perhaps due to the proximity of measurement locations to a wrinkle on hBN (pink dashed line in Fig 7a). Indeed, a recent study of the hBN dielectric strength found that such air gaps develop in about a third of devices, substantially decreasing the overall effective dielectric strength of the capacitor [52]. The air gap behaves as a capacitor in series, i.e. 1 Cg = d ϵr = 1 C g,hBN + 1 C g,air . By plugging in dielectric strengths of 3.4 [52] for hBN and 1 for air, and a thickness of d hBN = 15nm (see panel b), we estimate that an 8nm air gap is present between the hBN and our sample, producing an overall effective dielectric strength of ϵ r = 2, resulting in doping behavior presented in Fig 7c. Overall, the spatial inhomogeneity present in our twisted graphene/hBN heterostructure in the form of an air gap results in an inaccurate estimate of the carrier density, which is resolved using nanoARPES. Finally, we comment briefly on the quantum capacitance contributions, which occur at low carrier density. Quantum capacitance is general given by the density of states i.e.

Supplementary
C q = e 2 D(E), and in a massless Dirac fermion model of twisted bilayer graphene this is 2 · 10 6 m/s is the Fermi velocity, and n d = 1.0 · 10 1 0cm −2 is the excess carrier density present due to disorder [53]. This quantum capacitance  The data reported in the main manuscript provide evidence for a method of tuning the band velocities and band gaps in twisted bilayer graphene in operando. As we will argue below, we believe that these effects can be best explained by a combination of the substrate interaction and the spatially inhomogeneous Hartree-Fock interaction, which is controlled within different layers using a displacement field.
First we address the K point dispersions, which can be modified by several different mechanisms, including single particle effects from band displacement and strain, as well as many-body effects such as electron-electron interactions.
On a single particle level, band displacement in twisted bilayer graphene brings the Dirac point of one layer closer in energy to the valence van-Hove-singularity E vHs , and the opposite for the other layer [7,54], while keeping the overall bandwidth constant. At large displacement (D ≈ E vHs ), the valence band velocity at the Dirac point which is energetically closer to the van Hove singularity decreases, and the band velocity at the Dirac point which is farther from the van Hove singularity increases [54]. However, this behavior is qualitatively inconsistent with our results, as the band velocity at both Dirac points is modified in the same direction upon changing the doping.
Both uniaxial and heterostrain can also change the bandwidth in tw-BLG [43,44]. Uniaxial strain brings the Dirac points of different layers in closer proximity in momentum space [55], producing an effect similar to that produced by decreasing the overall twist angle [56,57]. Heterostrain, on the other hand adjusts the interlayer hopping such that the first magic angle condition for flat bands occurs at larger twist angle, and notably resulting in a shift of the Dirac cones from the mini Brillouin zone corners [44]. However, typical strain present in twisted graphene samples is ≈ 1% [58], which is an order of magnitude smaller than the 16% strain required to observe the band structure modifications of ≈ 40% in our experiment [55]. Additionally, the Dirac cone momentum separation does not change with application of a gate voltage in our experiment, suggesting that any strain present in our sample is not gate-tunable, and therefore cannot explain the doping behavior of the band velocity.
As mentioned earlier, electron-electron interactions are known causes of significant band velocity modification both in single layer [9,13,14,[59][60][61] and twisted bilayer [10,[16][17][18]62] graphene. In the case of single layer graphene, the band velocity enhancement is strongest at the the charge neutrality point, and weaker for electron and hole-dopings [9,59,61], which is inconsistent with our data. On the other hand, the spatially inhomogeneous Hartree Fock interaction in twisted blayer graphene qualitatively matches the doping dependence seen in our data at the K points: i.e. a steady valence band velocity decrease with holedoping [10,11]. We therefore attribute the doping behavior of the K point dispersions to such Hartree-Fock interactions.
We now address the gap present in our data. Several mechanisms can produce a band separation in graphene-based systems, including many body interactions, interlayer potentials, and inversion symmetry breaking. The electron-plasmon interaction, for example, has been demonstrated to generate a plasmaron satellite in highly doped graphene systems [63][64][65][66], which can take the appearance of a gap. However, the doping dependence of the satellite separation for a plasmaron is linear in E F , and therefore it cannot explain the presence of the gap at the neutrality point at zero doping. Additionally, the plasmaron has been predicted [67] and demonstrated [66,68] to compete with gaps, as in the case of graphene aligned to hBN [66].
Interlayer potentials may also produce band gaps in twisted bilayer graphene. At particular commensurate twist angles, the K points from each layer become equivalent points in the mini Brillouin zone , resulting in an umklapp scattering gap whose magnitude scales positively with the twist angle [69][70][71][72][73][74]. However, our sample has a twist angle of 3 • , rendering it among the twist angles where K upper and K lower are not connected by a reciprocal lattice vector and such umklapp scattering is not predicted [71,72].
Inversion symmetry breaking, is well known to produce band gaps in graphene [8,75,76] and commonly occurs in the presence of a substrate with a different lattice constant [8,[76][77][78][79], even when the substrate is crystallographically misaligned [76,78,79] as in the case of our sample. Though a gap from a misaligned hBN substrate may be small [78,80], interactions may enhance them significantly [80,81]. For example, electron-electron interactions enhance band gaps in graphene with a magnitude that scales with the interaction strength [80][81][82], which in tBG scales linearly with doping [16][17][18]62]. Upon generating a charge imbalance between the two graphene layers, e.g. via a displacement field, this mechanism enables layer-dependent gap enhancements such as those seen in our experiment.
These observations make the substrate and Hartree-Fock interactions the most likely candidates for generating and enhancing the band gap in 3 • twisted bilayer graphene. Indeed these interactions can qualitatively explain both the presence of the gap at charge neutrality as well as the layer-dependent gap enhancement upon electron doping the sample.