Optical conductivity, Drude weight and plasmons in twisted graphene bilayers

Since the realization that neutral graphene monolayers display a universal optical absorption of πα ≈ 2.3% in the infrared and visible spectrum (with α ≈ 1/137 the fine structure constant) [1–4], interest in the interactions between light and graphene has grown rapidly [5–7]. A remarkable aspect of this absorption, other than its suggestive relation to the universal constant α, is its high value, given the ultimate thinness of a graphene monolayer. In fact, it can be even enhanced up to 10% for photon energies around 4.5 eV due to transitions close to the van Hove singularity [8, 9].

The unique feature of neutral graphene responsible for its strong optical response is the existence of low energy particle–hole interband excitations with zero momentum, by virtue of the gapless Dirac spectrum, that act as an effective sink for photons of arbitrarily low energy. By the same token, all thermally excited plasmons in undoped graphene are damped by this excitation density, but nevertheless remain well-defined [10]. The unique bipolar tunability of graphene further allows one to adjust these optical properties quite simply, by shifting the Fermi energy μ away from the Dirac point, be it electrostatically or through chemical doping. This opens a gap of energy 2μ in the particle–hole excitations, making graphene transparent below this energy, useful, e.g. for broadband optical modulation [11]. It also allows for the emergence of undamped plasmons [12–16].

The spectral features involved in the unusual optical properties outlined above are also present in graphene bilayers [17]. These are often classified according to their stacking, which in the absence of strain have three possible (meta) stable configurations: AB'-type (or BA', also known as Bernal stacking), AA'-type and twisted. The former two are minimal stackings (unit cell with only four carbon atoms) in which half or all sites in opposite layers are vertically aligned, respectively. The twisted bilayer is the most generic type, and is defined by a relative rotation of angle θ between the two layers. Each of the three stacking types exhibits a different low energy electronic structure, although all three react strongly to light like the monolayer. The twisted bilayer, however, and in particular the case of low twist angles (around and below 1°), has the richest optical response, a consequence of its complex low energy spectral properties. The response, moreover, is highly tuneable through the angle θ, an additional control parameter not available in the other graphenes. As a striking example, the strong frequency-dependence of twisted bilayer conductivity and thus reflectivity leads to a visible coloration under white light, that varies with angle θ, of flakes on 100 nm thick SiO₂, as recently observed experimentally [18].

In this work we study the optical response of twisted graphene bilayers at different twist angles and temperature responsible for this effect. We analyze its optical conductivity at and away from neutrality, the Drude weight and plasmon spectrum. Some of these aspects, specifically the real part of the conductivity, have been studied in recent works for large-angle twisted bilayers [19, 20]. We will make connection with these, and extend them into the low angle regime.

The paper is organized as follows. In section 2, we introduce the continuum model of twisted graphene bilayer and discuss the electronic structure. In section 3, we introduce the Kubo formalism and also comment on the regularization of the Kramers–Kronig relation for the conductivity, necessary in the case of unbounded Dirac fermions. In sections 4 and 5, we present our results for the optical conductivity and the related Drude weight for different angles, chemical potential and temperature. In section 6, we discuss the plasmonic spectrum based on the local conductivity by plotting the energy loss function and close with conclusions.

Twisted graphene bilayers are a peculiar structure from a crystallographic point of view, since for a generic (non-commensurate) angle the lattice is not exactly periodic, although crystalline order within each layer is extraordinarily robust. This is a consequence of the highly anisotropic binding forces in these systems, also characteristic of graphite. In principle only commensurate twist angles θ = θ_mn, where [21]

$$\begin{equation} \mathrm{cos}(\theta_{mn})=\frac{3m^2+3\,m n+n^2/2}{3m^2+3\,m n+n^2} \end{equation} \tag{ 1 }$$

for integer m,n are amenable to band structure calculations by application of the Bloch theorem.

It turns out, however, that the bandstructure at low energies depends continuously on the value of the angle θ_n,m itself, not on the specific n and m (note that the set of possible commensurate angles θ_mn is dense). This is a consequence of the fact that the physical interlayer distance d > 0.3 nm is greater than the intralayer carbon–carbon distance a_cc = 0.14 nm. In this regime, all properties associated to precise commensurability are efficiently averaged out using microscopic tight-binding models. Consequently, a simple, low energy continuum theory can be formulated that may be applied to an arbitrary twist angle θ, and obviates all crystallographic complications to a very good approximation [21, 22].

The first ingredient in the continuum model is the momentum shift between Dirac cones in different layers due to the rotation, ΔK = 2|K|sin(θ/2), where |K| = 4π/3a₀ and $a_0=\sqrt {3}a_{\mathrm { cc}}$ is the monolayer Bravais period. The second non-trivial ingredient is the smooth spatial modulation in the interlayer coupling. The twist θ produces a triangular Moiré pattern in the local stacking, of period L_M = a₀/[2 sin(θ/2)], see figure 1(a). Hence, the coupling between layers becomes position dependent. The modulation is smooth on the scale a₀, so the two valleys in each layer remain decoupled in the twisted bilayer. For each valley, the continuum model Hamiltonian reads [21, 23, 24]

$$\begin{equation} H=\left(\begin{array}{@{}cccc@{}} 0 & \hbar v_{\rm F} \Pi_+^\dagger & V^*(\boldsymbol{r}) & V^*(\boldsymbol{r}-\boldsymbol{\delta r}) \\ \hbar v_{\rm F}\Pi_+ & 0 & V^*(\boldsymbol{r}+\boldsymbol{\delta r}) & V^*(\boldsymbol{r}) \\ V(\boldsymbol{r}) & V(\boldsymbol{r}+\boldsymbol{\delta r}) & 0 & \hbar v_{\rm F} \Pi_-^\dagger\\ V(\boldsymbol{r}-\boldsymbol{\delta r}) & V(\boldsymbol{r}) & \hbar v_{\rm F}\Pi_- & 0 \\ \end{array}\right). \end{equation} \tag{ 2 }$$

This matrix is written in the local basis A,B,A',B' comprised of the two sites in the unit cell of each of the two layers. Here, the Fermi velocity is $v_{\mathrm { F}}=(\sqrt {3}/2)a_0\gamma _0/\hbar \approx 9\times 10^5\,{\mathrm { m}}\,{\mathrm { s}}^{-1}$ and $\Pi _\pm =[k_x+\mathrm {i}(k_y\mp \Delta K/2)]{\mathrm { e}}^{{\mathrm { i}}\pm \theta /2}$. The function $V(\boldsymbol {r})=\frac {1}{3}\gamma _\perp \sum _{i=1}^3\,{\mathrm { e}}^{{\mathrm { i}}\boldsymbol {g}_i\cdot \boldsymbol {r}}$ describes the periodic spatial variation of the interlayer coupling, with γ_⊥ ≈ 0.33 eV. The Moiré lattice vectors are denoted by a_1,2 (so L_M = |a_1,2|) and conjugate momenta are g_1,2, such that g_i·a_j = 2πδ_ij (we also define g₃ = 0). The AA' sublattice is centered at the origin, with the AB'/BA' sublattices offset given by $\pm\boldsymbol{\delta r}=(\boldsymbol{a}_1-\boldsymbol{a}_2)/3$. These offsets are used in equation (2) to write the sublattice-resolved interlayer coupling so that it alternates in space between AB', BA' and AA' types, following the Moiré pattern.

Zoom In Zoom Out Reset image size

We first review the bandstructure phenomenology of the twisted bilayer. The continuum model has two distinct regimes, denoted as large and small angles, each with very different features. The large angle regime, relevant for θ greater than a few degrees, is characterized by two low energy Dirac cones per valley at the superstructure K and K', with a reduced Fermi velocity [25–27]

$$\begin{equation} v_{\rm F}^*\approx v_{\rm F}[1-9\alpha^2+\mathcal{O}(\alpha^4)], \end{equation} \tag{ 3 }$$

where α = γ_⊥/6_M is an expansion parameter that represents the dimensionless interlayer coupling and _M(θ) = ℏv_FΔK/2 = (2π/3)ℏv_F/L_M is the energy scale associated to the Moiré period L_M. Also, a low-energy van Hove singularity at energy _vH ≈ _M − γ_⊥/3, see figure 1(b), is created [23, 28–31].

At an angle θ^m₁ ≈ 1.05°, known as the 'first magic angle', v*_F becomes zero, the van Hove singularity degenerates into a narrow miniband around zero energy. This is the onset of the small angle regime. From this point on, the electronic structure becomes considerably complex. The Fermi velocity becomes finite again, as the miniband moves away from zero energy. Higher energy van Hove singularities appear and approach zero at subsequent magic angles θ^m_2,3,4..., where the Fermi energy vanishes once more.

In figure 2, the density of states (DOS) is shown for various angles labeled by the index i

$$\begin{equation} \cos(\theta_i)=1-\frac{1}{2(3\mathrm{i}^2+3\mathrm{i}+1)}, \end{equation} \tag{ 4 }$$

obtained from equation (1) with m = n + 1 = i. The curves are normalized by $\rho _0=\frac {4\gamma _0}{\pi (\hbar v_{\mathrm { F}})^2}$ such that the DOS of two uncoupled graphene layers is just a straight line with slope 1 (dashed line). Note that the DOS is almost but not exactly symmetric for positive and negative energies as a consequence of the underlying (not electron–hole symmetric) Hamiltonian.

Zoom In Zoom Out Reset image size

In the upper panel, we show the DOS for 'large' angles θ with i = 10 (blue) and i = 15 (red). Around the neutrality point, the spectrum is characterized by a van Hove singularity followed by a step singularity, forming a closed shell. For small energies the system is thus approximately described by a hexagonal tight-binding model with renormalized band-width. This shell structure is somehow repeated for larger energies until it fades away and the DOS for decoupled bilayer graphene is recovered. Note that the weight of the van Hove singularity is slightly decreasing with increasing (intermediate) i.

For small angles around the first magic angle θ^m₁ = 1.05° (i = 31, lower panel in figure 2), more distinct van Hove singularities appear and the height of the first van Hove singularity is strongly enhanced with increasing i. A close-up of this energy regime in shown on the left hand side of figure 8.

The above phenomenology in the electronic bandstructure has an associated signature in the optical response of the system as the angle decreases. The basic quantity governing the linear optical response of any electronic system is the optical conductivity σ_αβ(ω) (where α,β = x,y), which gives the electric current arising in response to a weak and uniform, time-dependent electric field of frequency ω, E(t) = E(ω) e^iωt,

$$\begin{equation} j_\alpha(\omega)=\sum_\beta\sigma_{\alpha\beta}(\omega) E_\beta(\omega). \end{equation} \tag{ 5 }$$

The optical conductivity may be written, using linear response theory, in terms of the current response χ_jαjβ(q,ω) at zero-momentum,

$$\begin{equation} \sigma_{\alpha\beta}(\omega)=\mathrm{i}\frac{\mathrm{e}^2}{\omega+\mathrm{i}\delta}\chi_{j_\alpha j_\beta}(\boldsymbol{q}=0,\omega), \end{equation} \tag{ 6 }$$

where, as is conventional, δ denotes a positive infinitesimal. The current response function χ_jαjβ consists, in general, of two terms: the diamagnetic and the paramagnetic contribution [32]. The general formalism for obtaining the diamagnetic contribution in a tight-binding model is outlined in [33], however, in the continuum model for the twisted bilayer, as in the monolayer, the diamagnetic contribution vanishes due to the fact that ∂²H/∂k²_α = 0. Hence, we have simply the paramagnetic contribution, given by the Kubo formula

$$\begin{eqnarray} &&\fl \chi_{j_\alpha j_\beta}(\boldsymbol{q},\omega)=g_{\rm s} g_{\rm v}\int\frac{{\rm d}^2\boldsymbol{k}}{(2\pi)^2} \sum_{mn} \frac{n_{\rm F}(\epsilon_{m,\boldsymbol{k}})-n_{\rm F}(\epsilon_{n,\boldsymbol{k}+\boldsymbol{q}})}{\hbar\omega-\epsilon_{n,\boldsymbol{k}+\boldsymbol{q}}+\epsilon_{m,\boldsymbol{k}}+\mathrm{i}\delta}\langle m,\boldsymbol{k}|j_\alpha|n,\boldsymbol{k}+\boldsymbol{q}\rangle \langle n,\boldsymbol{k}+\boldsymbol{q}|j_\beta|m,\boldsymbol{k}\rangle.\nonumber\\ \end{eqnarray} \tag{ 7 }$$

Here, g_s = g_v = 2 are the spin and valley degeneracies. States |m,k〉 are eigenstates of H in subband m and of momentum k in the first Brillouin zone of the superstructure. Their eigenenergies are _mk and n_F is the Fermi function. The paramagnetic current operator is j_α = −∂H/∂k_α. In the continuum model with local interlayer coupling, this current does not mix layers and is thus independent of the coupling γ_⊥. For a system with time-reversal symmetry, the conductivity tensor is diagonal and due to the threefold rotational symmetry of (twisted) bilayer graphene proportional to the unity tensor. We will, therefore, drop the tensor indices in what follows.

The real part of the conductivity σ(ω) represents the power dissipation of the system upon incidence of light at frequency ω, while its imaginary part represents inductive retardation of the current with respect to the driving field. From equation (6), it is obvious that the conductivity is ill-defined at ω = 0 and δ → 0 if χ_jj(q = 0,ω = 0) ≠ 0. It is thus customary to split up the conductivity into a regular part and a term containing the delta-singularity

$$\begin{equation} \sigma(\omega)=\pi D\delta(\omega)+\sigma_{\mathrm{reg}}(\omega), \end{equation} \tag{ 8 }$$

$$\begin{equation} \sigma_{\mathrm{reg}}(\omega)=\frac{\mathrm{i}e^2}{\omega}\chi_{jj}(\boldsymbol{q}=0,\omega) \end{equation} \tag{ 9 }$$

with the Drude weight (or charge stiffness) D defined as

$$\begin{equation} D=e^2 \lim_{\omega \to0}\mathrm{Re}\,\chi_{jj}(\boldsymbol{q}=0,\omega)=\lim_{\omega\to 0} \omega \,\mathrm{Im}\,\sigma(\omega)\;. \end{equation} \tag{ 10 }$$

Above we used the fact that the imaginary part of the current response is an odd function. We thus have Im σ = Im σ_reg and will drop the subindex when confusion is unlikely.

3.1. Graphene monolayer

In the case of a graphene monolayer [34, 35], described within a continuum Dirac model, we have the famous zero temperature result (μ is the chemical potential with respect to the neutrality point)

$$\begin{eqnarray} &&\begin{array}{ll} \displaystyle \mathrm{Re}\, \sigma_{\mathrm{reg}}(\gamma_0\gg\omega>2\mu)=\sigma_0=\frac{\pi}{2}\frac{{e}^2}{h},\\ \displaystyle \mathrm{Re}\, \sigma_{\mathrm{reg}}(\omega<2\mu)=0. \end{array} \end{eqnarray} \tag{ 11 }$$

This universal value σ₀ for the Dirac gas includes valley and spin degeneracy factors, and is valid for energies ℏω ∼ 1 eV well below the γ₀ ∼ 3 eV van Hove singularity of the monolayer. As ω approaches this scale, a discrete model must be employed to compute the monolayer σ₀ (see [4] for general results).

The imaginary part of the conductivity, in the same limit of low frequencies, reduces to [3, 4]

$$\begin{equation} \mathrm{Im}\,\sigma(\omega\ll \gamma_0)=\frac{4}{\pi}\frac{\mu}{\hbar\omega}\sigma_0+\frac{1}{\pi}\log\left|\frac{\hbar\omega-2\mu}{\hbar\omega+2\mu}\right|\sigma_0, \end{equation} \tag{ 12 }$$

where μ is the Fermi energy, measured from the Dirac point (charge neutrality). Using equation (10), we have a Drude weight

$$\begin{equation} D_0=\frac{4\mu}{\pi\hbar}\sigma_0, \end{equation} \tag{ 13 }$$

which vanishes at neutrality. Note that the above expression only scales with the chemical potential and no other material constant (Fermi velocity) is involved. The Drude weight for twisted bilayer as function of μ will thus be independent of the twist angle in the regime where the electronic spectrum can be described by a (renormalized) Dirac model. This breaks down close to the first magic twist angle with the emergence of a flat band.

3.2. Twisted bilayer

In the more general case of a twisted bilayer, the computation of σ(ω) becomes a numerical task and it is usually simpler to compute Re σ_reg than Im σ, since, by virtue of equation (7), the former reduces to an integral constrained to eigenstates within ℏω of the Fermi energy

$$\begin{eqnarray} &&\fl \mathrm{Re}\, \sigma_{\mathrm{reg}}(\omega)=\frac{16\pi\sigma_0}{\omega}\int \!\!\frac{{\rm d}^2\boldsymbol{k}}{(2\pi)^2} \sum_{mn}\left[n_{\rm F}(\epsilon_{n,\boldsymbol{k}})-n_{\rm F}(\epsilon_{m,\boldsymbol{k}})\right]\left|\langle m,\boldsymbol{k}|j_x|n,\boldsymbol{k}\rangle\right|^2 \delta\left[\omega-(\epsilon_{m,\boldsymbol{k}}-\epsilon_{n,\boldsymbol{k}})/\hbar\right]\!.\nonumber\\ \end{eqnarray} \tag{ 14 }$$

The real part of the conductivity differs from that of two decoupled monolayers only at frequencies smaller than or comparable to the interlayer coupling ω < γ_⊥. Hence we have Re σ_reg(ω ≫ γ_⊥) ≈ 2σ₀. This result will be useful in the following discussion to obtain a regularized expression for the imaginary part of the conductivity.

To compute the imaginary part, and from that, the Drude weight D through equation (10), one could in principle employ equation (7). This is problematic, since it involves a summation, even for small ω, over the whole band, which in the continuum model extends to infinity. We will, therefore, make use of the Kramers–Kronig (KK) relation, applicable to any response function. In terms of the conductivity, this yields

$$\begin{equation} \mathrm{Im}\,\sigma(\omega)=\frac{2}{\pi\omega}\mathcal{P}\int_0^\infty {\rm d}\nu \frac{\nu^2\,\mathrm{Re}\,\sigma_{\mathrm{reg}}(\nu)}{\omega^2-\nu^2}, \end{equation} \tag{ 15 }$$

where $\mathcal {P}$ denotes the Cauchy principal value. There is apparently no advantage to using equation (15) as compared to equation (7), since the integral above extends to infinity, too. Moreover, the integral is ill defined, since at high frequencies the continuum model yields a constant Re σ_reg(ω ≫ γ_⊥) ≈ 2σ₀.

We will thus define the imaginary part relative to this divergent background following the usual procedure of ultraviolet field theories. In fact, this is one way of obtaining equation (17) and was discussed more generally in the context of the f -sum rule in [36]. The low energy properties are not modified by this regularization and we can simply subtract the ground-state expectation value of the two independent graphene layers at zero doping. Since the imaginary part for the neutral system is zero for the continuum model, the final regularized definition of Im σ simply reads

$$\begin{equation} \mathrm{Im}\,\sigma(\omega)=\frac{2}{\pi\omega}\int_{0}^\infty {\rm d}\nu\frac{\nu^2(\mathrm{Re}\,\sigma_{\mathrm{reg}}(\nu)-2\sigma_0)}{\omega^{2}-\nu^{2}}. \end{equation} \tag{ 16 }$$

Numerically, equation (16) can be evaluated by introducing a finite cutoff Λ for which Re σ^reg(Λ) ≈ 2σ₀. Integrating the second term involving the 'ground-state value' 2σ₀ then yields the final formula

$$\begin{equation} \mathrm{Im}\,\sigma(\omega)=\frac{2}{\pi\omega}\left(\mathcal{P}\int_0^\Lambda {\rm d}\nu \frac{\nu^2\mathrm{Re}\,\sigma_{\mathrm{reg}}(\nu)}{\omega^2-\nu^2}+2\sigma_0\Lambda\right)+\frac{2\sigma_0}{\pi}\log\left(\frac{\Lambda-\omega}{\Lambda+\omega}\right). \end{equation} \tag{ 17 }$$

Note that this expression should be independent of Λ as it was obtained from equation (16). The extra terms compared to the original KK-relation (15) can be interpreted as the diamagnetic contribution that naturally arises in the Dirac model with finite band cutoff.

Using this expression, the Drude weight D in equation (10) may be written as

$$\begin{equation} D=\frac{2}{\pi}\left(2\sigma_0\Lambda-\int_0^\Lambda {\rm d}\nu \mathrm{Re}\,\sigma_{\mathrm{reg}}(\nu)\right). \end{equation} \tag{ 18 }$$

This expression, like equation (17), becomes Λ-independent for large cutoffs. In terms of D, equation (17) is then finally rewritten as

$$\begin{equation} \mathrm{Im}\,\sigma(\omega)=\frac{D}{\omega}+\frac{2\omega}{\pi}\mathcal{P}\int_0^\infty {\rm d}\nu\frac{\mathrm{Re}\, \sigma_{\mathrm{reg}}(\nu)-2\sigma_0}{\omega^2-\nu^2}. \end{equation} \tag{ 19 }$$

From equation (18), we can also deduce the f -sum rule

$$\begin{equation} \int_0^\Lambda {\rm d}\omega \mathrm{Re}\,\sigma(\omega)=2\sigma_0\Lambda, \end{equation} \tag{ 20 }$$

which is independent of the angle θ, chemical potential μ or temperature T.

The conductivity at neutrality μ = 0 is shown in figure 3 for two 'large' twist angles. Much like the monolayer, the real part (black) exhibits a universal asymptotic value Re σ(ω → 0) = 2σ₀. It also exhibits signatures of the van Hove singularities, in form of a dip-peak structure as discussed in detail in [19]. The imaginary part (red) exhibits similar features, but vanishes at ω → 0, i.e. there is no Drude weight at μ = 0. This changes for finite temperature where the real part of σ is suppressed for $\hbar \omega \lesssim k_{\mathrm {B}}T$ [37], as a consequence of the cancelation of the n_F difference in equation (14). This results in a finite Drude weight, necessary to fulfill the f -sum rule equation (20). The Drude weight for finite temperature can be deduced from equation (13) with the substitution μ → 2 ln 2k_BT [38]. Note that increasing temperature merely changes the Drude weight (figure 7), not the width of the Drude peak, i.e. Re σ_reg(ω → 0) = 0 at finite T. Temperature can further cause a shift of spectral weight close to the transitions associated to the van Hove singularity (see the circled region), creating new peaks that are absent at zero temperature.

Zoom In Zoom Out Reset image size

In the case of finite doping μ ≠ 0, see figure 4, the imaginary part of the conductivity develops a  ∼ D(μ)/ω low frequency divergence, while the real part becomes suppressed below a threshold frequency 2μ for small enough μ (see black and red curves). Unlike for the monolayer, however, this is only true for μ below the van Hove singularity, since it relies on the linearity of the bandstructure (and DOS) around μ = 0. For chemical potentials close to the van Hove singularity, i.e. μ ≈ _M, the optical gap is again filled by spectral weight and the Drude weight reduces consistent with the f -sum rule (blue dashed curves). At finite temperature, the optical gap may persist for small twist angle and the Drude weight can thus even increase, see figure 7.

Zoom In Zoom Out Reset image size

For large chemical potentials μ ≫ _M, the electronic spectrum (DOS) of the continuum model is just that of two uncoupled graphene layers and one could thus expect to recover the conventional graphene conductivity. But the Moiré coupling between layers induces additional damping, since it opens up the possibility of more interband dissipation channels, i.e. any particle–hole excitation with momentum equal to a multiple of the momenta of the Moiré lattice. In order to understand the basic features of the conductivity on the basis of the uncoupled system, we can thus attribute a finite in-plane wavenumber q ∼ nΔK to the incoming photon (n > 0), see figure 5. The two-particle spectrum of interband and intraband excitations then predicts dissipation at low energies (black) and a broadened transition (blue) from the optical gap to the universal value 2σ₀ of the conductivity.

Zoom In Zoom Out Reset image size

In figure 6, the optical conductivity at the first magic angle i = 31 is shown for various chemical potentials (solid lines) and compared to Re σ_reg for i = 10 (dashed blue line). For low chemical potential, the optical response differs considerably for small and large angles, but for increasing μ the dependence on the twist angle becomes less and Re σ is characterized by a large Moiré-broadening of the optical gap around ℏω = 2μ. Still, there is the pronounced spectral weight at ℏω = _vH almost independent of μ, which is especially striking for larger angles. This feature cannot be explained from figure 5 and is twist angle dependent.

Zoom In Zoom Out Reset image size

Considering only intraband transitions of one unbounded conduction band, the Drude weight or charge stiffness is proportional to the chemical potential, independent of its dispersion relation. For parabolic band electrons, D is thus simply related to the charge density, D = e²n/m, whereas for Dirac Fermions one gets $D\sim \sqrt {n}$, see equation (13). In systems with a finite band-width, we expect a non-monotonous behavior expressing the fact that the Fermi energy is first related to electrons and later to holes when sweeping the chemical potential through the band. For the one-dimensional tight-binding model, one finds, e.g. $D\sim \sqrt {1-(\mu /2\gamma _0)^2}$ with 2γ₀ the band width. For the two-dimensional (2D) hexagonal tight-binding model, the system is characterized by a van Hove singularity at μ = γ₀ around which the Drude weight shows a maximum [37].

For twisted bilayer graphene described by the continuum model, we expect a combination of both features, because of the unbounded linear dispersion of Dirac Fermions and the superlattice structure of the Moiré-pattern. From figure 2 we see that the DOS is characterized by a van Hove singularity followed by a 'real' step-like singularity which is repeated to higher doping levels. This pseudo-gap structure divides the spectrum in several bands and will be followed by the Drude weight. Nevertheless, these modulations are superposed onto a constant linear background as suggested by the unbounded continuum model.

The μ dependence of the Drude weight D(μ) in units of 2D₀ is presented in figure 7 for two large twist angles. As is apparent, for small chemical potential D(μ) follows twice the monolayer value 2D₀ independent of the twist angle, as expected for Dirac systems. For larger μ, the Drude weight drops and remains less than 2D₀ until it converges to the uncoupled bilayer system for μ ≫ γ₀.

Zoom In Zoom Out Reset image size

As outlined above, D(μ) exhibits dips that accurately correlate with dips in the DOS, also shown in figure 7 (red). The resulting picture suggests the above mentioned approximate shell structure, albeit (in general) without gaps between shells, just as the suppressed DOS. The first shell can be fairly well approximated be the hexagonal tight-binding model.

At finite temperature T = 300 K (and also for disorder), the Drude weight is smeared out for large twist angle i = 10, thus reaffirming the interpretation of the repeated shell-structure. Nevertheless, the shell structure becomes less pronounced with smaller θ which is also reflected by its temperature response.

For small angles, the Drude weight provides a sensitive fingerprint for the breakdown of Dirac physics that occurs around the magical angles due to the emergence of flat bands. This can be seen on the right hand side of figure 8 for small angles above (i = 26) and below (i = 36) the critical angle. The Drude weight becomes larger than the Drude weight of two monolayers before it drops indicating the band width of the flat band. At the critical angle (i = 31), this band width becomes zero, and interestingly, the scaling of D₀ is recovered. Note that for i = 26, we find a transport gap D = 0 at finite chemical potential μ = 0.06γ₀.

Zoom In Zoom Out Reset image size

The increase of the band width for twist angles larger than θ^m₁ can be nicely appreciated from the Drude weight. This becomes considerably more difficult from the DOS, displayed on the left hand side of figure 8. The reason for this is that the Drude weight is obtained from integrating over the whole spectral region and not only from the low energy band structure. Nevertheless, the weight of the peak at zero energy is maximal at the magic angle (red curve), indicating its singularity.

A complementary aspect to the optical absorption is the plasmon excitation spectrum of the electron gas. Conceptually, a plasmon is an oscillating charge density mode (hence accompanied by an oscillating electric potential), which is sustained by the Coulomb interactions between electrons. The dispersion relation of these modes is given by the zeroes of the dielectric function (q,ω) = 0, that relates the screened electric potential V_sc(r,t) to an externally applied potential V_ext(r,t),

$$\begin{eqnarray*} &&V_{\rm sc}(\boldsymbol{q},\omega)=\frac{V_{\mathrm{ext}}(\boldsymbol{q},\omega)}{\epsilon(\boldsymbol{q},\omega)}. \end{eqnarray*}$$

A plasmon is, therefore, a finite solution V_sc(q,ω) that requires no external driving, i.e. it is self-sustained. A Landau-damped plasmon corresponds to a peak in 1/(q,ω), and decays due to dissipation into particle–hole excitations.

The dielectric function (q,ω) may be computed within the random phase approximation (RPA). In the long-wavelength limit, we may thus relate _RPA directly to σ(ω) using equation (6),

$$\begin{equation} \epsilon_{\mathrm{RPA}}(\boldsymbol{q},\omega)\approx 1+\frac{\mathrm{i}}{2\varepsilon_0\epsilon}\frac{|\boldsymbol{q}|}{\omega}\sigma(\omega)=1+\mathrm{i}\frac{\pi\alpha_{\mathrm{g}}}{2\epsilon} \frac{v_{\rm F} |\boldsymbol{q}|}{\omega}\frac{\sigma(\omega)}{\sigma_0}, \end{equation} \tag{ 21 }$$

where α_g = e²/(4π₀ℏv_F) ≈ 2.2 is graphene's fine structure constant. In equation (21) we have neglected corrections due to the finite interlayer distance, and consider that interlayer and intralayer electron–electron interaction are of the same form  ∼ 1/q. In this approximation all the collective modes correspond to charge excitations in which the charge oscillates in phase in the two layers [39, 40].

To compute ω_p for small frequencies, we may expand Im σ in equation (19) by the Drude term, $\mathrm {Im}\, \sigma (\omega )\approx \frac {D}{\omega }$. This results in the dispersion of a classical 2D plasmon

$$\begin{equation} \omega_{\rm p}\approx \sqrt{\frac{D|\boldsymbol{q}|}{2\varepsilon_0\epsilon}}. \end{equation} \tag{ 22 }$$

For two decoupled monolayers, one just replaces D = 8μ/πℏσ₀ above, which is shown as dotted curved line in figures 9 and 10.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

6.1. Plasmons in twisted bilayer

In twisted graphene bilayer there may exist four types of plasmonic modes or resonances: (i) undamped (conventional) graphene plasmons for chemical potentials with μ ≪ _M; (ii) interband 'plasmons' present also for μ = 0; (iii) damped plasmons for large chemical potential μ ≫ _M; and (iv) transverse plasmons. Below, we will briefly outline the main properties of these different types.

6.1.1. Undamped plasmons

At low doping, there exists a low energy window with Re σ = 0 that gives rise to undamped plasmons governed by the Drude weight D. In this regime, the twisted bilayer can be described by linear Dirac Fermions and D(μ) follows the Drude weight of the uncoupled system 2D₀, as shown in figure 7. This energy window becomes smaller for decreasing twist angle since the van Hove singularity moves closer to the neutrality point and is only relevant for large angles.

6.1.2. Interband plasmons

In monolayer graphene, a peak in the loss function associated to π → π* transitions around the van Hove singularity and with linear dispersion was first predicted by DFT-studies [41], and later experimentally observed in suspended graphene by electron energy-loss spectroscopy [42]. Within the hexagonal tight-binding model and RPA, i.e. without including correlation or renormalization effects, no zero of the dielectric function is obtained [43]. The absorption peak would thus be merely due to interband transitions enhanced by a band-structure effect; nevertheless, in bi- or multilayer, ^RPA(q,ω) becomes zero around the M point [44]. In the twisted graphene bilayer, there exist several van Hove singularities such that interband plasmons (or absorption enhancement) are to be expected. These are especially dominate for large twist angle and can lead to a broad optical gap in the interband excitations where Im σ < 0.

6.1.3. Damped plasmons

In the regime of large chemical potential, μ ≫ _M, the conventional intraband plasmon is recovered, albeit Moiré-damped. This can be deduced from figure 6 which will be used to discuss the plasmon dispersion in terms of the loss function. The dispersion only depends slightly on the twist angle and becomes well-defined for large μ, extending into the Landau damped region just as for the monolayer.

6.1.4. Transverse plasmons

There is also the possibility for transverse or transverse electromagnetic (TE) plasmons modes given by the condition [45]

$$\begin{equation} 1-\frac{\omega \mathrm{i} \sigma(\omega)}{c^2\sqrt{q^2-\omega^2/c^2}}=0. \end{equation} \tag{ 23 }$$

For TE modes, the plasmon dispersion is closely pinned to the light cone which justifies our local approximation (q → 0) and there are (damped) solutions to equation (23) if Im σ < 0. Since they can only exist in twisted bilayer surrounded by a homogeneous dielectric background [46], we will not discuss them any further.

6.2. Loss function

We have analyzed the DOS, optical conductivity, Drude weight and plasmon spectrum in a twisted bilayer for various angles, doping levels and temperature, and compared them to the uncoupled system. We have found a universal low-frequency conductivity and peaks at finite frequency that correlate to van Hove singularities in the DOS. These latter features, however, become scrambled at decreasing angles. We have also shown that the Drude weight in twisted bilayers is usually below that of two decoupled monolayers, 2D₀, and exhibits a shell-like structure, following again van Hove singularities in the DOS. Only for small angles close to the neutrality point, the Drude weight becomes larger than 2D₀, signaling the emergence of a flat band. The Drude weight resembles a sensitive quantity to numerically distinguish the different electronic structure at critical angles since it is obtained from integrating over the whole spectrum and is directly related to transport experiments. We have finally discussed how plasmons may become damped by the interlayer Moiré coupling which opens up new dissipation channels. Significant spectral structure survives, however, including undamped plasmons and multiple long-lifetime resonances (interband plasmons) within the extended interband particle–hole continuum.

We thank G Gómez-Santos and E Prada for discussions. This work has been supported by FCT under grants PTDC/FIS/101434/2008; PTDC/FIS/113199/2009, by MIC-Spain under grant FIS2010-21883-C02-02, FIS2011-23713 and FIS2012-33521, and by the European Research Council Advanced grant, contract 290846.