Electrical band flattening, valley flux, and superconductivity in twisted trilayer graphene

Twisted graphene multilayers have demonstrated to yield a versatile playground to engineer controllable electronic states. Here, by combining first-principles calculations and low-energy models, we demonstrate that twisted graphene trilayers provide a tunable system where van Hove singularities can be controlled electrically. In particular, it is shown that besides the band flattening, bulk valley currents appear, which can be quenched by local chemical dopants. We finally show that in the presence of electronic interactions, a non-uniform superfluid density emerges, whose non-uniformity gives rise to spectroscopic signatures in dispersive higher energy bands. Our results put forward twisted trilayers as a tunable van der Waals heterostructure displaying electrically controllable flat bands and bulk valley currents.


I. INTRODUCTION
The interplay between topology and correlations represents a highly fruitful area in condensed matter physics. However, exploring unconventional states of matter requires identifying systems where electronic correlations, topology and electronic dispersions can be realistically controlled. In this line, twisted van der Waals materials 1-8 provide a powerful solid state platform to realize exotic quantum phenomena. The tunability of twisted van der Waals materials stems from the emergence of a band structure that can be controlled by the twist between different two-dimensional materials. 9,10 In particular, the different quantum states in twisted graphene systems stem from the possibility of controlling the ratio between kinetic and interaction terms. In twisted graphene bilayers, such tunability allowed to realize superconducting, 1,2,4 correlated insulators, topological networks, 5,6 Chern insulators 11 and quasicrystals. [12][13][14][15] As a result, current experimental efforts are focusing on exploring new twisted van der Waals materials, with the aim of finding platforms that allow for an even higher degree of control. 16,17 From the quantum engineering point of view, applying a perpendicular bias between layers 18 provides a versatile way of tuning correlated states in twisted graphene multilayers. This has been demonstrated in paradigmatic examples of correlated states in twisted tetralayers (double bilayer) 19,20 and twisted trilayers (monolayer/bilayer). 16,17 Moreover, interlayer bias is known to generate internal valley currents in twisted graphene bilayers, 5,6,21,22 creating topological networks at low angles 5,6,21 and generating valley fluxes in flat bands regimes. 22 This interplay of correlations and topology in twisted graphene multilayers makes these materials a powerful platform to explore exotic states of matter [23][24][25][26] in a realistically feasible manner.
From the theoretical point of view, electronic structure calculation of twisted graphene bilayer conducted with real-space tight-binding models 10 or continuum Dirac  1. (a,b) Sketch of the structure a twisted graphene trilayer, in which the top and bottom layers are aligned, and twisted with respect the middle one. Panel (c) shows the color map of a DFT-based relaxed 2×2 supercell of twisted trilayer. Colored pattern corresponds to the difference ∆ = z − zav between the z coordinate of an atom in surface layer and the average coordinate zav of all atoms in the same layer. AA stacking regions are colored with darker tones whereas AB/BA stacking regions are colored with lighter tones. Superimposed is the unit cell of the trilayer. Panel (d) shows the histogram shows the number of occurrences of each distance to the average deviation value. Tiny deviations from flatness occur in the central layer.
descriptions 9 capture the fundamental features of the electronic dispersion. Nevertheless, internal coordinates optimization can quantitatively modify the electronic dispersion. [27][28][29][30][31][32][33][34] Well known examples of this are the growth of AB/BA regions in twisted bilayers. 29,35 It is important to note that studying twisted graphene multilayers from first-principles represents a remarkable challenge, due to the large amount of atoms present in a unit cell.
Here, by combining first-principles calculations and low-energy models, we show that twisted graphene trilayers 36 host flat bands whose bandwidth can be controlled electrically. We address the impact of an interlayer bias both from first-principles and effective models, showing that van Hove singularities can be merged electrically. We show that associated with the interlayer bias, bulk valley currents emerge, that are impacted by the existence of chemical impurities in the system. We finally address the superconducting states in these doped trilayers, showing that the non-uniform superfluid density has an impact in high-energy bands. Our manuscript is organized as follows, in Sec. II we show the electronic structure of twisted graphene trilayer both from firstprinciples and low-energy models, in Sec. III we explore in detail the impact of an electric field, in Sec. IV we explore the effect of chemical impurities, and in Sec. V we address the impact of an emergent non-uniform superfluid density. Finally, in Sec. VI we summarize our conclusions.

II. ELECTRONIC STRUCTURE OF TWISTED TRILAYER GRAPHENE
The electronic structure of twisted graphene trilayers [36][37][38][39][40][41] shows different features in comparison with twisted graphene bilayers. 10,42 The electronic structure of small "magic" angle tBLG features four flat bands lying around the Fermi level, with band splitting at the Γ point of the Brillouin zone of the emergent moiré superlattice. 10,42 Density functional theory (DFT) calculations have shown consistent results with those effective models in twisted graphene multilayers, yet quantitative modifications are observed when including relaxation of the atomic coordinates 29,30,43 and crystal-field effects 34,44 . Therefore, to benchmark the electronic properties of twisted graphene multilayers, it is essential to start from a correct description that takes into account the geometric corrugation and ab-initio electrostatics of the moiré system.
We consider a twisted trilayer structure in which the upper and lower layers are aligned, and the middle one is twisted with an angle θ with respect to those (Fig. 1ab). In particular, in the following we will consider a twisted trilayer whose middle layer has a twisting angle of 1.9 • with respect to the external layers. The system of 5514 atoms is fully relaxed allowing for lateral and vertical displacement of the C atoms in the structure. Figure 1c shows the color map of one of the two equivalent external layer relaxation. The color scheme represents the vertical variation of each C atom at the surface with respect to the average deviation within each layer. The darker areas indicate a displacement of atoms out of the surface which occurs predominantly in the AA stacking region. The histogram of Fig. 1d shows that the number of atoms in the upper layer whose vertical coordinate is above the average is twice as large as the atoms displaced in the opposite direction below the average.
The first-principles electronic band diagram of the fully relaxed structure of θ = 1.9 • TTG is shown in Fig. 2. In contrast with twisted bilayer graphene, two highly dispersive bands coexist with four low-dispersive bands grouped at the Fermi energy. The Dirac-like crossing above the charge neutrality point leads to a small charge transfer between the flat and dispersive bands even at half filling. Additional flattening of the localized states can be induced by application of an external electric field perpendicular to the TTG surface. A field of 0.03 eV/Å reduces the dispersion of all electronic states in the vicinity of the Fermi energy without inducing any inter-band charge transfer. Increasing the strength up to 0.25 eV/Å a disruption of the linear bands is observed and a hybridization of the flat bands with neighboring state increase their dispersion. It is observed that the application of an interlayer bias generates Dirac crossings above and below charge neutrality, besides a variety of anticrossings (Fig. 2). The first-principles calculations above show that the electronic structure of twisted trilayer graphene show strong differences with the one of twisted graphene bilayers. In particular, a highly dispersive set of bands coexists with the nearly flat bands at charge neutrality. In order to explore more in detail the physics and twisted bilayers, in the following we will exploit a lowenergy model. We take a single orbital per carbon atom, yielding a tight-binding Hamiltonian of the form where d is the interlayer distance and β controls the decay of the interlayer hopping. As a reference, for twisted graphene multialyers t ≈ 3 eV and t ⊥ ≈ 0.15t. 45 Similar realspace models were used to study a variety of twisted graphene multilayers, 10,36,46,47 providing a simple formalism to study the effect of dopands and impurities. 48,49 However, in contrast to continuum models, 9,42,50 measuring of valley related quantities with a real-space based formalism is non-trivial.
Twisted graphene multilayers own an approximate symmetry associated to the valley quantum number. 9 Valley physics in both DFT calculations and tightbinding are emergent symmetries, in the sense that valley are not easily defined in terms of real-space chemical orbitals. This limitation can be overcome by defining the so-called valley operator 21,48,51 in the tight-binding description. With the valley operator the expectation value of the valley can be computed in a real space representation. 52,53 The valley operator can be defined as 21,48,51 Electronic band diagrams of twisted trilayer graphene. Localized states at the Fermi energy coexist with highly dispersive states bands forming a Dirac cone in the conduction band. Applying a small electric field perpendicular to the TTG a reduction of all bands dispersion is induced. A stronger electric field removes the Dirac cones and induces hybridization between electronic states.
where i, j denotes second neighbor sites, η ij = ±1 for clockwise or anticlockwise hopping, and σ ij z is a Pauli matrix associated with the sublattice degree of freedom. With the previous operator, we can compute the valley flavor of each eigenstate of the twisted trilayer supercell within the real-space formalism as V z = Ψ|V z |Ψ |. It is worth noting that this operator can be easily defined in the tight-binding basis but not in the DFT basis.
With the previous formalism, we now compute the electronic structure of the low-energy model at the same angle θ = 1.9 • as in the first-principles calculations (Fig.  3). We note that, besides some additional splittings observed in the first-principles calculations (Fig. 2), the band-structure obtained with the low-energy model (Fig.  3ab) gives comparable results. As shown in Fig. 3a, in the absence of an interlayer bias the system shows nearly flat bands coexisting with highly dispersive states. Whereas the nearly flat bands are degenerate in valleys in the path shown, the dispersive states belong to different valleys in different parts of the Brillouin zone. It is also observed that the dispersive Dirac cones are slightly displaced from charge neutrality, as observed in the first principles results.
When an interlayer bias is turned on (Fig. 3b), the valley degeneracy in the G − K path is lifted, and the dispersive bands show two crossings above and below the nearly flat bands, analogously to the first-principles results. It is also observed that the nearly flat bands are slightly modified by the interlayer bias. As we will show below, such effect becomes stronger for even smaller angles.
A last interesting point is related with the the localization of the different states in the supercell. This can be characterized by means of the inverse participation ratio (IPR), defined as IP R = i |Ψ k (i)| 4 . Large values of the IPR correspond to states localized in the moiré supercell, whereas low values correspond to delocalized in the moiré supercell. Focusing now in a structure with θ = 1.6 • twisting angle, it is clearly observed that the nearly flat bands show a substantially higher degree of localization than the dispersive bands (Fig. 3c). In particular, the flat band states are associated to an emergent triangular lattice in the supercell, as shown in the local density of states (LDOS) of Fig. 3d. In the following we will see how these flat band states can be electrically controlled, and how the bias creates bulk valley currents associated with the valley splittings.

III. BAND FLATTENING AND VALLEY CURRENTS BY AN INTERLAYER BIAS
We now move on to systematically analyze the effect of an interlayer bias in the twisted graphene trilayer, and in particular its effect on the low energy density of states. Electric biases in twisted graphene multilayers are known to give rise to valley Hall currents, 5,6,21,54-57 and represent an effective knob to control the low-energy electronic structure. 18 An interlayer bias can be easily included in the low-energy tight-binding model by means of where z i is the z-position of site i, d is the interlayer distance and V is the strength of the interlayer bias. The full Hamiltonian in the presence of interlayer bias is thus In the following we focus on the structure at θ = 1.6 • . As shown in Fig. 4a, the application of an interlayer bias merges the original two van Hove singularities of zero bias in a single one, dramatically enhancing the low-energy density of states. Interestingly, this merging of van Hove singularities is similar to the evolution with the twist angle in twisted bilayers, 10,42,58 with the key difference that in the present case the merging is electrical, turning the process highly controllable in-situ.
We now focus on the flattest regime with finite bias, when the two van Hove singularities get merged. Fig.  4bcd shows the band structure of that regime . It is observed that almost perfectly flat bands are present at the Fermi energy (Fig. 4bcd), accounting for the van Hove singularity at charge neutrality (Fig. 4a). It is also observed that the states remain highly localized in the unit cell as highlighted by the IPR, whereas at higher energies the dispersive bands are delocalized in the moiré unit cell (Fig. 4b). It is also observed that the bands show perfect valley polarization (Fig. 4c), so that the interlayer bias does create any intervalley scattering in agreement with continuum models. Given that the interlayer bias breaks the symmetry between the top and bottom layers, it is interesting to look at the layer polarization of the states L , defined by the layer polarization operator L = i,s z i c † i,s c i,s , with z i the z-component of the site i. As shown in Fig. 4d, the Dirac cones at the K, K points above and below the flat bands have a slightly opposite layer polarization. However, such layer polarization is reversed at the Γ point (Fig. 4d), highlighting that the states remain highly entangled between all the layers.
We now explore how the band flattening is accompanied with the emergence of bulk valley currents. The emergence of valley currents associated with interlayer biases is a well known effect in aligned graphene bilayers 54-57 and tiny-angle twisted bilayers. 5,6,21 Such valley currents arise due to the emergence of valley a non-zero local valley Chern numbers, 59 whose quantization is associated to the emergent valley conservation. 60 In twisted systems, the local Chern number is expected to change from region to region, 6 due to the locally modulated Hamiltonian. To address this local topological property, in the following we compute the local valley flux by means of the Berry flux density χ V . The realspace valley flux χ V defines the valley Chern number as C V = C K − C K = χ V (r, ω)d 2 rdω. The real-space valley flux χ V can be computed as 52,53,61 (4) where αβ denotes the Levi-Civita tensor, G V the valley Green's function G V = [ω −H k +i0 + ] −1 V z , H k the Bloch Hamiltonian, and V z the valley polarization operator of Eq. 2. Fig. 4f shows the spatial profile of the valley flux density at the Fermi energy χ V (r) ≡ χ V (r, ω = F ), for a chemical potential with one hole per unit cell. It can be observed clearly that sizable valley currents appear in regions in the complementary regions to the low-energy states, similarly to other moiré systems. 52,53 . Since the emergence of such currents relies on valley conservation, terms in the system creating intervalley mixing are expected to substantially impact them, 6,62-64 as we address in the next section.

IV. IMPACT OF CHEMICAL DOPING
In this section we address the modification of the electronic structure under both electrical and chemical dop- ing. First, we start with the study of the electric doping from first-principles, and afterwards we move on to consider its effect in the low-energy effective model.
We first address the impact of doping from firstprinciples. First-principles DFT calculations are employed to study the modification of the band diagram of TTG upon doping with both an extra charge and one of the layers with foreign species. Uniform doping of the TTG is realized by introducing an extra electron which is compensated with an equally uniform background charge of opposite sign. This can be realized experimentally via field effect. Figure 5 shows that the extra electron fills the flat states and align the Dirac point with the Fermi energy. Small modifications in the electronic band structure as a result of new electrostatic contributions is observed, similarly to other twisted bilayer systems. [65][66][67] A similar effect can be induced with chemical doping, namely adding one N atom in substitution of a C atom. The effects of chemical substitution have been extensively studied in graphene, [68][69][70][71][72][73][74][75][76][77] including twisted bilayers. 49,78 Figure 5 shows the effect on the bands structure of one N atom in one of the surface layers. The extra charge supplied the N atom shifts the chemical potential, increasing the filling of the flat bands. Major difference with respect to the electrostatic doping is the opening of a meV large band gap as a result of the symmetry breaking imposed by the impurity which also induces a mixing of the linear states with the less dispersive states. These new anticrossings are due to the intervalley scattering created by the chemical impurity, and they do not appear in the case of electrostatic doping. A similar hybridization of electronic states is also observed when one electron is removed by means of doping with a B atom, as shown in Fig. 5. The chemical potential shift occurs in the opposite direction and the filling of flat band states

decreases.
We now consider the effect of a chemical impurity in the tight-binding model. We model the addition of a chemical impurity by adding to the Hamiltonian H D = w s c † i,s c i,s , where i in the site that has been chemically replaced. We take w = −2t, which is the typical energy scale expected for a N dopand, and yields results comparable with the first principle calculations. For the sake of concreteness we will focus on the case with interlayer bias, so our full Hamiltonian will be H = H 0 + H V + H D . With the previous Hamiltonian, we now compute the electronic band structure and project each eigenstate onto the valley operator. The result is shown in Fig. 6a, where we see that small anticrossings appear as in the first-principles calculations. The valley projection clearly shows that such anticrossings are associated with intervalley mixing. The effect of the impurity in terms of intervalley mixing can also be readily seen in the bulk valley currents. In particular, the existence of the impurity is expected to strongly perturb the original valley fluxes in the unit cell. This is verified in Fig.  6b, where we observe that the original valley fluxes are impacted by the presence of the impurity.
Interestingly, despite the effect in terms of intervalley mixing, the low-energy bands remain relatively flat, retaining their associated large density of states. This suggests that correlated states can still appear in chemical doped twisted trilayers. It is worth to emphasize that despite this large density of states, correlated phases relying on valley coherent states will be strongly suppressed due to the impurity-induced intervalley mixing. In particular, valley ferromagnet states, and valley triplet superconducting states will be depleted due to chemical dopands. Nevertheless, conventional spin-singlet valleysinglet states are not affected by the presence of intervalley scattering. Motivated by this, in the next section we address the emergence of spin/valley singlet supercon-ductivity, and show how the superfluid density impacts the high energy dispersive states.

V. SUPERCONDUCTING STATE
The large density of states close to charge neutrality suggests that the twisted graphene trilayer can have superconducting instabilities, similarly to twisted bilayers and tetra-layers. As shown above, both with electrostatic and chemical doping the system shows divergent density of state close to charge neutrality. For the sake of concreteness we will now focus on the electrostatically doped system, yet we have verified that our results remain qualitatively unchanged with chemical doping.
An emergent superconducting state is associated with a Fermi surface instability, yet its effect can give rise to second order perturbations above the Fermi energy. In order to understand the potential impact at highenergies, we first briefly analyze the structure of the low-energy states. This can be done by comparing the spatial distribution of the low-energy states with respect to the Fermi surface states. In particular, we define the projection over the Fermi surface states as R = d 2 r|Ψ(r)| 2 ρ F (r) with ρ F (r) the local density at the Fermi energy ρ F (r) ∼ r|δ(E F − H)|r . The quantity R allows to qualitatively distinguish which states are localized in the same region at the Fermi surface states.
By computing the Fermi surface state projector R, it observed that the states of the flat band are localized in similar regions. In contrast, as one departs from the Fermi energy, the states start to delocalize to other regions of the unit cell (Fig. 7a). This highlights the different orbital nature of the nearly flat bands (in purple) and low-energy dispersive bands (in green). Importantly, and in strike contrast with twisted graphene bilayers, the flat bands of this system are not decoupled from the dispersive states, suggesting that the superconducting states of this system will have a genuine multiorbital nature.
Given the unavoidable entanglement between the flat and dispersive bands, an effective model description of this twisted trilayer cannot be easily performed. Therefore, in the following we will study the emergent superconducting state by exploiting the full atomistic model, including all the bands in our calculation. A variety of mechanisms have been suggested to give rise to attractive interactions in these systems, including phonon, [79][80][81] Coulomb 82 and magnon fluctuations. 83,84 For the sake of concreteness, we will focus on effective local attractive interactions, as employed in other twisted graphene systems 79,80,[85][86][87][88][89][90] that we solve at the mean field level H M F H V +H M F I , that we solve with the Bogoliubov-de Gennes (BdG) formalism. The local attractive interaction g will give rise to a net superfluid density, with a moiré momentum structure of s-wave symmetry. In the following we take g = t, and we verified that our results remain qualitatively similar with smaller interaction strengths. We note that although interactions are local, the would lead to a non-trivial multiorbital structural in the moiré orbital space.
By solving the previous selfconsistent problem, we find that the superfluid density is non-uniform in the moiré unit cell (Fig. 7b), stemming from the non-uniformity of the low-energy states. By projecting the BdG eigenstates in the electron sector via the electron projector P e , we observe that the net non-uniform superfluid density ∆ ↑↓ (r i ) ∼ c † i,↑ c † i,↓ gives rise to a full gap in the Brillouin zone in the superconducting state. This phenomenology is similar to the one found in twisted graphene bilayers. More interestingly, besides the gap opening at the chemical potential, anticrossings appear in high energy bands when the selfconsistent pairing is included (Fig.  7c). It is worth to emphasize that, in the presence of a uniform pairing artificially imposed, the high energy anticrossings disappear, leading only to the gap opening at charge neutrality (Fig. 7d). The emergence of gap openings away from charge neutrality is associated to the intrinsically multiorbital nature of the superconducting state, and stems from the non-unitarity of the superconducting matrix in orbital subspace. 91,92 Interestingly, this shows that signatures of the superconducting state can be obtained by analyzing the system away from the chemical potential, and could provide powerful spectroscopic signatures 93,94 of the superconducting state.

VI. CONCLUSIONS
By combining first-principles calculation and lowenergy effective models, we have shown that twisted graphene trilayers realize tunable electronic systems. In particular, is was shown that nearly perfect flat bands can be electrically controlled, that coexist with highly dispersive states. Interestingly, such electric flattening of the bands is accompanied by the emergence of bulk valley currents. We have found both from first principles and low-energy calculations that chemical doping does not destroy the flat bands, yet it substantially impacts the bulk valley currents. This suggests that chemical doping of twisted graphene bilayers could provide an intrinsic way of providing the necessary electronic doping required for the emergence of a superconducting state. We finally demonstrated that an emergent superconducting state would give rise to spectroscopic changes int he high energy bands, associated to the non-uniform superfluid density. Our results highlight the rich physics of twisted graphene trilayers, and provide a starting point to explore the interplay between flat bands, correlation and dispersive states in twisted graphene multilayers. Description of the coupling between the graphene layers was conducted through self-consistent calculations with the SIESTA code 95 within a localized orbital basis set scheme. Paramagnetic calculations were conducted using a double-ζ basis set, and the local density approximation (LDA) approach 96 for the exchange-correlation functional was used. Atomic positions of systems formed by over 5514 atoms were fully relaxed with a force tolerance of 0.02 eV/Å. The integration over the Brillouin zone (BZ) was performed using a Monkhorst sampling in Γ point. The radial extension of the orbitals had a finite range with a kinetic energy cutoff of 50 meV. A vertical separation of 35Å in the simulation box prevents virtual periodic parallel layers from interacting.