Electrostatic effects, band distortions, and superconductivity in twisted graphene bilayers.

Bilayer graphene twisted by a small angle shows a significant charge modulation away from neutrality, as the charge in the narrow bands near the Dirac point is mostly localized in a fraction of the Moiré unit cell. The resulting electrostatic potential leads to a filling-dependent change in the low-energy bands, of a magnitude comparable to or larger than the bandwidth. These modifications can be expressed in terms of new electron-electron interactions, which, when expressed in a local basis, describe electron-assisted hopping terms. These interactions favor superconductivity at certain fillings.

The discovery of strong repulsive interactions and superconductivity in twisted graphene bilayers (1–3) [see also previous experimental results (4, 5) and recent work (6, 7)] has led to an increased interest in the role of the electron–electron interaction in these systems. Previously, extensive experimental studies had been carried out in graphene superlattices on a boron nitride (BN) substrate; for example, refs. 8–11. The mismatch in lattice spacing between graphene and BN limits the maximum wavelength of the Moiré lattice. As the lattice constant in the two graphene layers in a twisted bilayer is the same, the periodicity of the Moiré structure diverges at small twist angles (12–16). For sufficiently small angles, almost flat bands arise near the charge-neutrality point (14, 17). The effects of the intrinsically small interaction effects in graphene are expected to be enhanced for those “magic” angles where the width of the low-energy bands is smallest. Novel magnetic phases become possible when the lowest band is half filled (18). Layer-dependent strains can also lead to Moiré structures and narrow bands (19, 20).

In the following, we consider the effect of the long-range Coulomb interaction on the bands nearest to the neutrality point, as a function of their filling. The occupation of the low-energy bands in a twisted bilayer graphene leads to inhomogeneous electrostatic potentials of order , where is the length of the Moiré lattice unit, and is the dielectric constant of the environment. We expect this interaction to be comparable, or larger, than the bandwidth, and thus it will give rise to a significant distortion of the bands. We analyze the screening of the Coulomb interaction within the Hartree–Fock approximation. As this is a variational method, it can be reliably used to treat screening effects, even when the interaction energy is comparable to the kinetic and crystal field contributions.

We find that the electrostatic potential due to the inhomogeneous charge distribution in the Moiré unit cell gives rise to significant distortions in the band. Using general properties of the Wannier functions which describe this band (21–23), we identify local, filling-dependent couplings which change the band dispersion away from the neutrality point. Finally, we discuss the relation between these interactions, which can be defined as electron-assisted hopping terms, and superconductivity. Our work presents a rather different angle and thus complements the extensive corpus of work on the role of local interactions (24–33). We will not consider the role of lattice relaxation (34), since this is a complex problem and will be discussed separately. We also do not address the role of the electron–phonon interaction on superconductivity; see refs. 35 and 36.

Estimates of the Electron–Electron Interactions

We evaluate electron–electron interactions in the low-energy bands of a twisted Moiré bilayer with twist angle to leading powers in , where is the lattice unit of graphene, and is the Moiré lattice unit. The charge distribution of the states in the bands near the neutrality point is mostly concentrated in regions with stacking. This inhomogeneous charge leads to electrostatic potentials of strength . For an angle of the order , , and a Moiré pattern of wavelength nm, we obtain a Coulomb potential in the order of .

The overlapping orbitals of the carbon atoms also generate an intraatomic Hubbard interaction, , which prevents double occupancy. This potential, when projected onto wavefunctions of extensions similar to , generates an effective interaction , where is the number of atoms within the Moiré lattice unit, and the factor originates from the normalization of the wavefunctions. As a result, . Thus, the intraatomic Hubbard term is comparable to the electron bandwidth at the magic angles where meV, and the long-range Coulomb interaction can be larger. In the following, we only consider the effect of the long-range part of the Coulomb interaction. A detailed analysis of the short-range Hubbard repulsion can be found in ref. 18.

Band Structure and Charge Distribution

For small twist angles, the wavelength of the Moiré superlattice is much larger than the graphene unit cell. It therefore makes sense to analyze the band structure by using the continuum approximation; refs. 12–14 and 16. We specifically use the parametrization of ref. 22 for a Moiré pattern with a twist angle . Numerical values are given in . The bands and the charge density at high symmetry points in the Brillouin zone are shown in Fig. 1. As extensively discussed in the literature, the charge density associated with states at the edges of the Brillouin zone, the and points, is concentrated in regions with stacking. On the other hand, at the center of the Brillouin zone, the point, the charge density vanishes in these regions and is much more homogeneously distributed throughout the rest of the Moiré cell. This asymmetry between charge densities in different regions of the Brillouin zone has been pointed out in ref. 37. As the bands in the two valleys are related by time-reversal symmetry, we will consider in most of the following analysis a single valley, taking into account the fourfold spin and valley degeneracy when necessary.

Fig. 1.

(A) Low-energy bands of a twisted bilayer with twist angle plotted over the full Brillouin zone (see also Fig. 2). (B) Charge densities of states at the point in the Brillouin Zone (t for top layer; b for bottom). (C) The charge density at the point summed over the two degenerate states. (Scale bars: .)

Fig. 2.

(A) Dependence of the density of states on the value of the Hartree potential for a twisted bilayer graphene with a twist angle . The blue line is for and the red one for . The vertical scale is truncated: The red peak extends more than twice as high as shown. The dotted lines are the Fermi energies at neutrality. (B) Fourier transform of the charge density in the lowest state above the Fermi energy, with the dominant peak at zero momentum suppressed. The area of the circles is equal to the magnitude of the Fourier components for each of the two layers. The momentum distribution is concentrated at the six lowest wavevectors, . (C) The value of the charge density at wavevector as a function of the filling relative to charge neutrality. The blue line gives the real part, and the yellow line gives the imaginary part.

Hartree Potential

The Hamiltonian of twisted bilayer graphene breaks electron-hole symmetry, and charge inhomogeneities exist, even at half filling. These fluctuations are expected to be small, as, for low twist angles, the Moiré pattern can be decomposed into regions with , and stacking, which are locally neutral. In the following, we neglect these charge inhomogeneities. The charge distribution which arises away from the neutrality point in inhomogeneous, and it is peaked in the regions, so that a change in the occupancy of these bands will lead to charge fluctuations with the periodicity of the Moiré unit cell, , as confirmed by numerical calculations described below. Fluctuations of the charge density and electrostatic potential are dominated by the first star of reciprocal lattice vectors, the six points with .

If we consider partially filled conduction bands, their contribution to the charge density can be parameterized by a dimensionless number, , a complex quantity which encodes the symmetric (real) and antisymmetric (imaginary) parts of the charge density. From the value of , we can easily find the electrostatic potential in Fourier space, , and . The amplitude of the fluctuations of the potential in real space are given by , where is the area of the Moiré unit cell. We finally obtain . In the following, we neglect the dependence of on sublattice and layer, as both the density and the potential vary over length scales of order . The value of depends on the extent of the wavefunctions in momentum space, and it is bounded by , as it is a cross-product of amplitudes times a spin-layer degeneracy factor 4,Here, is a constant which takes into account the contribution to the total density from all of the bands not included in the calculation. We fix its value by imposing a homogeneous state at charge neutrality, , where is the density at half filling. The coefficient is the amplitude that the wavefunction of band state is represented by a Bloch state with momentum , and a sum over sublattice and layer indices is omitted. The wavefunctions of the lowest bands in twisted bilayer graphene are delocalized in momentum space, and we find . In all cases considered, . This imaginary part arises from the small breaking of the symmetry between the and the and regions in the Hamiltonian (22). We therefore describe the Hartree potential in terms of its real part only.

We approximate the fully self-consistent Hartree potential by describing it in terms of the six shortest reciprocal lattice vectors. Our calculations indicate that the next set of higher Fourier components are an order of magnitude smaller. The amplitude of the Hartree potential, , determines the parameters of the Hamiltonian. As function of these parameters and the filling, , we determine the charge density, . The self-consistency equation for the Hartree potential can be written aswhere is the band filling.

In Fig. 2, we show the density of states for the original problem and also for a Hartree potential of amplitude meV. The Hartree potential distorts the bands considerably. The Fourier transform of the charge as function of band filling, shown in Fig. 2, looks very similar for both cases—the difference is just in a normalization factor, and we only give the example for .

As we increase the Hartree potential, the energies at the and points are raised in comparison with those at the point, and the upper band becomes flatter (as can be seen in , this effect is largely described by first-order perturbation theory). We show a selection of bands for a few choices of the Hartree potential in Fig. 3. The Hartree potential is determined by the dimensionless parameter , which by definition vanishes at the neutrality point. Its value shows an approximately linear a dependence on the fractional band filling, , , as can be seen in Fig. 4.

Fig. 3.

Bands calculated for different electrostatic potentials. The blue lines are the bands at charge neutrality; the green lines are at ; and the orange, red, and dark-red lines are at , and 2 meV, respectively.

Fig. 4.

Dependence of the real part of on band filling for a twisted bilayer graphene with Hartree potentials described by different values of the amplitude . The band filling is normalized such that corresponds to the half-filled case, and to the case where the four upper bands are filled. We use values of from 0 to 24 meV.

For and nm, we obtain meV. Using this value, the solution to Eq. is shown in . The value of grows approximately linearly with . The value of changes from 0 at the neutrality point, to when all of the low-energy bands are occupied.

Exchange Potential

We now consider the Fock part of the Hartree–Fock approximation. Exchange coupling takes place between states in the same valley with identical spin. The exchange potential plays an important role in monolayer graphene, as it is the origin of the renormalization of the Fermi velocity (38, 39). Unlike the case of the Hartree term considered above, the effect of the exchange potential will lead to nontrivial changes in the bands at the neutrality point. The small value of the Fermi velocity, , in twisted bilayer graphene implies that the effective fine structure constant associated with the Dirac cones, is very large, . The renormalization of has a small value for the high-energy cutoff, which is of the order of the bandwidth . The renormalization of depends logarithmically on and linearly on , so that the corrections to will be important. The main effect will be a depletion of the density of states at the Dirac point.

Besides the renormalization of the Fermi velocity, the exchange potential leads to an overall widening of the bands. We can estimate this effect by considering the shift of the states at the band edges, located at the point. As the Coulomb potential is singular at zero momentum, the leading contribution to the shift of the states at the point comes from the interaction with occupied states with similar momenta. The occupied states at the bottom of the band have a similar internal structure to the neighboring ones, differing only in the crystal momentum. If we assume that this approximation is valid up to some cutoff in momentum space at a distance from the point, the energy shift has the valueAssuming that , we find that the shift of the lower state and the bandwidth scale as . We have analyzed numerically the decay of the overlap as a function of the distance between and . Results are shown in . They show that . We thus obtainThis energy is an approximation for the increase in bandwidth. As the unoccupied state at the top of the conduction band has an internal structure which is orthogonal to the states at the lower point, it will not experience a significant exchange shift.

Description of Electrostatic Effects in Terms of Local Interactions

The Hartree potential significantly distorts the lowest band in a twisted bilayer graphene, but it leaves the bands which lie further away from the Dirac point mostly unchanged. Hence, it can be expected that the Hartree potential does not significantly alter the Wannier wavefunctions for these bands. These Wannier functions have been extensively discussed (e.g., refs. 21–23). For a given valley, the two Wannier wavefunctions can be approximated by functions centered in the and regions of the Moiré pattern. These regions form a honeycomb lattice. Each Wannier function has a three-lobe structure with lobes peaked at the regions. The Hartree potential can be projected onto these wavefunctions, and it can be written in terms of diagonal and off-diagonal matrix elements. Each Wannier function is a four-component spinor, with weight on both sublattices in each of the two layers; that is,where label the sublattice and valley degrees of freedom. Each of the components has a different symmetry around the lobes, with angular momentum ; ref. 22. The Hartree potential is diagonal in sublattice, layer, and spin space. The matrix elements where label the six sites around a given region, can be written asThe functions belong either to the or the sublattice in the emergent honeycomb lattice with the Moiré pattern as a building block. We choose the sites to belong to the sublattice, and the sites to belong to the sublattice. The angular momentum of the different components of each Wannier function determines the phase of the amplitudes in Eq. . Matrix elements between functions which belong to different sublattices require the evaluation of integrals between states with angular momentum and states with . The phases shown in ref. 22 suggest that these elements vanish. The remaining matrix elements describe a second nearest-neighbor hopping term in the Moiré lattice. Hence, we can add the Hartree potential to the simplified two-parameter tight binding model described in ref. 22 to obtainwhere and are real-valued tight binding parameters, parameterizes the Hartree potential, and and are first and second nearest neighbors. This Hamiltonian ignores a number of the hopping terms considered in refs. 21–23, but the terms included are sufficient for a discussion focused on the role of the Hartree potential. Results for different values of are shown in Fig. 5. The bands are in reasonable agreement with the bands obtained using the continuum model, shown in Fig. 3.

Fig. 5.

Band structure of a simple tight binding model based on the Wannier functions of twisted bilayer graphene, including the electrostatic potential. Black lines, bands at charge neutrality; red and blue lines, bands as function of increasing Hartree potential.

The description of the Hartree potential in terms of effective hopping parameters allows us to define an Wannier-based model for the long-range Coulomb interaction in graphene. The occupation of Wannier functions leads to charge accumulation at the sites. The strength of the potential at a given site is proportional to the occupancy of the six Wannier functions centered at the neighboring and sites. The ensuing potential leads both to a change in the on-site energies and to the creation of effective couplings between Wannier orbitals which are second nearest neighbors in the honeycomb lattice, but overlap in the same region. We obtain the effective Hamiltonianwhere the sum over the index implies a sum over the centers of the hexagons in the honeycomb lattice, which define the sites. The first term in describes the local repulsion between charges placed at the sites, and it has been discussed in ref. 32. The second term in describes a hopping which depends on the charge state at the regions. The description of the Hartree potential depends only on the Coulomb interaction, parameterized by . Hence, we expect that both terms in scale as . The relation between and can be inferred from the band structures shown in Fig. 3. The Hartree potential used to calculate the bands, , gives the scale of the on-site term in in Eq. , parameterized by , while the relative shift between the and points is due to the assisted hopping term, . The difference in the Hartree potential between the and regions in the Moiré cell is , which gives an upper bound to the on-site energy described by (note that the weight of a Wannier on a given region is ). On the other hand, the shift of the point with respect to the point induced by the hopping term in Eq. is . From the results in Fig. 3, we conclude that .

Electron-Assisted Hopping and Superconductivity

The existence of assisted hopping terms (40) in the effective Hamiltonian of twisted bilayer graphene seems natural, as different Wannier functions overlap at regions where charging effects are maximal (21–23). The effect of assisted hopping interactions in high- superconductivity has been studied in refs. 41 and 42. The most likely instability favored by such an interaction is superconductivity, since (i) the effective interaction is attractive for some range of fillings, and (ii) an assisted hopping term does not lead to many of the other typical broken-symmetry phases, such as a magnetic state, or a charge-density wave.

We analyze the effect on superconductivity of the interaction term in Eq. using the Bogoliubov–de Gennes approximation. For simplicity, we consider here only superconducting states with spin zero and invariant under time-reversal symmetry. The competition between the repulsion described by the term proportional to and the assisted hopping term, , can lead to many other superconducting phases. The local description derived from Eq. is not strictly necessary for this step; . The electron and hole Hamiltonians are described in momentum space by two matrices, one for electrons with a given spin and valley index and another for holes with opposite spin and valley index. The mean field decoupling of the interaction term leads toHere, and are the identity matrix and a Pauli matrix in sublattice space, and operate on the electron hole index, and is the Fermi energy. The vectors are the lattice vectors of the honeycomb lattice. The values of have to be determined self-consistently; .

The existence of superconductivity is determined by the functions and in defined in . The function is always positive. The terms in Eq. proportional to describe a repulsive interaction which cannot lead to superconductivity in an isotropic s-wave-like channel. This repulsive interaction is somewhat underestimated, as interactions over length scales longer than are ignored (see, however, the estimate of the screening length in ). On the other hand, the function is positive near the point and negative in a large region near the edges of the Brillouin zone. If , superconductivity becomes possible. Fermi surfaces for different fillings are shown in Fig. 6 and compared with the region in the Brillouin zone, where .

Fig. 6.

Self-consistent Fermi surface for different band fillings. The gray areas denote the filled states. The function in is negative in the blue region, where superconductivity is favored. (A), (B), (C), (D), (E), and (F). The magnitude of gives the fraction of bands filled. The value denotes the neutrality point, and describes four filled bands. The Fermi surfaces are calculated by using the full continuum Hamiltonian, including the Hartree term.

The results in Fig. 6 show that, typically, there are two pockets at the Fermi surface and that, for a large range of fillings, superconductivity is possible in one pocket but not in the other. As a result, the superconducting state will show two different gaps, or a gap-less pocket coexisting with a gapped one. The order of magnitude of the higher superconducting gap will be , where is the bandwidth. Repulsive interactions at the atomic scale not considered here, such as an on-site Hubbard term, can suppress the superconducting phase at integer fillings (18).

Conclusions

The occupation of the low-energy bands in a twisted bilayer graphene leads to inhomogeneous electrostatic potentials of magnitude of order of , where is the Moiré unit length. This estimate is comparable, or larger, than the width of the band.

Electrostatic effects, induced by charging the system away from the neutrality point, distort the bands significantly. The states at the edges of the Brillouin zone, at and , are shifted with respect to the states near the center of the Brillouin zone, the point . The exchange term, on the other hand, leads to an increase in the bandwidth, approximately a fraction of .

The band distortion induced by the electrostatic potential can be described in terms of induced assisted hopping couplings. These terms fit naturally with the complex overlapping Wannier functions which give a local description of twisted bilayer graphene. Assisted hopping interactions favor generally superconductivity, and we explicitly show that -wave pairing is possible at certain fillings.

Note.

After the submission of this manuscript, ref. 40 was posted, also analyzing assisted hopping terms in graphene.