Intrinsically Undamped Plasmon Modes in Narrow Electron Bands

TL;DR

This paper predicts that in narrow-band 2D electron systems like moiré graphene, plasmon modes can be intrinsically undamped due to strong coupling effects, leading to enhanced optical coherence and observable interference patterns.

Contribution

It introduces the concept of intrinsically undamped plasmon modes in narrow electron bands, highlighting their emergence when plasmon dispersion exceeds the particle-hole continuum.

Findings

Landau damping is quenched in narrow-band systems with large coupling constants.

Undamped plasmon modes extend into the energy gap, avoiding electron-hole pair excitation.

Enhanced optical coherence and interference signatures are predicted in moiré graphene at magic angles.

Abstract

Surface plasmons in 2-dimensional electron systems with narrow Bloch bands feature an interesting regime in which Landau damping (dissipation via electron-hole pair excitation) is completely quenched. This surprising behavior is made possible by strong coupling in narrow-band systems characterized by large values of the "fine structure" constant $\alpha=e^2/\hbar \kappa v_{\rm F}$. Dissipation quenching occurs when dispersing plasmon modes rise above the particle-hole continuum, extending into the forbidden energy gap that is free from particle-hole excitations. The effect is predicted to be prominent in moir\'e graphene, where at magic twist-angle values, flat bands feature $\alpha\gg1$. The extinction of Landau damping enhances spatial optical coherence. Speckle-like interference, arising in the presence of disorder scattering, can serve as a telltale signature of undamped plasmons directly accessible in near-field imaging experiments.

Equations

$${\omega }_{\mathrm{p}}(q)=\sqrt{{\beta }_{q}q},$$

$${\omega }_{\mathrm{p}}\sim \sqrt{\alpha }W\gg W,\alpha \sim 20-30,$$

$$S(\mathbf{r})=J_0\int d^2{\mathbf{r}}^{\prime }{G}_E(\mathbf{r}-{\mathbf{r}}^{\prime })\eta ({\mathbf{r}}^{\prime })G_E({\mathbf{r}}^{\prime }-\mathbf{r}),$$

$$|S_{\mathbf{k}}{|}^2\sim \frac{|{\eta }_{\mathbf{k}}{|}^2}{|k^2-4(q_0-i\delta {)}^2|},$$

$$H_{\mathrm{toy}}=\left(\begin{array}{cc}0& h_{\mathbf{k}}\\ h_{\mathbf{k}}^{\ast }& 0\end{array}\right),h_{\mathbf{k}}=\frac{W}{3}\sum _{{\mathbf{e}}_j}{e}^{i\mathbf{k}\cdot {\mathbf{e}}_j},$$

$$E_{s,\mathbf{k}}=s|h_{\mathbf{k}}|,{\Psi }_{s,\mathbf{k}}=\frac{1}{\sqrt{2}}\left(\begin{array}{c}s{e}^{i{\phi }_{\mathbf{k}}}\\ 1\end{array}\right),$$

$$\epsilon (\omega ,\mathbf{q})=1-V_{\mathbf{q}}\Pi (\omega ,\mathbf{q}).$$

$$\Pi (\omega ,\mathbf{q})=4\sum _{\mathbf{k},s,s^{\prime }}\frac{(f_{s,\mathbf{k}+\mathbf{q}}-f_{s^{\prime },\mathbf{k}})F_{\mathbf{k}+\mathbf{q},\mathbf{k}}^{s{s}^{\prime }}}{E_{s,\mathbf{k}+\mathbf{q}}-E_{s^{\prime },\mathbf{k}}-\omega -i0}.$$

$$F_{\mathbf{k}+\mathbf{q},\mathbf{k}}^{s{s}^{\prime }}=|⟨{\Psi }_{s,\mathbf{k}+\mathbf{q}}|{\Psi }_{s^{\prime },\mathbf{k}}⟩{|}^2=\frac{1+s{s}^{\prime }\cos ({\phi }_{\mathbf{k}+\mathbf{q}}-{\phi }_{\mathbf{k}})}{2},$$

$$\Pi (\omega ,\mathbf{q})=8\sum _{\mathbf{k},s,s^{\prime }}{f}_{s^{\prime },\mathbf{k}}\frac{F_{\mathbf{k},\mathbf{k}+\mathbf{q}}^{s{s}^{\prime }}(E_{s^{\prime },\mathbf{k}}-E_{s,\mathbf{k}+\mathbf{q}})}{(E_{s,\mathbf{k}+\mathbf{q}}-E_{s^{\prime },\mathbf{k}}{)}^2-(\omega +i{0}^{+}{)}^2}.$$

$$F_{\mathbf{k}+\mathbf{q},\mathbf{k}}^{s=s^{\prime }}\approx 1,F_{\mathbf{k}+\mathbf{q},\mathbf{k}}^{s=-s^{\prime }}\approx \frac{1}{4}{\left(\mathbf{q}\cdot {\nabla }_{\mathbf{k}}{\phi }_{\mathbf{k}}\right)}^2$$

$${\Pi }_1(\omega ,\mathbf{q})\approx \frac{4}{{\omega }^2}\sum _{\mathbf{k},s}{f}_{s,\mathbf{k}}(\mathbf{q}\cdot {\nabla }_{\mathbf{k}}{)}^2{E}_{s,\mathbf{k}}$$

$${\Pi }_2(\omega ,\mathbf{q})\approx 4\sum _{\mathbf{k},s}{f}_{s,\mathbf{k}}\frac{E_{s,\mathbf{k}}{\left(\mathbf{q}\cdot {\nabla }_{\mathbf{k}}{\phi }_{\mathbf{k}}\right)}^2}{4E_{s,\mathbf{k}}^2-(\omega +i0{)}^2}.$$

$${\beta }_q={\beta }_0+{\beta }_{1}q+O(q^2)$$

$$F_{\mathbf{k}+\mathbf{q},\mathbf{k}}^{nm}={\left|{\int }_{\Omega }{d}^{2}r{\Psi }_{n,\mathbf{k}+\mathbf{q}}^{†}(\mathbf{r})e^{i\mathbf{q}\cdot \mathbf{r}}{\Psi }_{m,\mathbf{k}}(\mathbf{r})\right|}^2,$$

$$\epsilon (\omega ,\mathbf{q})\approx 1+A(\mathbf{q})-\frac{B(\mathbf{q})}{{\omega }^2},$$

$$A(\mathbf{q})=\frac{8\pi e^2}{\kappa q}\text{sideset}{}^{\prime }{\sum }_{\mathbf{k},n,m}{f}_{m,\mathbf{k}}\frac{F_{\mathbf{k}+\mathbf{q},\mathbf{k}}^{nm}}{E_{n,\mathbf{k}+\mathbf{q}}-E_{m,\mathbf{k}}}$$

$$B(\mathbf{q})=\frac{8\pi e^2}{\kappa q}\text{sideset}{}^{\prime \prime }{\sum }_{\mathbf{k},n,m}{f}_{m,\mathbf{k}}{F}_{\mathbf{k}+\mathbf{q},\mathbf{k}}^{nm}(E_{n,\mathbf{k}+\mathbf{q}}-E_{m,\mathbf{k}}).$$

$${\omega }_{\mathrm{p}}^2\approx \frac{B(\mathbf{q})}{1+A(\mathbf{q})},$$

$$S(\mathbf{r})=J_0\int d^2{\mathbf{r}}^{\prime }{G}_E(\mathbf{r}-{\mathbf{r}}^{\prime })\eta ({\mathbf{r}}^{\prime })G_E({\mathbf{r}}^{\prime }-\mathbf{r}).$$

$$G_E(\mathbf{r})\approx \frac{e^{i{q}_0|\mathbf{r}|}}{\sqrt{2\pi |\mathbf{r}|}}{e}^{-\delta |\mathbf{r}|},$$

$$S(\mathbf{r})=J_0\int d^2{\mathbf{r}}^{\prime }\eta ({\mathbf{r}}^{\prime })\frac{e^{i2q_0|\mathbf{r}-{\mathbf{r}}^{\prime }|}{e}^{-2\delta |\mathbf{r}-{\mathbf{r}}^{\prime }|}}{2\pi |\mathbf{r}-{\mathbf{r}}^{\prime }|}.$$

$$S_{\mathbf{k}}=\int d^2\mathbf{r}{e}^{-i\mathbf{k}\cdot \mathbf{r}}S(\mathbf{r}).$$

$$S_{\mathbf{k}}=\int d^2{\mathbf{r}}^{\prime }\eta ({\mathbf{r}}^{\prime })e^{-i\mathbf{k}\cdot {\mathbf{r}}^{\prime }}\int d^2\mathbf{r}\frac{e^{-i\mathbf{k}\cdot \mathbf{r}+i2q_0|\mathbf{r}|}{e}^{-2\delta |\mathbf{r}|}}{2\pi |\mathbf{r}|},$$

$${\int }_0^{2\pi }d\theta \frac{1}{a+b\cos \theta }=\frac{2\pi }{\sqrt{a^2-b^2}}.$$

$$F_{\mathbf{k}+\mathbf{q},\mathbf{k}}^{s=-s^{\prime }}\approx \frac{1}{4}{\left(\mathbf{q}\cdot {\nabla }_{\mathbf{k}}{\phi }_{\mathbf{k}}\right)}^2\approx \frac{1}{4}\frac{q^2}{k^2}{\sin }^2\theta ,$$

$${\Pi }_2(\omega ,\mathbf{q})=-\frac{8q^2}{{\omega }^2}\sum _{\mathbf{k},s}{f}_{s,\mathbf{k}}s{v}_{F}k\frac{{\sin }^2\theta }{k^2}.$$

$${\Pi }_2(\omega ,\mathbf{q})\approx -\frac{2E_F}{\pi }\frac{q^2}{{\omega }^2}+\frac{2W}{\pi }\frac{q^2}{{\omega }^2}=\frac{2}{\pi }(W-E_F)\frac{q^2}{{\omega }^2},$$

$${\lambda }_2=2(W-E_F)/\pi ,$$

$${\Pi }_1(\omega ,\mathbf{q})$$

Full text

Landau damping, a process by which collective mode decays into electron-hole pairs, is often taken to be an integral attribute of graphene plasmon excitations wunsch ; hwang ; jablan ; koppens2012 ; basov2012 . Here, we predict extinction of this dissipation mechanism in materials with narrow electron bands, such as twisted bilayer graphene (TBG) mele2010 ; San-Jose2012 ; bistritzer2011 ; santos2007 ; fang2016 . Intrinsically undamped plasmons in narrow-band materials arise due to large fine structure parameter values $\alpha=e^{2}/\hbar\kappa v_{\rm F}$: strong interactions push plasmon dispersion into the energy gap above the particle-hole (p-h) continuum as illustrated in Fig. 1. In this region, plasmons become decoupled from p-h pair excitations. Dissipation quenching, which is a surprising manifestation of strong coupling physics, is a robust effect that persists up to room temperature and is insensitive to disorder (Figs. 1 and 2). Collective charge modes, which are damping free, are of keen interest for quantum information science as a vehicle to realize dissipationless photon-matter coupling, high-Q resonators, single-photon phase shifters and other missing components for photon-based quantum information processing toolbox gullans2013 . Although extinction of Landau damping is a general effect present in all narrow electron bands, our analysis will focus on TBG flat bands, a system of high current interest cao2018_1 ; cao2018_2 ; yankowitz2019 ; stauber2016 ; hu2017 , in which undamped plasmons can be directly probed.

Fig. 1 depicts plasmon mode for a narrow-band model that mimics the key features of the TBG band. Mode dispersion (red line) and its damping are of a conventional form at energies less than the bandwidth, $\omega\lesssim W$. At lowest energies, plasmon mode is positioned outside the p-h continuum, as expected; this suppresses the $T=0$ Landau damping, but does not protect the mode from decaying into p-h excitations through disorder scattering or from the conventional $T>0$ Landau damping hwang ; wunsch ; emani2012 ; low2014 ; principi2013 ; polini2008 ; yan2013 . At higher energies, $\omega\sim 2E_{\rm F}$ (marked by arrows in Fig. 1), the mode plunges into p-h continuum and is Landau-damped at $2E_{\rm F}\lesssim\omega\lesssim 2W$, even at $T=0$. However, an interesting change occurs after the mode rises above the p-h continuum. In the forbidden gap region, $\omega>2W$, it becomes damping-free, since at these energies there are no free p-h pairs into which plasmon could decay. This behavior is manifest in the $T$ dependence of the resonances, which are washed out with increasing temperature at $\omega\lesssim W$ but remain sharp at $\omega>W$ even at $T\sim E_{\rm F}$ (Fig. 1 b and c).

As we will see, mode dispersion has a square root form characteristic of 2-dimensional (2D) plasmons stern1967 ; principi2009 ,

[TABLE]

with a weak $q$ dependence in $\beta_{q}$ [Eq. (14)]. This expression, however, is valid not just at low energies, $0<\omega\lesssim W$, but also at higher energies, $\omega\gg W$, where the mode is undamped. While the dispersion in Eq. (1) is of the conventional 2D plasmon form, we emphasize that here it takes on a different role, as it describes the plasmon mode at frequencies much higher than the carrier bandwidth, extending to

[TABLE]

where the high-$\alpha$ values correspond to flat bands in magic-angle moiré graphene. Also, unlike the conventional plasmons, the dispersion in Eq. (1) is not limited to longest wavelengths. Indeed, as illustrated Fig. 1a, it extends to fairly high wavenumbers on the order of the mini Brillouin zone size.

The wavelengths of these plasmons are only 2 to 3 times greater than the moiré superlattice period. Such short wavelengths are of considerable interest for plasmonics and are within resolution of the state-of-the-art scanning near-field microscopy techniques basov2012 ; koppens2012 (currently as good as $10$ nm cohen2014 ; duan2012 ). In addition to measuring plasmon dispersion, these techniques can be used to directly visualize the qualitative change in the damping character and strength. Enhanced optical coherence will manifest itself in striking speckle-like interference as illustrated in Figs. 1d and 2.

Indeed, because of the absence of Landau damping at the energies of interest, $\omega>W$, and also because these energies are smaller than carbon optical phonon energies, the dominant dissipation mechanism is likely to be elastic scattering by disorder. At low energies, where plasmon mode coexists with p-h continuum, disorder scattering merely assists Landau damping, allowing plasmons to decay into p-h pairs by passing some of their momentum to the lattice. However, at the energies above p-h continuum, $\omega>W$, since the decay into pairs is quenched, disorder will lead to predominantly elastic scattering among plasmon excitations. Such scattering preserves optical coherence and is expected to produce speckle patterns in spatial near-field images as illustrated in Fig. 1d.

To model this behavior we consider the signal $S({\boldsymbol{r}})$, excited by the scanning tip and measured at the same location. Monochromatic plasmon excitation at energy $E$ is scattered by impurities or defects and on returning to the tip, produces signal

[TABLE]

where $\eta({\boldsymbol{r}})$ is the disorder potential, $J_{0}$ is excitation amplitude, and $G_{E}({\boldsymbol{r}})$ is the Green’s function of the plasmon excitation (see supplemental information). The spatial signal (Fig. 1d) exhibits a characteristic speckle pattern familiar from laser physics. In graphene plasmonics, speckle-like interference provides a direct manifestation of optical coherence enhancement in the absence of Landau damping. Accordingly, the Fourier transform of the image, $S_{{\boldsymbol{k}}}=\int d^{2}rS({\boldsymbol{r}})e^{-i{\boldsymbol{k}}{\boldsymbol{r}}}$, yields power spectrum $|S_{{\boldsymbol{k}}}|^{2}$ that features a ring-like structure; the ring radius is $k=2q_{0}$, where $q_{0}$ is the plasmon excitation wavenumber (Fig. 1d inset). Simple calculation, described in supplemental information, predicts power spectrum that sharply peaks at the ring:

[TABLE]

where $\delta$ is a parameter characterizing extrinsic damping due to phonon scattering and other inelastic processes. In the fully coherent regime ($\delta=0$) the quantity $|S_{\boldsymbol{k}}|^{2}$ exhibits a power law singularity at the ring, $k=2q_{0}$. As the amount of incoherent scattering increases, the peak is gradually washed out. This behavior is illustrated in Fig. 2.

We note that recent work stauber2016 analyzed interband plasmon excitations in TBG, which are dominated by polarization of the bands above the flat band and are distinct from the flat-band plasmons analyzed here. Recent experiment hu2017 reported observation of plasmons in TBG; however, their appeal for constructing intrinsically protected collective modes remained unnoticed in graphene literature. Also, plasmons in narrow bands were analyzed in the context of high-$T_{c}$ superconductivity pashitskii1993 , finding that plasmon mode can rise above the flat band. However, in cuprates, unlike moiré graphene, the narrow band is not separated from higher bands by a forbidden energy gap, and thus, the mode studied in ref. pashitskii1993 will plunge into a higher band before acquiring an undamped character.

Next, we present analysis of the hexagonal-lattice toy model that mimics the key features of Landau-damped and intrinsically undamped modes in TBG. The hexagonal-lattice tight-binding model possesses the same symmetry and the same number of subbands as the flat band in TBG. We match the energy and length scales by choosing the width of a single band $W$ and the hexagonal lattice period $L_{M}$ identical to the parameters in TBG: $W=3.75$ meV and $L_{M}=a/2\sin(\theta/2)$ is the moiré superlattice periodicity. For the magic angle value $\theta=1.05\degree$, using carbon spacing $a=0.246$ nm, this gives $L_{M}=13.4$ nm. To ensure that a unit cell of the toy model can accomodate 4 electrons just as the moiré cell does in TBG, we make the toy model 4-fold degenerate. Comparison with plasmons for the actual TBG model, presented below, will help us to identify the features that are general as well as those which are a specific property of TBG.

Our nearest-neighbor tight-binding Hamiltonian is

[TABLE]

with the hopping matrix element $W/3$ to nearest neighbors at positions ${\boldsymbol{e}}_{j}=(\cos(2\pi j/3),\sin(2\pi j/3))L_{M}/\sqrt{3}$, $j=0,1,2$. Here, $W$ is the bandwidth measured from Dirac point, and the nearest neighbor distance $L_{M}/\sqrt{3}$ is chosen such that the lattice period of the hexagonal toy model matches the moiré superlattice period. Corresponding energies $E_{s,{\boldsymbol{k}}}$ and eigenstates $\Psi_{s,{\boldsymbol{k}}}$ are then

[TABLE]

where $\varphi_{{\boldsymbol{k}}}={\rm arg}\,h_{{\boldsymbol{k}}}$ and the band index $s=\pm$ labels the conduction and valence band.

Plasmons can be obtained from the nodes of the complex dielectric function, describing the dynamical response of a material to an outside electric perturbation:

[TABLE]

Here, $V_{{\boldsymbol{q}}}=2\pi e^{2}/\kappa q$ is the Coulomb interaction in a medium with a background dielectric constant $\kappa$, and $\Pi(\omega,{\boldsymbol{q}})$ is the electron polarization function. The relation in Eq. (7) is exact as long as the polarization function is defined as an exact microscopic density-density pair correlator given by a sum of all irreducible bubble diagrams. As such, this relation can yield useful information about plasmon dispersion, even when electron interactions are strong.

Similar to the conventional analysis of plasmons in 2D systems, here a simplification occurs in the small-$q$ limit, regardless of whether the random-phase approximation (RPA) is used to evaluate $\Pi(\omega,{\boldsymbol{q}})$. Indeed, since the Coulomb potential diverges at small $q$, zeros of $\varepsilon(\omega,{\boldsymbol{q}})$ are found when the polarization function is small. However, at small $q$, this quantity vanishes as $\lambda q^{2}/\omega^{2}$, a behavior that is a consequence of the general symmetry requirements (namely, gauge invariance demanding that spatially uniform external potential does not perturb density) mahan . This immediately yields a $q^{1/2}$ scaling for plasmon frequency at small-enough $q$.

Below, we use the RPA approach to estimate the prefactor $\lambda$ and to demonstrate that the mode $\omega\sim q^{1/2}$ extends far above the TBG p-h continuum. To compare with other systems, we recall the familiar “classical acceleration” behavior found for particles with parabolic dispersion: $\Pi(\omega,{\boldsymbol{q}})=nq^{2}/m\omega^{2}$, where $n$ is the charge density and $m$ is the electron band mass mahan . For a more general band dispersion, the ratio $n/m$ is replaced by the band Fermi energy, $\lambda\sim{E_{F}}/\hbar^{2}$ wunsch ; hwang ; jablan . Interactions have no impact on the behavior of $\Pi(\omega,{\boldsymbol{q}})$ for the parabolic band case; however, for nonparabolic bands, the band mass $m$ must change to an effective value $m^{*}$ described by Landau Fermi-liquid renormalization levitov_plasmons .

In our case, the scaling relation $\Pi(\omega,{\boldsymbol{q}})\approx\lambda q^{2}/\omega^{2}$ features different values of $\lambda$ for low and high energies, $\omega\lesssim E_{F}$ and $\omega>2W$. To see this, we start with the RPA expression for polarization function

[TABLE]

Here, summation $\sum_{{\boldsymbol{k}}}$ denotes integration over the Brillouin zone (BZ), the indices $s$,$s^{\prime}$ run over the electron bands and the factor of $4$ in front of the summation accounts for the 4-fold degeneracy of the toy model. Here, $f_{s,{\boldsymbol{k}}}$ is the equilibrium distribution $1/(e^{\beta(E_{s,{\boldsymbol{k}}}-E_{\rm F})}+1)$, and $F^{ss^{\prime}}_{{\boldsymbol{k}}+{\boldsymbol{q}},{\boldsymbol{k}}}$ describes band coherence factors. For our toy model,

[TABLE]

where $\Psi_{s,{\boldsymbol{k}}}$ are pseudospinors given in Eq. (6).

As we now show, an analytic expression for plasmon dispersion can be obtained, describing both the Landau-damped and the undamped cases in a unified way. We first rewrite Eq. (8) by performing a standard replacement ${\boldsymbol{k}}+{\boldsymbol{q}}\to-{\boldsymbol{k}}$ in the term containing $f_{s,{\boldsymbol{k}}+{\boldsymbol{q}}}$ followed by $-{\boldsymbol{k}}-{\boldsymbol{q}},-{\boldsymbol{k}}\to{\boldsymbol{k}}+{\boldsymbol{q}},{\boldsymbol{k}}$ justified by the ${\boldsymbol{k}}\to-{\boldsymbol{k}}$ time-reversal symmetry. This gives

[TABLE]

The behavior of this expression at small ${\boldsymbol{q}}$, which will be of interest for us, can be found in a closed form. In the small-$q$ limit the coherence factors behave as

[TABLE]

The values $O(1)$ for intraband transitions and $O(q^{2})$ for interband transitions might suggest that the polarization function is dominated by the intraband transitions. However, as we now show, the interband and intraband contributions are of the same order of magnitude.

Indeed, the intraband contributions, $s=s^{\prime}$, can be rewritten by noting that, upon integration over ${\boldsymbol{k}}$, only the even-${\boldsymbol{k}}$ part of series expansion $E_{s,{\boldsymbol{k}}+{\boldsymbol{q}}}-E_{s,{\boldsymbol{k}}}$ survives, giving $\Pi_{1}(\omega,{\boldsymbol{q}})\approx\frac{4}{\omega^{2}}\sum_{{\boldsymbol{k}},s}f_{s,{\boldsymbol{k}}}\left(E_{s,{\boldsymbol{k}}+{\boldsymbol{q}}}+E_{s,{\boldsymbol{k}}-{\boldsymbol{q}}}-2E_{s,{\boldsymbol{k}}}\right)$. Expanding in small $q$, we have

[TABLE]

As a sanity check, for parabolic band $E_{{\boldsymbol{k}}}=k^{2}/2m$ this yields the familiar “classical acceleration” result $\Pi(\omega,{\boldsymbol{q}})=\frac{nq^{2}}{m\omega^{2}}$.

The interband contributions, $s=-s^{\prime}$, can be simplified by noting that $E_{s,{\boldsymbol{k}}+{\boldsymbol{q}}}\approx-E_{s^{\prime},{\boldsymbol{k}}}$, giving

[TABLE]

As a sanity check, at $T=0$ the imaginary part of $\Pi_{2}$, describing interband transitions, is nonzero only for $2E_{\rm F}<\omega<2W$, as expected. The real part of $\Pi_{2}$ is negative at small $\omega$ and positive at large $\omega$ because the valence band contribution dominates over that of the conduction band.

Plasmon dispersion $\omega_{\rm p}$ is given by the solution of the equation $\varepsilon(\omega,{\boldsymbol{q}})=0$ with $\Pi=\Pi_{1}+\Pi_{2}$. Comparing the $\omega$ dependence of $\Pi_{1}$ and $\Pi_{2}$, we see that at small frequencies, $\omega<2E_{\rm F}$, the intraband contribution $\Pi_{1}$ dominates. This gives the dispersion in Eq. (1) with

[TABLE]

where the leading term $\beta_{0}=4\alpha v_{\rm F}E_{\rm F}/\hbar$ originates from $\Pi_{1}$ (see supplemental information), and the subleading $q$-dependent contribution is due to $\Pi_{2}$. Negative sign of $\Pi_{2}$ translates into $\beta_{1}<0$, softening the dispersion at low frequencies. This behavior, which holds the limit $\omega<2E_{\rm F}$, agrees with refs. wunsch ; hwang ; principi2009 .

In the same manner, we can obtain the dispersion at high frequencies, $\omega>2W$ (the intrinsically undamped regime). The analysis is again simplified by noting that, since $\alpha=e^{2}/\hbar\kappa v_{\rm F}\gg 1$, the relevant values of $q$ are small compared to the Brillouin zone size, and thus, the small-$q$ limit considered above is sufficient to describe this behavior. Taking both the intraband and interband contributions in the asymptotic form $\Pi_{1}=\lambda_{1}q^{2}/\omega^{2}$, $\Pi_{2}=\lambda_{2}q^{2}/\omega^{2}$ where $\lambda_{1}\approx 2E_{\rm F}/\hbar^{2}\pi,\lambda_{2}\approx 2(W-E_{\rm F})/\hbar^{2}\pi$ (see supplemental information), yields Eq. (1) with $\beta=\frac{2\pi e^{2}}{\kappa}(\lambda_{1}+\lambda_{2})$. The first term is identical to $\beta_{0}$ found at low frequencies, the second term is of a positive sign, $\lambda_{2}>0$, describing stiffening of the plasmon dispersion due to interband transitions.

In the undamped regime, plasmon frequency peaks at $q$ values on the order of Brillouin zone scale. The peak value of $\omega_{\rm p}$, given in Eq. (2), can be found by estimating the energy differences $E_{s,{\boldsymbol{k}}+{\boldsymbol{q}}}-E_{s^{\prime},{\boldsymbol{k}}}$ in Eq. (10) as $W$ and noting that the coherence band factor for large $q$ is in general non-vanishing and of order $1$. This gives, for the practically interesting case of $E_{\rm F}\sim W$, the result $\omega_{\rm p}\sim\sqrt{\alpha}W$, which agrees with the dispersion $\omega_{\rm p}=\sqrt{\beta q}=2\sqrt{\alpha v_{\rm F}Wq/\hbar}$ provided that $\hbar v_{\rm F}q$ saturates at $W$. Indeed, the estimated values of $\beta_{0}$,$\beta$ compared with the fitted curve in Fig. 1a (see supplemental information) indicate that $\omega_{\rm p}=\sqrt{\beta_{q}q}$ relation from Eq. (1) is a good approximation for the plasmon dispersion at both small and large $q$.

The dielectric function of the 2-band toy model faithfully reproduces all of the qualitative features expected for the TBG bandstructure. However, we find that, despite matching the bandwidth $W$ and lattice period to those of TBG, the resulting plasmon dispersion extends to much higher energies then those that will be found below for the actual TBG bandstructure. This is simply because the 2-band model does not account for the effects of interband polarization of higher electron bands, which renormalize the dielectric constant down and soften the plasmon dispersion. We account for this in the toy model case by rescaling the effective fine structure constant such that the resulting plasmon dispersion is comparable in magnitude with the TBG result. Specifically, in Fig. 1a, we use an effective background dielectric constant $\kappa=12.12$, which is 4 times larger than the dielectric constant $\kappa=3.03$ corresponding to an air/TBG/hexagonal boron nitride (hBN) heterostructure.

Next, we turn to the analysis of plasmons in TBG flat bands at an experimentally relevant magic angle value $\theta=1.05^{\circ}$ cao2018_1 ; cao2018_2 ; yankowitz2019 . To accurately describe the TBG band structure and eigenstates, we use the effective continuum Hamiltonian $H_{{\mathrm{T}BG}}$ introduced in ref. koshino2018 . The full discussion of the band structure details can be found in the supplemental material; here, we only discuss 2 relevant energy scales: flat-band bandwidth $W$ and the gap $\Delta$ between the flat bands and the rest of the band structure. With regard to $W$ value, we note that, technically, the bandwidth of the flat-bands, as predicted by the continuum mode $H_{{\mathrm{T}BG}}$, is on the order of $W\approx 3.75$ meV. However, the bandwidth scale relevant for the interband and intraband excitations is actually closer to $\tilde{W}\approx 2$ meV, because most of the states in the band lie below $2$ meV. In addition, since the states with energies outside $-2\,{\rm meV}<E<2\,{\rm meV}$ are small $k$, their contribution to polarization function [Eqs. (12) and (13)), evaluated at small $q$, is small. We also note that, while the bandgap as predicted by the continuum model is $\Delta\approx 11.75$ meV, the actual gap is still a subject of debate tomarken2019 .

The definition of the polarization function for the TBG continuum model is essentially identical to that of the tight binding toy model [Eq. (8)]. Now, however, we must account explicitly for the valley and spin degrees of freedom, for a larger number of electron bands, and for different coherence factors. Accordingly, we promote the band indices $s$,$s^{\prime}$ in Eq. (8) to composite labels $n$,$m$, which label all electron bands, spins $\sigma$ and valleys $\xi$; this makes the additional factor of $4$ in front of Eq. (8) redundant. The toy model coherence factors are replaced by the TBG coherence band factors $F^{nm}_{{\boldsymbol{k}}+{\boldsymbol{q}},{\boldsymbol{k}}}$, which are given by

[TABLE]

where $\Psi_{n,{\boldsymbol{k}}}({\boldsymbol{r}})$ are the Bloch wavefunctions for momentum ${\boldsymbol{k}}$ and band/valley/spin composite label $n$, which diagonalize the continuum Hamiltonian (see supplemental materials). The integral in Eq. (15) is carried over the moiré unit cell $\Omega$.

After the polarization function is evaluated, we can determine the dielectric function and identify TBG’s collective modes from poles of $1/\varepsilon(\omega,{\boldsymbol{q}})$ as above. An example of a TBG’s dielectric function at approximately half-filling of the electron band, $E_{\rm F}=0.289$ meV, is shown in Fig. 3; fixed $q$ line cuts and zeros of $\varepsilon(\omega,{\boldsymbol{q}})$ are illustrated in the supplemental materials. In discussing the figure, it is helpful to contrast it with the calculation for the hexagonal-lattice toy model shown in Fig. 1a. We again see a well-defined intrinsically undamped plasmon mode $\omega_{\rm p}$ (red in Fig. 1a) positioned above the p-h continuum; the mode resides inside the band gap $2W<\omega_{\rm p}<W+\Delta$, which peaks at $\hbar\omega_{\rm p}\approx 8.5$ meV before decreasing and becoming almost flat $\hbar\omega_{\rm p}\approx 6.5$ meV at large momenta. In agreement with the analytic considerations above, we see the interband continuum extending from $2E_{\rm F}$ to $2W$, but since $E_{\rm F}=0.289$ meV is extremely small, it makes the conventional (Landau-damped) part of plasmon dispersion $\omega<2E_{\rm F}$ invisible on the figure.

There are several unique aspects of the TBG plasmon dispersion compared with the behavior of generic narrow-band plasmons discussed above. To analyze the dispersion at $\omega_{\rm p}>2W$, we proceed just as in the toy model case, rewriting the TBG polarization function in a slightly different form of Eq. (10), where the indices $n$, $m$ and the band coherence factor are modified as described above.

To proceed further analytically, we need to analyze Eq. (10) in the long-wavelength limit. However, unlike the 2-band toy model, where the only characteristic energy scale was the bandwidth $W$, the TBG band structure features an additional energy scale, namely, the gap between the flat bands and the rest of the energy spectrum. This impacts the small-$q$ series expansion of the polarization function, as now the energy difference $E_{n}-E_{m}$ between the occupied and unoccupied states can be larger than $\omega$. To account for such contributions in the series expansion, we split the summation over TBG bands into 2 parts, depending on whether $\omega$ or the energy difference $E_{n}-E_{m}$ is the largest energy scale in the denominator of Eq. (10). This yields an approximate expression for the dielectric function

[TABLE]

where we defined 2 auxiliary functions:

[TABLE]

and

[TABLE]

Here the band summations $\sum^{\prime}$ and $\sum^{\prime\prime}$ run over bands such that $\omega^{2}>(E_{n,{\boldsymbol{k}}+{\boldsymbol{q}}}-E_{m,{\boldsymbol{k}}})^{2}$ and $\omega^{2}<(E_{n,{\boldsymbol{k}}+{\boldsymbol{q}}}-E_{m,{\boldsymbol{k}}})^{2}$, respectively: for example, at large momenta, as seen in Fig. 3, the plasmon mode lies in the gap between the flat and non-flat bands, and hence, the $B({\boldsymbol{q}})$ summation extends only over the flat bands, whereas the summation in $A({\boldsymbol{q}})$ includes all of the remaining combinations of band indices. This allows us to write a closed form expression for the plasmon dispersion as

[TABLE]

which must hold for both small and large $q$. We consider these 2 limits separately.

At small $q$, the matrix element of the Bloch wavefunctions, just as in the toy model case, favors the overlap between states from the same band. At the same time, there are fewer states in the $A({\boldsymbol{q}})$ satisfying the condition $\omega^{2}>(E_{n,{\boldsymbol{k}}+{\boldsymbol{q}}}-E_{n,{\boldsymbol{k}}})^{2}$, and hence, $A({\boldsymbol{q}})$ vanishes for small $q$. This amounts to the plasmon dispersion $\omega_{\rm p}$ from Eq. (19) reducing to $\omega_{\rm p}^{2}\approx B({\boldsymbol{q}})$, and by comparison with Eq. (12), we similarly expect a conventional 2D plasmon dispersion $\omega_{\rm p}=\sqrt{\beta_{q}q}$ with $\beta_{q}$ given by the series from Eq. (14). As we see in Fig. 3, the $\omega_{\rm p}=\sqrt{\beta_{q}q}$ dispersion is a valid description only at very small $q$ compared to the Fig. 1a, which can be traced back to higher bands softening the plasmon dispersion through the $A({\boldsymbol{q}})$ term in Eq. (19).

To determine how high the plasmon mode rises above the p-h continuum, we consider large $q$ values comparable to the reciprocal lattice vector. The arguments similar to those in the toy model show that, since $\alpha\gg 1$, we have $A({\boldsymbol{q}})\gg 1$. The dependence on the $e^{2}/\kappa q$ ratio, therefore, cancels between the $A({\boldsymbol{q}})$ and $B({\boldsymbol{q}})$ functions, resulting in the value of the plasmon dispersion $\hbar\omega_{\rm p}\approx\sqrt{B({\boldsymbol{q}})/A({\boldsymbol{q}})}\sim\sqrt{W\Delta}\approx 6.6$ meV being dictated only by the continuum model’s band structure parameters. This lack of explicit dependence on $\alpha$ suggests that, after the doping is such that $\alpha\gg 1$, the large-$q$ value of $\hbar\omega_{\rm p}\approx\sqrt{W\Delta}$ becomes insensitive to doping (and hence, Fermi velocity). This behavior is different from that in the toy model, where $\omega_{\rm p}\sim\sqrt{\alpha}W$ at large $q$. The relatively more weak dependence on $\alpha$ in the TBG case is due to interband polarization involving higher bands, which significantly alters the effective dielectric constant. The weak $q$ dependence at large $q$ is in agreement with the properties of interband plasmons described in ref. stauber2016 .

We also note that, although plasmons above the p-h continuum are kinematically protected from p-h excitation, which makes them undamped at the RPA level, there exist relaxation pathways through higher-order pair production in which several electron-hole pairs are emitted with total energy exceeding $\tilde{W}$, as well as phonon-assisted processes. For conventional plasmons these processes were analyzed in ref. glazman_damping . The role of these effects for plasmon lifetimes in TBG will be a subject of future work.

Before closing, we note that suppressing damping has always been central to the quest for tightly-confined low-loss surface plasmon excitations. An early approach utilized surface electro-magnetic modes traveling at the edge of an air/metal boundary Raether , in which dissipation is low because most of the mode field resides outside the metal; however, the field confinement scale, set by optical wavelength, was fairly large. Next came surface plasmons propagating in high-mobility 2D electron gases in semiconductor quantum wells and monolayer graphene woessner2014 , which can provide deep-subwavelength confinement jablan . However, plasmons in these systems are prone to a variety of dissipation mechanisms, with Landau damping usually regarded as the one that sets the fundamental limit on possible plasmon wavelengths and corresponding lifetimes. The possibility to overcome this fundamental limitation in narrow-band systems, such as moiré graphene, creates a unique opportunity for graphene plasmonics. Damping-free plasmons can enable novel interference phenomena, dissipationless photon-matter coupling, and other interesting behaviors. It is also widely expected that low-dissipation plasmons can lead to unique applications for photon-based quantum information processing gullans2013 . Furthermore, reduced damping has more immediate consequences, as it translates into enhanced optical coherence that can be directly probed by scanning near-field microscopy, as discussed above, providing a clear signature of the undamped collective modes.

We thank Ali Fahimniya for useful discussions. This work was supported, in part, by the Science and Technology Center for Integrated Quantum Materials, NSF Grant No. DMR-1231319; and Army Research Office Grant W911NF-18-1-0116.

Supplemental Information

I Spatial speckle patterns in near-field optical microscopy

Here we elaborate on the analysis connecting Eq. (3) and the speckle patterns shown in Fig. 1d and Fig. 2a-d. For an in-depth discussion of the near-field optical microscopy measurement technique and quantitative modeling of the detected signal we refer the reader to refs. koppens2012 ; basov2012 ; woessner2014 .

As argued in the main text, we can estimate the strength of the measured signal in the near-field opical microscopy by evaluating the equal-point correlation function Eq. (3), which here we restate for convenience:

[TABLE]

It describes an amplitude of plasmon excitation, which traveled from the tip at position ${\boldsymbol{r}}$ to a disorder at position ${\boldsymbol{r}^{\prime}}$ and was then reflected back towards the tip at ${\boldsymbol{r}}$. Here the Green’s function $G_{E}(r)$ of the plasmon excitation of wavenumber $q_{0}$ is taken in the limit $rq_{0}\gg 1$ as

[TABLE]

which describes radially propagating waves in 2D. The factor $e^{-\delta|{\boldsymbol{r}}|}$ describes damping due to extrinsic effects such as phonons and other inelastic processes. Upon substitution of the Green’s function into (S.1), the measured signal $S({\boldsymbol{r}})$ is given by

[TABLE]

This expression, which is a convolution of two functions, will generate a product under Fourier transform.

For purposes of Fig. 1d and Fig. 2a-d we evaluate the above convolution numerically by using the convolution theorem, that is first performing a fast Fourier transform of both terms individually, multiplying them and then carrying out an inverse Fourier transform. The inset of Fig. 1d and Fig. 2a-d is the intermediate step of this process, but we can also determine it analytically by evaluating the Fourier transform of the signal $S({\boldsymbol{r}})$

[TABLE]

As expected by the convolution theorem the expression factorizes into a product of two separate factors

[TABLE]

where the first factor is nothing but the Fourier harmonic of $\eta({\boldsymbol{x}})$ and the second factor is the $r$-$r^{\prime}$ influence function, simplified by performing a variable change ${\boldsymbol{r}}-{\boldsymbol{r}}^{\prime}\to{\boldsymbol{r}}$. To evaluate the integral over $d^{2}{\boldsymbol{r}}$ we first integrate over $|{\boldsymbol{r}}|$ and then carry out angular integration using the identity

[TABLE]

After substituting $a=k$, $b=2k_{0}+2\delta i$ this gives Eq. (4) of the main text.

II Behavior of the intraband and interband polarization functions

Here we discuss the behavior of the intraband and interband polarization functions $\Pi_{1}(\omega,{\boldsymbol{q}})$ and $\Pi_{2}(\omega,{\boldsymbol{q}})$ of the toy model, defined in Eqs. (12), (13) of the main text. In particular, we estimate the coefficients $\lambda_{1}$ and $\lambda_{2}$ describing the small-$q$ behavior of $\Pi_{1}$ and $\Pi_{2}$, defined in the paragraph beneath Eq.(14). We are mostly interested in high frequency values $\omega>2W$ describing the intrinsically undamped regime.

We start with the quantity $\lambda_{2}$ describing the contribution of intraband transitions. At small $q$, the interband coherence factor from Eq. (9) is non-negligible only in proximity of the points $K$ and $K^{\prime}$. Near these points a linear dispersion $E_{s,{\boldsymbol{k}}}=sv_{F}k$, with $s=\pm 1$, is a good approximation for the bandstructure. In that limit, the small-$q$ interband coherence band factor from Eq. (11) becomes

[TABLE]

where $\theta$ is the angle between ${\boldsymbol{k}}$ and ${\boldsymbol{q}}$. The quantity $\Pi_{2}(\omega,{\boldsymbol{q}})$ is therefore given by:

[TABLE]

In the above we used the linear dispersion approximation $E_{s,{\boldsymbol{k}}}=sv_{F}k$ for the whole band and accounted for the $K$ and $K^{\prime}$ points through an additional factor of $2$. This gives

[TABLE]

with the first and second terms originating from the conduction band and the valence band respectively. This gives

[TABLE]

which takes positive values since $-W<E_{F}<W$.

Next, we proceed to estimate the $\lambda_{1}$. Without loss of generality, we place the Fermi energy in the conduction band. In this case, the interband contribution to the polarization function is non-vanishing only in the conduction band. This can be seen by going back to the Eq. (12), which for $s=s^{\prime}=-1$ and small $q$ vanishes:

[TABLE]

since $f_{-1,{\boldsymbol{k}}}=1$ for all ${\boldsymbol{k}}$ in the valence band. It is therefore sufficient to focus on the contribution of the partially filled (conduction) band. To be consistent with the $\lambda_{2}$ analysis above we replace the dispersion energy as $E_{1,{\boldsymbol{k}}}=v_{F}k$. The intraband contribution to the polarization function $\Pi_{1}(\omega,{\boldsymbol{q}})$ is then

[TABLE]

giving

[TABLE]

As argued in the main text [see discussion below Eq. (14)], this result remains unchanged for frequencies $\omega<2E_{F}$ and, therefore,

[TABLE]

Going back to the $\omega>2W$ regime, and using $\lambda_{1}$ and $\lambda_{2}$ derived above, gives the square-root plasmon dispersion $\omega_{\rm p}=\sqrt{\beta q}$ with

[TABLE]

Therefore, at small $\omega<2E_{F}$ the dispersion behaves as $\omega_{\rm p}=2\sqrt{\alpha v_{F}E_{F}q}$, becoming enhanced at high energies $\omega>2W$, $\omega_{\rm p}=2\sqrt{\alpha v_{F}Wq}$ by a factor $\sqrt{W/E_{F}}$.

To complete the analysis of the polarization function behavior, now we focus on frequencies in the region $2E_{F}<\omega<2W$. Working again in the small-$q$ limit we find that, as pointed out earlier, only the interband contribution to the polarization function $\Pi_{2}(\omega,{\boldsymbol{q}})$ develops an imaginary part, whereas the intraband polarization function $\Pi_{1}(\omega,{\boldsymbol{q}})$ is real-valued, given by the Eq. (S.13). To determine the form of $\Pi_{2}(\omega,{\boldsymbol{q}})$ in the interband p-h continuum energy range, we approximate the coherence factor as in Eq. (S.7) to obtain

[TABLE]

Here we used the linear approximation to the energy dispersion $E_{s,{\boldsymbol{k}}}=sv_{F}k$, accounting for the fact that, because of the behavior of the coherence factors, only the states near the Dirac point contribute to $\Pi_{2}$. As always, we account for the $K$ and $K^{\prime}$ points by an additional factor of $2$. After carrying out integration over $d^{2}k$ we arrive at:

[TABLE]

Here $\Theta(x)$ is the Heaviside function, which ensures that the imaginary part is non-zero only in the particle-hole continuum region $2E_{F}<\omega<2W$. The dielectric function in this region is therefore

[TABLE]

which shows that the collective mode $\omega_{\rm p}$ in the $2E_{F}<\omega_{\rm p}<2W$ region has a damped square-root dispersion

[TABLE]

The imaginary part, which scales linearly with $q$, describes damping due to particle-hole pair production.

We finish the discussion of the collective modes by comparing the analytically predicted dispersion with the numerical result in Fig. 1a. While the simulated dispersion closely follows the square-root dependence $\omega_{\rm p}\propto\sqrt{q}$, the agreement between the simulation and $\omega_{\rm p}=\sqrt{\beta_{q}q}$ dispersion is drastically improved if the two first terms $\beta_{0}$ and $\beta_{1}$ from the series expansion in Eq. (14) are used for a fitting. Although the terms $\beta_{0}$ and $\beta_{1}$ could in principle be computed by carrying out an expansion of the polarization function in Eq. (10) in powers of $q$ and then evaluating the resulting integrals numerically, we instead treat $\beta_{0}$ and $\beta_{1}$ as free parameters and fit them to the simulated dispersion. This approach yields values

[TABLE]

The best-fit $\beta_{0}$ value is close to $\beta_{0}=4\alpha v_{F}E_{F}\approx 0.86\times 10^{3}$ ${\rm meV}^{2}~{}{\rm nm}$ predicted from Eq. (S.15). We also see that, since $\beta_{1}$ is negative, the plasmon dispersion is indeed softened by interband polarization, in agreement with the argument given in the main text [see Eq. (14)].

III Twisted bilayer graphene - details of the model

Here we describe in detail the model for twisted bilayer graphene (TBG) bandstructure used in the main text. We use the effective continuum Hamiltonian introduced in ref. koshino2018 , adopting notations and numerical values used in ref. koshino2018 .

The continuum approach is made possible by the small values of the twist angle $\theta$ by which the two graphene layers in TBG are rotated relative to one another. We start by taking two AA-stacked graphene layers and rotating the layer 1 and the layer 2 around the B-sites by $-\theta/2$ and $\theta/2$ respectively. For the “magic” value of $\theta=1.05^{\circ}$, the moiré real-space lattice constant is $L_{M}=a/2\sin(\theta/2)\approx 13.4$ nm. This is two orders of magnitudes greater than the graphene’s lattice constant $a=0.246$ nm, justifying the use of the continuum approach.

In momentum space this real-space rotation translates into two graphene Brillioun zones rotated by angle $\theta$ relative to each other. Both BZs are centered at the same $\Gamma$ point but the $K$ (and $K^{\prime}$) points of the two layers are separated by a small momentum $4\pi/(3L_{M})$. As the moiré periodicity $L_{M}$ is much greater than the lattice constant $a$, we can ignore the intervalley mixing between the two valleys $K$ and $K^{\prime}$ of the original graphene layers - here labeled by $\xi=-1,1$. The total Hamiltonian of the system becomes therefore block diagonal in the valley index. The blocks $H^{(\xi)}$ describing each of the two valleys take the form

[TABLE]

in the basis of $(A_{1},B_{1},A_{2},B_{2})$ sites. The matrices $H_{l}$ ($l=1,2$) correspond to the intralayer Hamiltonians of the layers. The latter, due to the lengthscale separation between $L_{M}$ and $a$, can be approximated by performing the standard $kp$ expansion around the points $K$ and $K^{\prime}$.

This procedure gives $2\times 2$ Dirac Hamiltonians centered at the ${\boldsymbol{K}}_{\xi}^{(l)}$ points

[TABLE]

where ${\boldsymbol{k}}$ is a momentum in the BZ of the original graphene layers, and $R\left(\varphi\right)$ is the $2\times 2$ rotation matrix

[TABLE]

that accounts for rotation of the BZ of the original graphene layers. The signs $\pm$ in Eq. (S.23) correspond to the layers $l=1$ and $2$, respectively.

The energy scale for the Hamiltonians $H_{l}$ is $\hbar v/a=2.1354$ eV. The vectors ${\boldsymbol{K}}_{1}^{(l)}$, ${\boldsymbol{K}}_{-1}^{(l)}$, which denote the Dirac points $K$ and $K^{\prime}$ of the layers, are given by

[TABLE]

respectively. We stress that, while ${\boldsymbol{k}}$ alone has length close to $\sim 4\pi/3a$, the difference ${\boldsymbol{k}}-{\boldsymbol{K}}_{\xi}^{(l)}$ is small, since ${\boldsymbol{k}}$ is always located near the vicinity of the ${\boldsymbol{K}}_{\xi}^{(l)}$ points. This makes the linear expansion from Eq. (S.23) a well defined approximation.

More quantitatively, the expressions in Eq. (S.23), found by Taylor expanding the graphene tight-binding Hamiltonian, are valid for momenta close enough to the Dirac points of the two layers, $|{\boldsymbol{k}}-{\boldsymbol{K}}_{\xi}^{(l)}|a\ll 1$. For $\theta\ll 1$ this condition is obeyed in the entire mini Brillouin zones of the TBG superlattice.

In the analysis below the moiré superlattice BZ is defined as in the inset of Fig. 3, with the two reciprocal lattice vectors

[TABLE]

We denote the reciprocal lattice vector length as $G_{M}=|{\boldsymbol{G}}_{1}^{M}|=|{\boldsymbol{G}}_{2}^{M}|=4\pi/\sqrt{3}L_{M}$. Matrix $U$ is the effective moiré interlayer coupling given by:

[TABLE]

where we introduced a notation for the phase factor $\nu=e^{i2\pi/3}$. The interlayer couplings $u$ and $u^{\prime}$ are taken as $u=0.0797$ eV and $u^{\prime}=0.0975$ eV to match values in ref. koshino2018 .

To determine the energy bands and the eigenstates we take the Bloch wavefunction ansatz for a valley $\xi$ as

[TABLE]

with $X$ labeling the spinor components $X=A_{1}$, $B_{1}$, $A_{2}$, $B_{2}$. The band index, labeled by $n$ and ${\boldsymbol{k}}$, is the Bloch wave vector in the BZ of the original graphene layers. Here ${\boldsymbol{G}}$ runs over all possible integer combinations of the reciprocal lattice vectors, ${\boldsymbol{G}}=m_{1}{\boldsymbol{G}}_{1}^{M}+m_{2}{\boldsymbol{G}}_{2}^{M}$ with integer $m_{1}$ and $m_{2}$. As discussed in ref. koshino2018 , the low-energy states are expected to be dominated by states near the original Dirac points. We therefore take only not-too-large indices $m_{1}$ and $m_{2}$ that satisfy the condition

[TABLE]

where $z$ is a conveniently chosen number of order one [ref. koshino2018 uses $z=4$], and ${\boldsymbol{M}_{\xi}}$ are the “mean” Dirac point locations

[TABLE]

given by the midpoint between the $K$ (or $K^{\prime}$) points of the two layers.

IV Electron loss function for the TBG bandstructure

Fig. 4 details the behavior of the electron loss function for TBG, depicted in Fig. 3 of the main text. Panels a and b show constant-momentum ${\boldsymbol{q}}$ linecuts of the real and imaginary parts of the dielectric function $\varepsilon(\omega,{\boldsymbol{q}})$. The finite width of the plasmon resonance in the loss function in Fig. 4c is due to the infinitesimal imaginary part of $\omega+i0$ in the polarization function in Eq. (8) replaced with $\omega+i\gamma$, with a suitably chosen small $\gamma$ introduced for illustration purposes.

Strong electron-electron interactions in the narrow electron bands lead to large dielectric function values, as can be seen in Fig. 4. For energies $\hbar\omega<2W$ the dielectric function imaginary and real parts take values a few orders of magnitude higher than those of graphene monolayer. The origin of these large values can be traced to the high effective fine structure constant (or, equivalently, low Fermi velocity) in the flat electron bands, as discussed in the main text. To see this in more detail, we recall the Thomas-Fermi expression for the long-wavelength static dielectric function of graphene hwang

[TABLE]

with the Thomas-Fermi momentum $q_{TF}=N\alpha k_{F}$, where $N$ is the degeneracy factor $N=8$ ($2$ spins, $2$ layers, $2$ valleys). For illustration purposes, taking a fine structure constant $\alpha\sim 30$ and Fermi momentum $k_{F}\sim{K}$, for the momentum $q\sim K/2$ (red line in the Fig. 4) Eq. (S.31) predicts a dielectric function value $\varepsilon\sim 480$, which is in good agreement with the simulation results. Above $\hbar\omega>2W$ the dielectric function rapidly decreases until $\hbar\omega>20$ meV where the contributions of higher electron bands start to dominate.

At these energies, plasmon dispersion is strongly affected by the presence of higher electron bands. At small $q$ plasmon dispersion is predominantly due to intraband transitions, and is thus insensitive to other electron bands. At large $q$ the situation changes. In the absence of higher electron bands the zeros of the dielectric function would occur at much larger energy scales $\hbar\omega_{\rm p}\sim 40$ meV. However, as argued in Eq. 19 in the main text, higher electron bands push plasmon dispersion down with the large-$q$ zeros of the dielectric function on the order $\hbar\omega_{\rm p}\approx\sqrt{W\Delta}\approx 6.6$ meV. Here $W$ is the flat-band bandwidth and $\Delta$ is the band gap as defined in the main text. The independence of this value of $\alpha$ is the behavior to be expected for large enough $\alpha$, such that plasmon dispersion extends above the p-h continuum. The independence of $\omega_{\rm p}$ of $\alpha$ at $\alpha\gg 1$ is a characteristic feature of interband plasmons.

Lastly, we note that our simulation is expected to be accurate only for ${\boldsymbol{q}}$ inside the TBG Brillouin zone. When ${\boldsymbol{q}}$ approaches zone boundary it is necessary to consider local field effects adler1962 ; wiser1963 . Although these effects are often small, they require careful examination and thus will be a subject of a future work.