Defect-induced band restructuring and length scales in twisted bilayer graphene

We investigate the effects of single, multiple, and extended defects in the form of non-magnetic impurities and vacancies in twisted bilayer graphene (TBG) at and away from the magic angle, using a fully atomistic model and focusing on the behavior of the flat low-energy moir\'e bands. For strong impurities and vacancies in the $AA$ region we find a complete removal of one of the four moir\'e bands, resulting in a significant depletion of the charge density in the $AA$ regions even at extremely low defect concentrations. We find similar results for other defect locations, with the exception of the least coordinated sites in the $AB$ region, where defects instead result in a peculiar band replacement process within the moir\'e bands. In the vacancy limit, this process yields a band structure misleadingly similar to the pristine case. Moreover, we show that triple point fermions (TPFs), which are the crossing of the Dirac point by a flat band, appearing for single, periodic, defects, are generally not preserved when adding extended or multiple defects, and thus likely not experimentally relevant. We further identify two universal length scales for defects, consisting of charge modulations on the atomic scale and on the moir\'e scale, illustrating the importance of both the atomic and moir\'e structures for understanding TBG. We show that our conclusions hold beyond the magic angle and for fully isolated defects. In summary, our results demonstrate that the normal state of TBG and its moir\'e flat bands are extremely sensitive to both the location and strength of non-magnetic impurities and vacancies, which should have significant implications for any emergent ordered state.


I. INTRODUCTION
Twisted bilayer graphene (TBG) has attracted considerable attention as both a versatile tunable experimental platform and a host of a plethora of ordered states , including both superconductivity and correlated insulating states [2,7,23].Intriguingly, the wealth of ordered states in TBG is intimately connected to its remarkable unordered normal state electronic structure, where the Fermi velocity is suppressed with decreasing twist angle and even vanishes at so-called magic angles.The resulting flat dispersion has a large density of states and quenched kinetic energy that dramatically increases the importance of interactions and favors ordered states [29,30].As a consequence, around the magic angle, the ordered states of TBG depend crucially on the four spin-degenerate, emergent, low-energy, flat bands at the charge neutrality point (CNP), the moiré bands [31][32][33].Among their distinctive properties, the moiré bands have been shown to have a topological obstruction from nonlocal symmetries that impose a lower bound on the localization of their associated Wannier orbitals [4,34].Their unique energetics and spatial extent across the large emergent moiré pattern naturally prompts the question of how the moiré bands are impacted by explicitly local perturbations, such as atomic size lattice defects or impurities that are always present in any material.An answer to this question is important both in itself but also as it directly relates to the role of the moiré bands as a host of ordered states and as a versatile experimental probe.
When it comes to atomic size defects in single layer graphene, vacancies have been extensively studied and known to generate critically localized zero-energy states [35], exhibiting magnetism [36,37].Although isolated vacancies are not thermodynamically stable due to high formation energy, several mechanisms can still lead to their formation [38].Focused electron beams, for example, allow for the creation of vacancies with close to atomic precision [39].Stone-Wales reconstructions and double vacancy structures on the other hand, have lower formation energies and are therefore even naturally ubiquitous [38].Moreover, adatoms and substitutional impurities constitute other common and well-studied types of defects in graphene [38,40].As such, defect studies have become an integral part of studying graphene.Atomic size lattice defects have also been studied in TBG.Many of these studies have only focused on the large-angle regime, including both more comprehensive studies [41] and more focused studies into for example fluorination process [42,43], intercalation by lithium [44], and more general charged defects [45].Other, somewhat related, works include the study of impurity-induced Friedel oscillations [46], Raman spectroscopy of TBG samples with defects induced by ion beam irradiation [47], and the study of vacancies and their migration [48].Moreover, vacancies in TBG have been found generate similarly localized states as in graphene, including also leading to Yu-Shiba-Rusinov (YSR) magnetically induced sub-gap states in the superconducting phase of magic-angle TBG [49].Studies have also considered the interplay between defects and the pairing symmetries of the superconducting phase in magic-angle TBG [50,51].Finally, particularly relevant for this work is the finding of so-called triple point fermions (TPFs) [52][53][54][55][56][57][58][59], characterized by a triple band crossing, and associated valley polarization induced by single, weak, periodic impurities in TBG [60].
In this work we go beyond previous studies by providing a comprehensive investigation of the changes induced in the low-energy electronic structure of TBG by single, multiple, and extended atomic size lattice defects.By explicitly focusing on the impact low-energy electronic structure we both extract the inherent behavior of defects in TBG, which is likely notably different from graphene, and form a basis for understanding the implication of defects for emergent ordered states.To achieve this, we study TBG both near and away from the magic angle regime, for both weak and strong impurity strengths, and for both periodic and isolated defects.We primarily consider two types of lattice defects: non-magnetic potential impurities and vacancies, which model adsorbates, atomic replacements, and true vacancies and therefore effectively capture a wide range of different impurities and defects.We do so by employing fully atomistic tightbinding calculations including all carbon atoms and establish both the evolution of the band structure and the accompanied changes in the charge density as a function of impurity strength.
In the case of periodic defects, we introduce one or more defects to the moiré unit cell, the emergent unit cell of TBG, thus preserving the translational invariance such that the band structure is still well defined.To access instead the effects of isolated defects, we effectively separate the defects by using supercells comprised of many moiré unit cells, still within a fully atomistic approach.The analysis of the periodic and isolated defects therefore complement each other with different experimental relevance.While understanding the isolated defect case is important in itself and for example, for quasiparticle interference studies which can probe symmetries of ordered states [61][62][63], the periodic case has particular experimental and practical relevance based on contemporary impurity deposition techniques which enable the engineering of defect patterns [64,65], as has been also recently illustrated by the synthesis of periodic molecular arrays on graphene with both atomic precision and tunable periodicity [66].Additionally, the moiré unit cell itself defines an emergent periodic structure with a corresponding energy landscape that should intrinsically favor certain defect lattices, produce self-assembly of defects, or aid in the engineering of defect patterns.For instance, atomic hydrogen has been shown to preferentially adsorb following the moiré pattern produced between graphene and an Ir(111) substrate [67].The same moiré pattern has similarly been shown to produce regular Ir and Pt clusters [68][69][70].For TBG, ab-initio methods have already shown that atomic hydrogen preferentially absorb with a higher binding energy to the AA regions of the TBG moiré lattice [71].This suggests the possibility of structured patterning also of TBG and the engineering of periodic lattices of defects.
Our results show that atomic size impurities and vacancies have a profound and special effect on the lowenergy moiré band structure and thus directly on the properties of TBG.This is best illustrated by defects in the AA region, where we find a complete removal of one entire moiré band from the low-energy band struc-ture, even for only a single defect per moiré unit cell.At the magic angle this corresponds to an impurity concentration of only ∼ 0.01%, thus highlighting how extremely sensitive the low-energy electronic structure is to defects.This band removal results in a concomitant depletion of the AA lattice regions, where the moiré bands are primarily concentrated [72,73].With any interacting many-body ground state heavily dependent on the normal state low-energy moiré band structure, this defectinduced band removal will have severe consequences for the physics of TBG.We also observe an important dependence on the defect location.The defect-induced band removal also occurs for defects in the domain wall (DW ) region and for the higher coordinated sites of the AB region, while we instead find a band replacement occurring for defects in the least coordinated sites of the AB region.Thus, even if the degeneracy of the lowenergy moiré band structure is preserved in the last case, the properties of the bands are still completely altered.Moreover, we establish that the impurity strength necessary for an impurity to start behaving as a vacancy varies substantially for different impurity locations.Beyond the strong effect on the moiré bands we also find a localized defect state manifesting on the atomic length scale, which we trace back to the well-established defect state in monolayer graphene [35,37].The presence of this localized defect state is however easily obscured by the depletion of the AA region if the defect is also in this region.Thus there exist two length scales for defects in TBG: the atomic scale hosts a graphene-like defect state and the moiré scale controls the low-energy band structure and the change of charge density of the AA regions.Furthermore, we find that the defect-induced triple degeneracy at the Dirac point, generating TPFs [52][53][54][55][56][57][58][59], earlier found for single defects in TBG [60] is fundamentally not stable, but that the degeneracy is easily lifted with the introduction of either extended or multiple defects in the unit cell.We use degenerate perturbation theory to attribute this sensitivity of TPFs to an assumption of rank-1 perturbations, while more complex defect configurations generally violate such an assumption.
Our results establish that the low-energy electronic structure of normal-state TBG changes drastically with the introduction of non-magnetic defects, changes that should be taken into account when considering the influence of the moiré bands on any electronic ordering achieved at low temperatures.As we report both band structure and charge density, our results are experimentally easily verified using for example angle-resolved photoemission spectroscopy (ARPES) or scanning tunneling spectroscopy (STS) and transport measurements, and they can also be directly extended to quasiparticle interference experiments.Our results can also be straightforwardly used to engineer altered band structures, including changing the number of moiré bands, even completely removing all moiré bands, or introducing flat bands at the Dirac point, and thereby possibly generating very different electronic orders at low temperatures.This work is organized in the following way.In Sec.II we explain how we model TBG and its defects, including the observables used to examine the resulting electronic structure.Our results are presented in Sec.III.In Sec.III A we show that the electronic structure of magic angle TBG is generally extremely sensitive to the presence of periodic defects.We then classify the defect locations into two distinct cases in Sec.III B. In Sec.III C we complement earlier results by showing how the TPFs of TBG are broken by most realistic defects.In Sec.III D we consider isolated defects, which allows us to identify two distinct length scales of defect behavior.We subsequently show how this behavior is not restricted to the magic angle and hence showcases universal physics of TBG.Finally, in Sec.IV we draw our conclusions and propose ways in which our results could be tested experimentally.We also suggest implications of our work and pose questions for future investigation.

II. MODEL AND METHOD
TBG consists of two sheets of graphene stacked and rotated with respect to one another by a twist angle θ.As a consequence, the lattice is modulated by an emerging length scale, creating a moiré pattern.For small twist angles and if the rotation axis goes through a graphene lattice site in both layers, then the relative alignment between the layers around this site remains, locally and approximately, as AA layer stacking.Away from this site, the alignment transforms into AB or BA stacking depending on spatial direction.In between these directions, there is a transition region forming a domain wall (DW ), where the stacking does not belong to either of these classifications.These regions are schematically depicted in Fig. 1(a) where we plot one moiré unit cell.Note that because of periodic boundary conditions the four corners in Fig. 1(a) are connected and thus the AA region is split into the four corners of the moiré unit cell.The moiré unit cell forms a triangular lattice with lattice constant L m = a/(2 sin(θ/2)) [74], also called the moiré length, where a is the graphene lattice constant and θ is the twist angle.This lattice can be used to study the infinite system with atomic resolution, assuming appropriate commensuration conditions [74].We do this through a fully atomistic tight-binding model given by the Hamiltonian [75,76] where t ij are hopping matrix elements under Bloch boundary conditions and the sum is taken over the sites of the moiré unit cell, M. Near the magic angle, θ m ≈ 1.1 • , this amounts to considering on the order of 10 4 individual atoms.The operator c † ik creates an electron on a site i of the moiré unit cell.The layer, sublattice, and unit cell position of such a site are given by l i , s i , and ρ i , respectively.With δ li as the vector connecting A and B sites of layer l i and d 0 = 3.35 Å the interlayer distance, the site positions are given by where δ a,b is the Kronecker delta.We thus opt to use a rigid lattice model, ignoring lattice relaxation effects.In terms of lattice effects, this relaxation has been shown to cause the AA region to shrink, while in terms of the band structure it mainly rescales the magic angle and increases the gap between the moiré and the remote, conduction, and valence bands [77][78][79].But, because the AA area still remains a significant portion of the unit cell and, as we show, our main results are independent of twist angle, we do not expect lattice relaxation to strongly impact the effects of defects in TBG.Moreover, our implementation of the band structure, see below, achieves a finite gap isolating the moiré bands which is within experimentally measured bounds at the magic angle [33], leading to a quantitatively correct capturing of the pristine moiré bands.
Using the above stated position vectors, the hopping elements of Eq.( 1) can be explicitly calculated.For intralayer hopping we only include next neighbor processes with t ij equal to the graphene hopping t π .This is an approximation made in order to preserve the sparsity of the Hamiltonian matrix for computational efficiently.The result, compared to the full hopping model, is a rescaling of the twist angle, which we can simply compensate for to still achieve the magic angle, and an enhanced band gap isolating the moiré bands at the magic angle, which is beneficial when ignoring lattice relaxation as stated above.For the interlayer elements we use the Slater-Koster form [75,80] (3) where r ij = r i − r j − R is the displacement between the carbon sites i and j of unit cells connected by the lattice vector R of the moiré lattice.The sum over the lattice vectors R is important for sites near the edge of the unit cell, where the hopping occurs between different unit cells.The parameters a c = a/ √ 3 is the intralayer carbon to carbon distance.Other parameters are fixed according to the electronic structure of single and ABstacked bilayer graphene [75], with t π = 2.7 eV, out of plane hopping amplitude t σ = 0.48 eV and overlap decay length λ = 0.184a c , which was shown to reproduce well the pristine band structure of TBG [27].Because of the exponential form of the Slater-Koster terms, they become negligible for distant sites and we introduce a cutoff for d > 6a in order to preserve the sparsity of the Hamiltonian matrix, with no notable impact on the results.We also remark here that because there are no spin active terms in the problem we consider spinless fermions throughout the whole work.Hence, all observables can be thought of as of a single spin species and a multiplicative factor of 2 leads to the values for the physical, spinful electrons.
From Eq. ( 1) we obtain the energy spectrum and eigenvectors of the system by diagonalization of the sparse Hamiltonian matrix using the eigenvalue solver PRIMME [81], focusing on the relevant energy regions near the charge neutrality point (CNP).The CNP, or half-filling, of the pristine system lies at the same energy as the Dirac point.The low-energy band structure is shown in Fig. 1(b) at a twist angle θ ≈ 1.2 • , approximating the magic angle for our Hamiltonian.The four central lowenergy bands are the so-called moiré bands, separated from the valence and conduction bands by a finite band gap and also with completely flat regions in the Brillouin zone at the magic angle [31,32,74,82].Most of our calculations are done for θ ≈ 1.2 • , modeling the magic angle regime, including the important finite band gap and flat band regions.More precisely, for the magic angle calculations, we use the parameters (p, q) = (1, 55) in the commensuration condition cos(θ) = (3q 2 − p 2 )/(3q 2 + p 2 ) [74], such that the moiré unit cell contains 9076 atoms and the moiré length is L m ≈ 47.6a.However, as we later show, our findings hold for a much wider range of angles.We label the moiré bands from top to bottom as m 1 , m 2 , m 3 , and m 4 , and label the conduction and valence bands as c and v, respectively.The energy range spanned by the moiré bands, from E 0 to E 0 + W , where W is the moiré bandwidth, is from here on referred to as the moiré energy range, E m .The boundaries of this range are outlined by horizontal dashed purple lines in Fig. 1(b), and in all later band structure plots, in order to create an explicit reference to the pristine case.By integrating the local density of states (LDOS) over the moiré energy range we obtain the charge density of the moiré bands, where e is the electron charge, N k a normalization factor equal to the number of k-points sampled, and ψ n,k are the eigenstates with energies E n,k .In each case we choose a grid density for reciprocal space sampling such that we observe a convergence of the main features in the charge densities.We show the moiré charge density for pristine TBG at θ ≈ 1.2 • in Fig. 1(a), where we see clearly that the moiré bands are primarily localized in the AA regions.
With the pristine tight-binding model established above, we now introduce defects into the lattice.Specifically, the defects we consider are non-magnetic potential impurities and vacancies.In order to introduce potential impurities we define a perturbing potential v which enters the Hamiltonian as an onsite energy term on the affected sites For most of this work we focus on the simplest type of perturbing potential, which is that of a perfectly localized impurity with , where E I is the impurity strength and d is the impurity site.In Sec.III C we additionally consider v having a Gaussian profile around a central site, in order to investigate extended impurities.A vacancy can be introduced on the site i by letting v i → ∞, such that this site effectively decouples from the rest of the lattice.For numerical stability, however, we use an equivalent approach of simply restricting the sum of Eq. ( 1) such that no hopping is allowed into or out of the vacancy sites V, For a perturbing potential v or a set of vacancies V we obtain the energy spectrum of the total Hamiltonian H(k) = H V (k) + H I (k) near the CNP, where the moiré bands are located.Note that in the absence of vacancies, H V simply reduces to the pristine Hamiltonian H 0 .These perturbing terms are effective models of both actual vacancies and chemisorbed adatoms in the case of H V and physisorbed adatoms in the case of H I .By allowing the potential v i to have a finite spatial extent in H I and thereby creating extended impurities, we can even model larger physisorbed adatoms or even small molecules.
Because the creation and annihilation operators of Eq. ( 1) are of Bloch electrons, we effectively model a periodic lattice of defects, repeated in each moiré unit cell when solving H(k), even though the defect concentration is only 1 in 10 4 for a single defect per unit moiré unit cell.We are also interested in the case of completely isolated defects.In order to study these within an atomistic model we turn to the use of supercells.In this approach we enlarge the unit cell of our lattice, by considering a supercell consisting of an m × n array of moiré unit cells, with only a single defect per such supercell.The periodicity of the defect is then that of the supercell and thus by choosing large enough supercells we can completely isolate the Bloch copies of the defects from one another, enabling us to study the isolated defect limit.

III. RESULTS
In order to perform an analysis of the effects that defects have on the low-energy electronic structure of TBG we begin with the most simple type of defect: a singlesite impurity at site d, with the impurity potential given by v i = E I δ ,id .In Sections III A and III B we explore both the influence of the impurity strength, including the vacancy limit, and the defect location.We choose the defects to be always in the top layer since for single defects both layers are equivalent due to symmetry.Because of the sheer number of possible defect sites, we choose representative sites in the AA, AB, and DW regions as candidates for the defect location.For the more computationally intensive calculations of the charge density, we focus on the vacancy limit and use the band structure results to guide our interpretations.Then in Section III C we turn our attention to putative TPFs at the Dirac point created by defects, and also extend our study to multiple defects per unit cell as well as defects with extended spread.Finally, in Section III D we study the length scale behavior of the effects of defects.In Sections III A, III B, III C we stay approximately at the magic angle using θ ≈ 1.2, while in Section III D we study the behavior away from this regime.

A. Extreme electronic structure sensitivity to atomic size lattice defects
Starting with a single defect per moiré unit cell, we find that for most defect locations the overall effect of an impurity or vacancy on the band structure is similar.The main features are most easily seen for the case of a defect in the AA region, illustrated in Fig. 2. For an impurity strength up to E I = 0.1t, see Fig. 2(a), we find that the band structure changes very little even at the small energy scale of the moiré energy range.The most significant change seen is a breaking of the exact fourfold degener-acy of the Dirac point at K.However, as E I increases, the topmost of the four moiré bands, m 1 , detaches from the others, except at the Γ point, and lifts in energy.This band lifting is significant already at E I = 1t, see Fig. 2(b) where the m 1 band has already been removed almost entirely from the moiré energy range.This result shows how extremely sensitive the low-energy spectrum of TBG is with respect to defects, especially since the defect concentration here is only of the order of ∼ 0.01%.We note that this band removal behavior can be captured by band structure measurements, such as ARPES, or measurements of the density of states, such as STM or transport measurements, which would show a high peak from the relatively flat band m 1 at a much higher energy than the peaks due to the moiré bands in the pristine system.We also briefly note that m 1 never fully disconnects from the other moiré bands at Γ, while one of the valence bands v lifts in energy just enough to touch m 4 , closing the energy gap also from below and allowing m 3 to detach from m 4 at this point in the process.These degeneracies are discussed in more detail in Sec.III B.
By further increasing the impurity strength we reach around E I = 6t a behavior asymptotic in impurity strength, see Fig. 2(c) and Fig. 2(d) for a single vacancy.The resulting band structure contains three bands (m 2 , m 3 , m 4 ) within the moiré energy range, with the missing band (m 1 ) having acquired a parabolic character and joined the conduction bands.This results in a three-band moiré band structure that is not gapped neither from below or above.However, because m 1 and the highest lying valence band have both a very strong curvature, this band touching still results in a very small density of states when compared to the one from the leftover moiré bands.Beyond the moiré bands only containing three bands we also find that the middle band (originally m 3 ) becomes much flatter compared to pristine TBG and lies essentially at the CNP.This directly exemplifies the possibility to engineer new flat band structures in magicangle TBG by using impurities or vacancies.
In Fig. 3(a) we show the accompanying change in the moiré charge density, ∆n m , induced by the same vacancy as in Fig. 2(d).Each site in the unit cell is represented by colored dots, with the intensity of the color representing the magnitude of ∆n m .Blue (red) represents a decrease (increase) in charge density due to the vacancy, with the vacancy site encircled in magenta.Here the contributions from each layer are superimposed, with a finite transparency of the dots allowing for better visualization since the dots representing different sites often overlap, with the sites with greater magnitude in ∆n m brought to the foreground.We see that a vacancy in the AA region induces a strong depletion of the charge density in the AA region of the unit cell.This might at first seem counter-intuitive with respect to what is known about vacancies in graphene, where a vacancy is known to induce a zero-energy state centered on the vacancy, which causes a positive and localized change in charge density [35,37].Instead, the AA depletion observed here has to be un- derstood from the band structure of Fig. 2(d).Because one of the moiré bands, m 1 , has been removed from the moiré energy range, and because these bands are localized in the AA regions, see Fig. 1(a), the vacancy must cause a depletion of states over the entire AA region, not just locally around the vacancy.A similar effect occurs for impurities with strengths above 1t, since m 1 is already then removed from the moiré energy range, see Fig. 2(b).This extended effect of just single-site defects in TBG is in sharp contrast with how the same defect behaves in monolayer graphene, where the defect-induced state affects a much smaller region, and only locally around the defect.The connection between the depleted AA regions and the lifting of m 1 can be further verified by integrating ∆n m over the whole unit cell, where we find that there is 1e less charge in the integrated energy range, exactly corresponding to the removed moiré band.This overall charge depletion in the AA region easily overshadows the monolayer graphene defect state in Fig. 3(a), which we discuss in Sec.III B.

B. Role of defect location
The results in Sec.III A were for a particular choice of a defect site in the AA region.Now we also explore how the band structure and LDOS change for different defect locations.First of all, we find that different sites in each region within the same sublattice and in the same layer behave similarly to each other.Moreover, we find that the AA and DW regions show a sublattice symmetry with respect to the defect site, such that single defects in different sublattices lead to approximately the same band structure and charge density, assuming the sublattice components of the charge density are also switched (n m ) A ↔ (n m ) B .In all these cases we observe the removal of the m 1 band from the moiré energy range and the consequent depletion of the DOS in the AA regions, as well as a migration of the m 3 band towards the CNP, where it becomes very flat in the vacancy limit.We note, however, that a larger impurity strength is needed in order to observe band structure changes for defects outside the AA region, as expected since the moiré bands have much smaller presence outside the AA region.In the DW case we also observe that a valence band v becomes flatter and rises in energy and, for some defect sites, even partially enters the moiré energy range.In this case, the m 1 band detaches at Γ, leading to an energy gap above the moiré bands.
We next point out an observation we did not comment on in the last subsection.In the charge density for AA defects, we observe that the charge depletion of the AA region is reduced for sites immediately surrounding the defect that sit on the opposite sublattice to the defect site.This is a local effect, which we attribute to an induced, localized defect state akin to the one of monolayer graphene which appears on the sublattice opposite to the defect [35,37].In TBG, this state ends up partially canceling the change in charge density due to the removal of the m 1 band from the moiré energy range.Technically, for this defect location and at the magic angle, we cannot isolate the localized graphene-like defect behavior from the moiré depletion because the two effects are occurring in the same spatial location.However, this localized graphene-like defect state can be directly observed when the defect is in regions other than the AA region.We find that it possesses a C 3 symmetry and decays almost within a few atomic sites, just as in monolayer graphene [35,37].For a DW defect, in addition to the important AA region depletion and the graphene-like localized defect state, we also observe a slight depletion of states within the DW region itself, but it is contained to sites within the same sublattice as the defect site.For more details we refer to Appendix A.
Moving on to single-site defects in the AB region, we observe a behavior very similar to that of a DW defect for half of the sites of the AB region.More precisely, these are sites that have a neighboring site in the same position in the other layer.For this reason we call these sites higher-coordinated (HC) sites, while the other sites, which lie at the center of the carbon rings of the other layer, we refer to as the lower-coordinated (LC) sites.While a vacancy in an AB-HC site leads to a very clear charge depletion of the AA region, very similar to the behavior so far discussed (see Appendix A for details), the introduction of a vacancy in an AB-LC site leads to a drastically different charge density, as illustrated in Fig. 3(b).In this case we find a much smaller ∆n m in the AA regions.Instead, the most pronounced feature is the graphene-like localized defect state with its C 3 rotational symmetry.A calculation of the band structure in this case yields a moiré structure with four bands, which explains why the charge density is similar to the pristine case.
At first sight it might look as if TBG is thus insensitive to a vacancy on an AB-LC site, given the largely unmodified charge density in AA regions and the same number of bands, but this is not true.This precarious pitfall is revealed by the more careful analysis in Fig. 4, where we interpolate between the pristine case and the vacancy limit by tuning the impurity strength E I .In Fig. 4(a) we show that even at the large value of E I = 20t very little change from the pristine band structure is present among the moiré bands.The main difference is the attachment of a valence band v to the moiré bands at Γ, closing the band gap.By further increasing E I up to 40t, see Fig. 4(b), we find that this valence band v starts a process of inverting its curvature.In this process, v is initially outside of the moiré energy range, except at Γ, but as E I increases, the energy of v, especially around at the M point, is brought into the moiré energy range, see Fig. 4(c).Then finally, with further increasing E I , v becomes less dispersive and completely joins the moiré bands, see Figs. 4(d-e).At the same time we observe the detaching and removal of the m 1 band from the moiré energy range.Thus, the similarity of the vacancy band structure to the pristine band structure is deceiving and FIG. 4. Band replacement process induced by a single-site impurity in an AB-LC site.(a-e) Band structure for impurities with strengths 20t, 40t, 60t, 100t, and a vacancy, respectively.The vacancy limit is achieved for an impurity strength of around 500t.Although initial and final band structures and their corresponding density of states are very similar, they do not contain the same bands, but a band replacement process takes place for increasing EI , where the valence band v enters the moiré energy range as the m1 band is removed.Here θ ≈ 1.2 • .
the system has in fact undergone a band replacement process in between: one band joins from the valence bands, while one of the pristine moiré band is lost to the conduction band in the vacancy limit.We notice here that the impurity strength necessary to drive this band replacement process is quite large, with the vacancy limit only being reached at around E I = 500t, to be compared to defects in the AA region, where the vacancy limit is achieved already around E I = 6t.We further point out that the AB-LC sites correspond to half of the sites in the AB regions, which gives experimental relevance to these results.
Another interesting observation to be made about the role of the defect location is what happens to the degeneracies of the moiré bands at the Γ-point.In the pristine system the m 1 (m 3 ) and m 2 (m 4 ) bands are degenerate along the path from K to Γ (see Fig. 1).As defects are introduced, the degeneracy along this path is lifted as m 1 rises in energy.However, in the cases of AA and AB-LC defects shown in Figs. 2 and 4, the degeneracy remains at the Γ-point for all impurity strengths.This is because the wavefunctions of the moiré states at Γ have nodes at these atomic sites.These nodes must be present because these states each constitute non-trivial irreducible representations (IRREPs) of the C3 symmetries of the lattice with rotation axes at the center of the AA, AB or BA regions.These IRREPs have eigenvalues exp(±i2π/3) for the generator of C3 and show a phase winding around the rotation axis.This leaves the phase at the central site undetermined and hence the wavefunction amplitude must exactly vanish there.As a consequence, the moiré states at Γ are completely insensitive to perturbations added to these lattice sites, so the observed degeneracy is in these cases protected by the lattice symmetry.As the defect site is moved away from the center site of these regions, the Γ-point degeneracy is lifted, as illustrated by the cases of defects in DW and AB-HC sites (see Appendix A).We also note that in the AB-LC case of Fig. 4 the m 3 and m 4 bands separate at Γ, but only after the latter is joined by a valence band, such that a twofold degeneracy is always present, although after this process the degeneracy is between m 4 and v.
Finally, we also point out that pristine TBG has an approximate valley symmetry, such that its bands can be classified according to a valley quantum number.Atomic size defects, however, lead in general to intervalley mixing, such that valley is no longer a good quantum number.In Ref. [60] the authors have shown that defects in different positions couple to the valleys in different ways.Defects in the AB and BA regions preserve valley polarization, such that the TPF structure shows an intact Dirac cone from one valley and a split Dirac cone from the other valley.On the other hand, defects in the AA and DW regions lead to valley-unpolarized states, which can be understood as a consequence of intervalley mixing due to the defect.Due to these extensive intervalley processes, we refrain from further discussions regarding the valley quantum number.

C. Destroying triple point fermions
For all the defects we considered so far, a triple degeneracy of the moiré bands exists at the Dirac points K, K , see e.g.red arrows in Fig. 2(d).This gives rise to a TPF [52][53][54][55][56][57][58][59][60], since there is a Dirac spectrum crossed by a flat band, giving rise to a triply degenerate Dirac point.This is in contrast with the pristine TBG case, where, because of the valley degree of freedom, there are two degenerate Dirac cones at K, K , leading to a fourfold degeneracy of the Dirac point.Hence, the TPF scenario represents a reduction of the degeneracy of the Dirac point.The presence of TPFs in TBG has already been discussed in Ref. [60], where their robustness with respect to the impurity strength has been highlighted.Here we extend these results by showing that this robustness actually relies on the assumption of a single-site perturbing potential.In fact, we show that the TPF is split upon the introduction of either multiple defects or an extended impurity.
Let us begin by investigating the three-fold degeneracy in more detail, focusing on the spectrum near K, as a similar argument holds for the other inequivalent Brillouin zone corner K .In the pristine case there is an exact fourfold degeneracy at this point due to Dirac cones from the two so-called valleys of TBG, which originate from the two layers.Moreover, because the conduction and valence bands are strongly dispersive in comparison with the moiré bands, there is a sizable gap ∆E K in this region of reciprocal space, much larger than the energy splitting of the moiré bands E where k is a lattice wave vector near K and v F is the Fermi velocity corresponding to the slope of the Dirac cones of TBG.Because of this, it is possible to treat the effect of a weak perturbation on the moiré bands around K by projecting it into the moiré band subspace.Particularly at the Dirac point, k = K, because of its fourfold degeneracy, we need to use degenerate perturbation theory.
With the above considerations, the energies of the moiré bands can be approximated by projecting the perturbed single particle Hamiltonian h = h 0 + h I into the moiré subspace, where h 0 and h I are the single particle versions of Eqs. ( 1) and ( 5) evaluated at K. Letting P m ≡ i |m i m i | be the projector into this subspace, where |m i are the degenerate eigenstates of the moiré bands at K, the projected Hamiltonian becomes where we use the fact that the pristine moiré bands are fourfold degenerate at the Dirac point, with energy E K and further define hI ≡ P m h I P m .We next recall that a rank-1 operator O has only one non-zero eigenvalue and can be written as an outer product O = | |.Using this, we note that for a single-site impurity at site x 0 , h I = E I |x 0 x 0 |, both h I and hI are rank-1 operators.The proof for h I follows directly from its definition.For hI we have where ψ mi (x) = x|ψ mi is the wavefunction of the moiré band m i at K and in the last line we simply defined Thus the operator hI has a single non-zero eigenvalue E I and is a rank-1 operator.As a consequence h has three degenerate eigenvalues E K and another eigenvalue E K + E I .This means that a single-site impurity reduces the degeneracy at the Dirac point from fourfold to threefold, leading to the formation of a TPF, exactly as earlier predicted [60].
Next, let us consider two distinct, but still weak and perfectly localized impurities, h I and h I .What we are interested in is whether the TPF found above remains for the total perturbation h (t) I .In order to do this, we again investigate the rank of the projected total perturbation, h(t) I P m .We can quickly verify that the unprojected perturbation h (t) I must be of a higher rank.One way to do this is to tentatively assume that h(t) I is rank-1.This means that it can be expanded as an outer product.By definition, this requires that there exist α i such that where h ij = m i |h|m j are the matrix elements of an operator h in the basis of the moiré bands |m j and in the last line we used the definition of an outer product.
For this last equality to hold h(t) I must be separable, that is, there must exist α i such that This constraint is however too restrictive and is, in the general case, not satisfied by any α i .Thus, in general, h I cannot be rank-1 and thus the introduction of a second weak impurity in the unit cell leads to a further change in the degree of degeneracy at the Dirac point and a subsequent splitting of the TPF.Moreover, although the argument above does not hold for arbitrarily large impurity strengths E I , since we assumed E I ∆E K at the start, we verify through extensive numerical calculations that the conclusion that h (t) I is not rank-1 and thus that multiple defects further reduce the degeneracy of the Dirac point still holds even in the vacancy limit.As an example, in Fig. 5(a) we show the band structure for the simple case of two nearby vacancies in the AA region at sites of opposite sublattice in the same layer.In this case we observe the removal of another moiré band from the moiré energy range as compared to the single vacancy case.This clearly splits the TPF by lifting the needed degeneracy, as there are now only two bands degenerate at K. We verify that this splitting occurs for many other defect locations, including for vacancies in different parts of the unit cell, such as one in the AA region and the other in the DW region.The only exception we find so far, where the TPF survives the inclusion of more than one single-site impurity, is tied to the lower-coordinated sites in the AB region, AB-LC, where the band replacement process discussed in see Sec.III B can in fact restore the TPF.
Fig. 5(a) also illustrates another interesting feature, namely that the effect of introducing two vacancies in the AA region is an additive process of the effects of the single vacancies, in the sense that each vacancy is responsible for removing a single band from the moiré energy range and thus with two vacancies, only two moiré bands are left.In exploring the possible combinations of two vacancies, we find that this is a common pattern.However, we find that it is not quite universal.For example, the combination of two vacancies in the DW region results in four bands in the moiré energy range, as shown in Fig. 5(b).The same result is also obtained when having two vacancies in the AB region or one in each type of the HC/LC site (not shown).Still, in all of these cases, the TPF is split.This non additivity of the effects of defects leads us to investigate the range of influence of each defect later in Sec.III D.
In order to further corroborate that the origin of TPFs in TBG is tied to a rank-1 perturbation, we next consider the case of a single but extended impurity.We model this by giving the impurity potential a Gaussian profile with a spread of σ centered around the impurity site.We here set the strength of the Gaussian profile such that the perturbing potential is 1t in the central site.This can be regarded as a simple model for an impurity that affects multiple sites around its binding center, realistic for molecule adsorbates or for an adatom absorbed in the honeycomb lattice hollow site.For simplicity, we cut off the Gaussian at a distance of 2.5a away from the binding center, where a is the graphene lattice constant, which we verify does not influence the results.In Fig. 5(c) we show the band structure resulting from such Gaussian impurity perturbation with the impurity center located in the DW region, using σ = 0.8a.We see that four bands remain mostly in the moiré energy range but notably the triple band crossing at K is clearly no longer present, meaning the TPF is split.We also verify that the splitting of the TPF holds for other impurity locations.In fact, when the extended impurity is centered in the AA region we find that all degeneracies of the moiré bands at K are lifted, see Fig. 5(d).We also observe in this case a lifting of all four moiré bands, here illustrated for σ = 0.5a.For a larger spread σ = 0.8a we find an even larger depletion as the moiré bands, which are lifted further in energy such that they leave behind only strongly dispersing states near Γ (not shown).This is a fascinating result showing how the entire moiré band structure can be destroyed by a single weak extended impurity.Taken together, our results in this subsection show that if the defects are not simple rank-1 perturbations, but instead, for example, multiple defects or extended impurities, then the TPF at K and K is generally split.Based on these results, we do not expect TPF to be likely observed in TBG.

D. Isolated defects and length scales
In the previous subsections we explored the effects of a defect lattice with same periodicity as the pristine system, a model which can be implemented in practice by defect engineering.We now turn to the limit of an isolated defect.In the fully atomistic framework we use, this amounts to reducing the periodicity of the defect lattice, such that the distance between the different Bloch copies of the same defect become large enough as to not influence each other.We explore different supercell sizes and defect locations, with the constraint that all defects are within a single moiré unit cell inside the supercell, which we refer to as the defective moiré unit cell.In all defect configurations we explore, we find that an asymptotic behavior is reached in a specific spatial direction when the two defects are separated by three or more moiré unit cells along that particular direction.This means that for 3 × 3 or larger supercells single defects are effectively isolated from each other.
We illustrate the result for a single isolated vacancy in Fig. 6, where we plot the change in the moiré charge density ∆n m for a 3 × 3 supercell with a vacancy in a DW region.As a guide to the eye we mark the boundaries between the individual moiré unit cells outlined by green lines.In the defective moiré unit cell (bottom left) we observe a similar charge density redistribution as in the 1 × 1 supercell case discussed in Sec.III B, with a clear depletion in the AA regions closest to the vacancy.This can again be understood as one of the moiré bands associated with this unit cell leaving the moiré energy range.However, for a 3×3 supercell, the band folding caused by the reduced periodicity results in (3 × 3) × 4 = 36 bands in the moiré energy range in the pristine case.Thus, out of these 36 bands, only one exits this energy range when the vacancy is introduced, causing a depletion of 1e distributed through the AA region of the defective moiré unit cell.Away from the defective moiré unit cell we observe that the charge density quickly recovers to its pristine value, even in the AA regions.We verify that a similar behavior is present also for larger supercells, which shows that we use a system size capable of modeling the asymptotic isolated defect limit.
We find that the separation required between vacancies for reaching the isolated defect limit is reduced when the vacancy is in the AA region.In this case, we find that even supercells as small as 2 × 2 give results converged to the isolated defect limit.We illustrate this in Fig. 7(a,b), where we plot the contribution to ∆n m on the sites of the top layer in sublattice A and B, respectively, with the vacancy being on an A sublattice site.Note how only the AA region with the vacancy has an altered charge density and thus the isolated defect limit is achieved already for 2 × 2 supercells.Moreover, it is clear that the charge density depletion is primarily in the same sublattice as the vacancy.For the vacancy-free layer we find a similar but smaller change in charge density in the AA regions.Taken together, Figs. 6 and 7(a,b) illustrate that a single vacancy generally influences the moiré pattern up to a distance of the moiré length L m .
In order to further corroborate that the vacancies affect a region with radius of the order of L m , we vary this length by changing the twist angle away from the magic angle, up to θ ≈ 6 • .We then study ∆n m along a line cut passing through a vacancy using a 3 × 3 supercell, depicted by the yellow lines in Fig. 7(a,b).This way, the vacancy concentration per unit area changes with twist angle, but the vacancy concentration per unit cell stays the same and can easily be compared.Also, the vacancies always stay isolated from each other and thus we stay within the isolated defect limit.We particularly choose a line cut direction that goes through a nearest neighbor site of the vacancy in the vacancy layer in order to also probe the graphene-like localized defect state that exists in this layer.The depletion of the AA regions is however visible for any line cut direction.We further choose a vacancy location in the center-most site of the AA region because this site is always present for all twist angles.This site is in sublattice A, while the layer index does not matter.
In Fig. 7(c) we plot change in charge density for the A sublattice sites along the line cut, where we know that graphene-like localized defect state does not contribute because it is primarily located in the opposite sublattice [35,37].Here we scale the x-axis with respect to the moiré length L m , which is dependent on the varying twist angle.We observe a clear trend where the depletion in the AA region recovers away from the vacancy with a length scale that approaches L m , but is smaller especially for smaller angles approaching the magic angle.In contrast, in Fig. 7(d) we show ∆n m for the B sublattice sites along the line cut.Here we let the x-axis be normalized with the atomic scale a, corresponding to the graphene lattice constant, as we find that the dominant contribution comes from the graphene-like localized defect state, giving a positive change in charge density on the atomic scale.This graphene-like localized defect state decays over a similar length scale, set by a, for all twist angles, further corroborating that its origin is due to graphene physics and not the moiré pattern.We thus conclude that isolated defects affect TBG on two different length scales.On the atomic length scale, a, it induces a localized defect state similar to that of in monolayer graphene.On the angle-dependent moiré length scale, L m , it induces a strong charge depletion in the AA regions, which near the magic angle can be understood from the removal of an entire moiré band from the low-energy region.We note that at larger angles, the moiré bands are not energetically separated from the conduction and valence bands, which means that the simple picture that one of the moiré bands leaves the moiré energy range breaks down.However, we verify the presence of a flat defect-induced band at the energy of the Dirac point in the vacancy limit, similar to the magic angle case, for all angles up to 6 • .The one exception we find to the latter behavior is for defects in AB-LC sites, since in this case there is a moiré band replacement instead of a band removal.

IV. CONCLUDING REMARKS
To summarize, in this work we show that the lowenergy band structure of twisted bilayer graphene (TBG) is extremely sensitive to atomic size lattice defects even at very low concentrations.In particular, we show that a single weak non-magnetic impurity in each AA region is able to cause a large depletion of charge in the lowenergy regime and in the whole AA region of the order of 1e per spin species due to the lifting of one of the low-energy moiré bands into the conduction bands.We investigate different defect locations and verify that this behavior is quite general and thus illustrates a direct way to manipulate the low-energy moiré band structure using impurities.The only notable exception we find is for a special set of defect sites in the AB region, where a band replacement process happens instead, where one moiré band is lifted to the conduction band, while another joins from the valence band, resulting in a reconstruction of the original low-energy band structure in the vacancy limit.We strongly suspect that this band replacement directly influences the topology of the moiré band structure although that remains to be verified.
We further find that the previously reported defectinduced TPFs in TBG [60], which represent a triple degeneracy at the Dirac point, rely on the rank-1 perturbation characteristic of single-site defects, and is thus not generally present.In fact, we show that the introduction of multiple defects or more realistic extended impurities easily split the TPF by lifting the degeneracy at the Dirac point, in some cases even completely removing the Dirac point.Finally, we use supercells to reach the isolated defect limit, where we find the previous results to hold locally in the moiré unit cells surrounding the defect.By varying the twist angle we are further able to identify two length scales, with the atomic scale displaying a graphene-like localized defect state, and the twist angle-dependent moiré length controlling the charge depletion of the AA region and accompanied moiré band restructuring.
Our results establish how non-magnetic impurities and vacancies drastically change the band structure and charge density of TBG at and near the magic angle, which can be experimentally verified with ARPES, STM, or transport measurements.Incorporating these profound changes of the moiré bands will further be important for analyses of quasiparticle interference patterns.These measurements should be performed at low enough temperatures for good energy resolution (we estimate around ∼ 10K), but above the critical temperatures of any emerging electronic orders of TBG.We expect that the changes induced by defects will also have profound impact on these electronic orders, including the superconducting and correlated insulator orders, since these orders depend not only on the interactions, but also heavily on the underlying normal-state band structure.This impact will be particularly large on any mechanism relying on the symmetries of the system or number of moiré bands, as defects strongly modify the low-energy moiré band structure through band lifting and band replacement processes.In fact, we show that the number of moiré bands can easily change from the pristine case of four to three or two, or even be completely annihilated, with only strongly dispersive bands left in the low-energy region.This sensitivity of TBG with respect to impurities and vacancies demonstrate the need to understand the disorder level before further analyzing any electronic ordered state.It also opens up the possibility of using defects to engineer the low-energy electronic structure of TBG in order to produce a desired number of flat bands and thereby possibly other electronic orders.
In the main text we discuss the effects of defects in AA and AB-LC sites extensively, while focusing less on defects in the DW and AB-HC.The reason for this is that these latter sites show a behavior quite similar to that of AA sites.Still, for completeness, we here provide relevant data for defects in the DW and AB-HC sites, supporting the conclusions in the main text.In Fig. 8(a,b) we show the change in moiré charge density ∆n m induced by a vacancy in DW and AB-HC sites, respectively, using a 1 × 1 supercell at θ ≈ 1.2 • and with the vacancy highlighted by a pink circle.Most importantly, and as mentioned in the main text, we find that the introduction of a defect induces a strong depletion of the AA region even in these cases where the defect is located far away from it.We additionally find a slight depletion of states in an extended region between the defect and the nearest AA regions.This provides further evidence of the effect that atomic size defects have on the moiré scale.Once again, the depletion of the AA regions can be understood from the evolution of the band structure as we interpolate between the pristine and vacancy limits with an impurity of finite strength, as we illustrate in Figs. 9 and 10 for DW and AB-HC impurities, respectively.In both cases we observe a band removal process akin to the the one discussed in the main text, with the m 1 band leaving the moiré energy range and joining the conduction bands.The impurity energies required to trigger the removal process is higher than for an AA impurity, with over 10t being required to remove the m 1 band from the moiré energy range.Just as in the AA defect case, this explains the charge depletion of the AA region, since it is there the moiré bands are located.We also note that the degeneracy of the moiré bands at the Γ-point, discussed in Sec.III B, is lifted, since the defect site in these cases is not on the rotation axis of a C3 symmetry.

FIG. 1 .
FIG. 1. Charge density and band structure of pristine magicangle TBG.(a) Charge density of the moiré bands of TBG, nm(x), with contributions from both layers shown.The AA, AB, and DW regions of the moiré unit cell are schematically highlighted and the moiré length Lm is marked.(b) Lowenergy band structure of TBG, along high symmetry points: K, Γ, M .Horizontal dashed purple lines mark the moiré bandwidth, W . Moiré, conduction, and valence bands are labeled by mi, c and v, respectively.Here θ ≈ 1.2 • and we use a 12 × 12 k-point grid to compute (a).

FIG. 2 .
FIG.2.Effects on the low-energy spectrum of magic-angle TBG from a defect in the AA region.(a-d) Band structure for potential impurities with strengths of 0.1t, 1t, 6t, and a vacancy, respectively.Around 1t the moiré band m1 leaves the moiré energy range, while the vacancy limit for the potential impurity is achieved around 6t. Red arrows in (d) highlight the TPF.Here θ ≈ 1.2 • .

FIG. 3 .
FIG. 3. Change in moiré charge density ∆nm from pristine magic-angle TBG due to a vacancy.(a) Vacancy in the AA region, showing a strong charge depletion of the AA region, due to the removal of one of the moiré bands.(b) Vacancy in an AB-LC site (defined in main text), where the main change in charge density is instead the graphene-like localized defect state.Pink circles mark the vacancy site and k-points were sampled in a 12 × 12 grid.In each panel the contributions from both layers are shown.Here θ ≈ 1.2 • .

FIG. 5 .
FIG. 5. Removal of the TPF of TBG.(a-b) Two vacancies per unit cell in the AA (a) and DW (b) regions.(c-d) Extended Gaussian impurity centered in the DW (c) and AA (d) regions.Gaussian spread of σ = 0.8a (c) and σ = 0.5a (d) with a cutoff of 2.5a.In all cases (a-d) the TPF is split.Here θ ≈ 1.2 • .

FIG. 6 .
FIG. 6. Change in moiré charge density ∆nm for a single DW vacancy in a 3 × 3 supercell.Depletion of the AA regions neighboring the vacancy.Pink circle marks the vacancy location and green lines mark the individual moiré unit cells.The contributions from both layers are shown.Here 12 × 12 k-point sampling was used and θ ≈ 1.2 • .

FIG. 7 .
FIG. 7. Two length scales associated with defects in TBG.(a-b) Change in moiré charge density ∆nm for a single AA vacancy in the top layer in sublattice A at θ ≈ 1.2 • , i.e. near the magic angle, for the A (a) and B (b) sublattice of the top layer.Green lines mark the individual moiré unit cells.(c-d) Change in moiré charge density ∆nm along the line cut in (a,b) marked by the yellow line, for different twist angles, for the A (c) and B (d) sublattice of the top layer.The x-axis is scaled with respect to the moiré length Lm (c) and atomic length a (d), respectively.Here 6 × 6 k-point sampling was used.

FIG. 8 .
FIG. 8. Change in moiré charge density ∆nm from pristine magic-angle TBG due to a vacancy in the DW region (a) and in an AB-HC site (defined in main text) (b), highlighting the strong charge depletion of the AA region due to the removal of one of the moiré bands, and also a smaller depletion in the larger vicinity of the defect.Pink circles mark the vacancy site and k-points were sampled in a 12×12 grid.In each panel the contributions from both layers are shown.Here θ ≈ 1.2 • .

FIG. 9 .
FIG. 9. Band removal process for a defect in the DW region.(a-d) Band structure for potential impurities with strengths of 1t, 10t, 100t, and a vacancy, respectively.Here θ ≈ 1.2 • .