Plasmons in realistic graphene/hexagonal boron nitride moiré patterns

van der Waals heterostructures employing graphene and hexagonal boron nitride (hBN) crystals have emerged as a promising platform for plasmonics thanks to the tunability of their collective modes with carrier density and record values for plasmonics figures of merit. In this paper we investigate theoretically the role of moiré-pattern superlattices in nearly aligned graphene on hBN by using a continuum-model Hamiltonian derived from ab initio calculations. We calculate the system's energy-loss function for a range of chemical potential values that are accessible in gated devices. Our calculations reveal that the electron-hole asymmetry of the moiré bands leads to a remarkable asymmetry of the plasmon dispersion between positive and negative chemical potentials, showcasing the intricate band structure and rich absorption spectrum across the secondary Dirac point gap for the hole bands.

Figure 1

Graphene-hBN superlattice minibands along the $\tilde{\mathrm{\Gamma }}\text{−}\tilde{K}\text{−}{\tilde{M}}^{\prime }\text{−}{\tilde{K}}^{\prime }\text{−}\tilde{\mathrm{\Gamma }}$ direction in the mBZ, for the set of parameters given in Eq. (9), corresponding to vanishing twist angle $\theta =0$. On the horizontal axis, the quantity $\xi$ indicates the total length along the path in the reciprocal space. The Fermi surfaces at the extrema of the shaded area are shown in Fig. 2.

Figure 2

Fermi surfaces at the extrema of the shaded area shown in Fig. 1, for two Fermi energies ${\varepsilon }_{\mathrm{F}}=-180\mathrm{meV}$ (red) and $-215\mathrm{meV}$ (blue). The gray shaded area corresponds to the moiré superlattice Brillouin zone. The Fermi surfaces are periodically repeated in the reciprocal space for clarity.

Figure 3

A 2D density plot of the RPA loss function $L(\mathit{q},\omega )$, as a function of the chemical potential $\mu$ and excitation energy $\hslash \omega$, for the parameters given in Eq. (9). The calculation is performed for a fixed wave vector of length $q=0.02{\mathrm{nm}}^{-1}$ along the $\tilde{M}\text{−}\tilde{\mathrm{\Gamma }}$ direction. Panel (b) shows a magnification of (a) around the gap at $\mu \simeq -200\mathrm{meV}$ (cf. the minibands in Fig. 1). In (a) the red dotted line corresponds to Eq. (16) and the red dashed line to Eq. (17). In both panels, and in the following figures as well, the range of the monochrome shades has been truncated to improve the visibility of the less intense features. For comparison, the graph of the loss function at fixed chemical potential is shown in Fig. 6.

Figure 4

A 2D density plot of the RPA loss function $L(\mathit{q},\omega )$, as a function of the wave vector length $q$ and the excitation energy $\hslash \omega$. The wave vector is taken along the $\tilde{M}\text{−}\tilde{\mathrm{\Gamma }}$ direction. Results for positive $\mu =200\mathrm{meV}$ and negative $\mu =-205\mathrm{meV}$ chemical potentials are shown in (a) and (b), respectively.

Figure 5

As in Fig. 4, but in a smaller wave-vector range around the origin of the reciprocal space. Results for two negative chemical potentials $\mu =-205\mathrm{meV}$ and $-215\mathrm{meV}$, slightly lower than the gap (cf. the minibands in Fig. 1), are shown in (a) and (b), respectively.

Figure 6

(a) The RPA loss function $L(\mathit{q},\omega )$ (shaded area), the real (solid line), and the imaginary (dashed line) part of the quantity ${\epsilon }_{\mathit{0},\mathit{0}}(\mathit{q},\omega )$, as functions of the transition energy $\hslash \omega$, at fixed wave-vector length $q=0.01{\mathrm{nm}}^{-1}$. The values of the loss (dielectric) function are reported on the left (right) vertical axis. (b) A 2D density plot of the imaginary part of the dielectric function, as a function of the wave-vector length $q$ and the excitation energy $\hslash \omega$. In both panels, the wave vector is taken along the $\tilde{M}-\tilde{\mathrm{\Gamma }}$ direction and the chemical potential is $\mu =-215\mathrm{meV}$, as in Fig. 4.