Quantum dot-like plasmonic modes in twisted bilayer graphene supercells

We aim to study finite-size twisted bilayer graphene supercells as depicted in figure 1. These are constructed such that for zero twist angle ($\theta = 0^{\circ}$) we get an AB (Bernal) stacking. The rotation axis is centered at the upper A sublattice (lower B sublattice), as indicated in the inset of figure 1. The outer boundaries are chosen to be of armchair type. We describe these supercells with a low-energy Hamiltonian for the p_z states,

$$\begin{align} H = \sum_{i,j,\sigma} t_{ij} c_{i\sigma}^\dagger c_{j\sigma} + \sum_{i} U_{ii} n_{i\uparrow} n_{i\downarrow} + \sum_{i>j,\sigma, \sigma^{^{\prime}}} U_{ij} n_{i\sigma} n_{j\sigma^{^{\prime}}}, \end{align} \tag{ 1 }$$

with i and j labeling atomic site positions of (twisted) bilayer graphene. $c_{i\sigma}$ ($c_{i\sigma}^\dagger$) and $n_{i\sigma} = c_{i\sigma}^\dagger c_{i\sigma}$ are p_z -orbital annihilation (creation) and orbital occupation number operators, respectively, at site i with spin σ. t_ij and U_ij are hopping and density-density Coulomb interaction matrix elements, respectively, which we derive in the following section from ab initio calculations. At this stage we do not explicitly differentiate between the upper and lower layer: i and j label lattice sites in both layers.

Zoom In Zoom Out Reset image size

To study plasmonic properties we utilize a real-space version of the RPA [38–42] to calculate the electron energy loss spectra (EELS). We start with the eigen decomposition of $\varepsilon(\omega)$:

$$\begin{equation*} \varepsilon(\omega) = \sum_n \varepsilon_n(\omega) \left|\phi_n(\omega)\rangle \right. \left.\langle \phi_n(\omega)\right|, \end{equation*}$$

where we order the eigenvalues ε_n according to $-\operatorname{Im}[1 / \varepsilon_n(\omega)]$. ε₁ is the eigenvalue with the largest $-\operatorname{Im}[1 / \varepsilon_n(\omega)]$ and we call it 'leading'. EELS is defined by:

$$\begin{equation*} \operatorname{EELS}(\omega) = -\operatorname{Im}\left[ \frac{1}{\varepsilon_1(\omega)} \right]. \end{equation*}$$

$\phi_1(r, \omega) = \langle r | \phi_1(\omega) \rangle$ is the corresponding eigenvector in real space which renders the plasmonic excitation pattern. Within a real-space tight-binding approximation, the RPA dielectric matrix is given by:

$$\begin{equation} \varepsilon_{ij}(\omega) = \delta_{ij} - \sum_{k} U_{ik} \Pi_{kj}(\omega), \end{equation} \tag{ 2 }$$

with the Coulomb interaction U_ik entering the Hamiltonian (1) and the polarizability function $\Pi_{ij}$ given by:

$$\begin{equation} \Pi_{ij}(\omega) = 2~\cdot \sum_{ab} \psi_{ia}^* \psi_{ib} \psi_{jb}^* \psi_{ja} \frac{f_a - f_b}{E_a - E_b + \omega + i\eta}. \end{equation} \tag{ 3 }$$

Here E_a , $\psi_{ia} \equiv \psi_{a}(i)$ and $f_a = f(E_a)$ are eigenvalues, eigenvectors and corresponding Fermi–Dirac distribution functions obtained upon diagonalization of the single-particle tight-binding Hamiltonian. i, j, and k label site indices while a and b label eigenstates of the Hamiltonian. η is a small positive constant which we set to $1\,$meV in our calculations. More details on the implementation can be found in [41]. In appendix C we furthermore describe how these real-space RPA calculations can be significantly accelerated by exploiting the sparsity of the involved matrices and making use of modern graphical processing units (GPUs). In the following we derive all the necessary model parameters from ab initio.

In order to derive all model parameters for arbitrary twist angles via down folding of first principles calculations we have to make one well justified approximation: only the p_z states will experience the moiré potential and thus the twist angle θ. To stress the validity of this assumption we show in figure 2(a) the sp²-projected density functional theory (DFT) DOS of pristine AA-stacked (blue) and AB-stacked (orange) bilayer graphene. Their difference is nearly invisible so that we can readily assume that the relative layer rotation does not affect these electronic states. This, of course, does not hold for the p_z states, as is obvious from p_z -projected DOS shown in figure 2(b).

Zoom In Zoom Out Reset image size

Based on this assumption we can derive the single-particle hopping matrix elements t_ij for the p_z states via a Wannier construction based on DFT calculations for the AB-stacked bilayer graphene (see appendix B for details). We calculate the hopping matrix elements via:

$$\begin{equation*} t_{ij}(\theta = 0^{\circ}) = \langle w_i | H_\textrm{DFT}^{AB} | w_j \rangle, \end{equation*}$$

for the untwisted ($\theta = 0^{\circ}$) geometry and using the DFT Hamiltonian together with the p_z -like ab initio Wannier functions $w_i(r)$. In table 1 we list the resulting intralayer and interlayer hopping matrix elements for an interlayer distance of $d = 3.35$Å. To account for finite twisting angles on the single-particle level, we utilize a Slater-Koster based interlayer hopping model [43] $t_\perp(r) = \gamma_0~\operatorname{exp}[-\alpha (r-r_0)]$, which we fit to the interlayer hopping matrix elements from table 1 to obtain $\gamma_0 = 0.29\,$eV and $\alpha = 5.63~{\unicode{x00C5}}^{-1}$. As mentioned above, for small twist angles one would additionally need to account for modulations in the interlayer distance [35–37], but here we are interested in moderate to large twist angles such that the assumption of a purely nominal twisting is adequate.

The description of fully screened, retarded and θ-dependent Coulomb interaction $W(\omega, \theta)$ requires more attention. $W(\omega, \theta)$ is defined by:

$$\begin{align} W(\omega, \theta) = \frac{\upsilon}{1 - \upsilon \Pi^\textrm{total}(\omega, \theta)}, \end{align} \tag{ 4 }$$

where v is the bare Coulomb interaction. $\Pi^\textrm{total}(\omega, \theta)$ renders all possible screening processes at a given rotation angle θ which can be separated into two terms:

$$\begin{align} \Pi^\textrm{total}(\omega, \theta) \approx \Pi^{p_z}(\omega, \theta) + \Pi^{\textrm{rest}}(\omega = 0), \end{align} \tag{ 5 }$$

with $\Pi^{p_z}(\omega, \theta)$ being the partial polarizability resulting from virtual excitations within the low-energy p_z sub-space. It can be evaluated using equation (3). The rest polarizability $\Pi^{\textrm{rest}}(\omega = 0)$ describes in principle instantaneous screening processes from virtual excitations from, to, and between non-p_z states (such as sp² and others). Due to the orthogonality of p_z and sp² states, the cross-polarization terms between p_z and non-p_z orbitals can be safely neglected [44]. This way, $\Pi^{\textrm{rest}}$ effectively describes the screening of a semiconductor with a gap of about $10\,$eV (see figure 2(a)), which renders the instantaneous approximation valid. Importantly, this renders the background (or rest) polarizability independent of the twist angle. Equation (4) now reads:

$$\begin{align} W(\omega, \theta) & = \frac{\upsilon}{1 - \upsilon \left[\Pi^{p_z}(\omega, \theta) + \Pi^\textrm{rest}(\omega = 0)\right]}\nonumber\\ & = \frac{U}{1 - U \Pi^{p_z}(\omega, \theta)}\,, \end{align} \tag{ 6 }$$

$$\begin{align} \textrm{with} \qquad U & = \frac{\upsilon}{1 - \upsilon \Pi^\textrm{rest}(\omega = 0)}, \end{align} \tag{ 7 }$$

being the θ-independent background-screened Coulomb interaction, as needed for the evaluation of equation (2). We calculate U within the constrained RPA [45] based on ab initio calculations for AB-stacked bilayer graphene (see appendix B for details). This yields discretized U_ij with i, j labeling AB bilayer graphene lattice positions. To check for the θ independence of U we compared the local and nearest-neighbor (in and out-of-plane) Coulomb interactions U_ij between AB and AA stacked bilayer finding modifications smaller than $3\%$. For the evaluation of equation (2) we, however, need to evaluate U_ij at other positions resulting from the finite twisting angle θ. To this end, we map the discretized U_ij to a continuous model $U(r = r_i - r_j)$. For the latter we choose the analytic image-charge model for the potential within a dielectric slab of height d reading [42, 46–48]:

$$\begin{equation} U(r) = \mathop{\underbrace{\frac{e^2}{\varepsilon_m z_0(r)}}}\limits_{\mathop{\textrm{Local}}\limits_{\textrm{screening}}} + \mathop{\underbrace{2~\sum_{n = 1}^\infty \frac{e^2~\beta_b^n}{\varepsilon_m z_n(r)}}}\limits_{\mathop{\textrm{Non-local}}\limits_\textrm{screening}}, \end{equation} \tag{ 8 }$$

with e being the elementary charge, ε_m the dielectric constant of the slab, $z_n(r) = \sqrt{r^2 + \delta^2 + (nh)^2}$, and $\beta_b = (\varepsilon_m - 1) / (\varepsilon_m + 1)$. An additional parameter δ allows us to also fit the on-site potential $U_{ii} = U(r = 0)$. In figure 3 we show the ab initio cRPA data together with the fit using equation (8), the locally-screened interaction $\frac{e^2}{\varepsilon_m z_0(r)}$ ($h \rightarrow \infty$), and the fully-screened interaction $W(\omega = 0,\theta = 0)$. For the fit we fixed $d = 6.7\,\unicode{x00C5}$ (twice the interlayer distance) and find $\varepsilon_m \approx 2.26$ (in good agreement with similar fits in momentum space [49]) and $\delta \approx 0.763\,\unicode{x00C5}$, which evidently interpolate the ab initio data well. From the comparison to the locally-screened interaction we see that the background screening, as described by the second term in equation (8), acts differently at each r due to its non-local character. The fully-screened interaction $W(\omega = 0,\theta = 0)$ behaves as expected from Thomas-Fermi screening theory in two dimensions [50] which predicts a strongly decaying potential with a r⁻³ asymptotic behavior, but in our calculations we also find a finite offset c. This offset c decays with the supercell size and vanishes in the infinite-size limit (not shown). We attribute this behavior to finite-size/boundary effects. In detail, although the polarization function $\Pi^{p^z}(r, r^{^{\prime}}, \omega = 0)$ is rather localized, as shown in figure 3, it still has some non-vanishing oscillating tails due to the finite Fermi surface. These tails in r are partially missing if $r^{^{\prime}}$ is fixed to an edge side, which likely induces the finite offset c.

Zoom In Zoom Out Reset image size

Equipped with the continuous representation of the background-screened Coulomb interaction U(r) and the tight-binding model described by the hopping matrix elements from table 1 we can evaluate the dielectric function from equation (2) and thus plasmonic properties for arbitrary twist angles θ. At this point it is important to note the limits of our material-realistic Coulomb modeling approach. In the 'low-energy' range ($\omega < 5\,$eV) we are mostly describing virtual excitations (i.e. screening processes) solely within the p_z manifold. For larger excitation energies transitions involving the rest space would become important, which we do not fully render here. Thus, we restrict our plasmon considerations to energies below $5\,$eV.

The momentum resolved EELS($q, \omega$) of pristine doped bilayer graphene is dominated by two plasmonic branches, $\omega_+(q)$ and $\omega_-(q)$, sketched in figure 4(a). The $\omega_+(q)$ branch is higher in energy with $\omega_+(q) \propto \sqrt{q}$ in the long-wavelength limit, while the $\omega_-(q)$ branch is linear in this limit, $\omega_-(q) \propto q$ [51–53]. The upper branch does not show any layer polarization with charge modulations only within the bilayer plane. Classically the corresponding plasmon can thus be understood as interacting charge monopoles governed by conventional monopole Coulomb interaction terms so that the usual $\sqrt{q}$ dispersion is expected. The lower branch describes plasmons with a layer polarization, which hence forms interacting charge dipoles in the bilayer plane. In contrast to the monopole interactions, the resulting dipole interactions yield a linear plasmonic dispersion. These long-wavelength-limit dispersion relations are indicated as dashed lines in the sketch of figure 4(a). For intermediate q we see that these trends are modified. The $\sqrt{q}$ dispersion of the $\omega_+(q)$ branch is significantly reduced (flattened) due to the non-local background screening from the sp² states (and, in reality, also due to substrate screening) [42, 54]. The $\omega_-(q)$ is also affected, but rather mildly. For pristine infinite systems, the resulting full spectral functions of these plasmonic branches are indicated in figure 4(b). They are continuous functions attaining their maxima before entering the inter-band continuum.

Zoom In Zoom Out Reset image size

In the finite-size bilayer graphene supercells under consideration here, these plasmonic properties are modified as sketched in figure 4(c). The finite diameter of the hexagonal supercells defines an upper plasmonic wavelength λ_max and a minimal wave vector $q_\textrm{min} \propto \lambda_\textrm{max}^{-1}$. Below this $q_\textrm{min}$ the system cannot host any of the plasmonic excitations of the $\omega_{\pm}(q)$ branches. Furthermore, the finite system size also discretizes the allowed momenta in multiples of this $q_\textrm{min}$. As a result the plasmonic dispersions become gapped between these q values and the full plasmonic spectrum becomes discretized as well, as sketched in figures 4(c) and (d). The previously continuous plasmonic spectrum now shows well defined discretized levels.

This is in full analogy to the formation of electronic quantum dot states: the finite system size discretizes the continuous electronic spectrum into well defined electronic quantum dot states. Based on this analogy, we call the resulting plasmonic states in the finite-size cells plasmonic quantum dot states. Their tunability is, however, much richer than in the case of electronic quantum dots. The energies of the plasmonic quantum dot states are controlled by (a) the system size as this affects $q_\textrm{min}$, (b) the doping level and substrate screening as these control the energies of the underlying continuous branches $\omega_{\pm}(q)$, as well as by (c) the twisting angle, as we will discuss in more detail in the following. Also the discretization of the linear $\omega_-(q)$ branch yields a very different level spacing than the discretized $\omega_+(q)$ branch as indicated in figure 4(d).

In the following sections we discuss how this scheme is realized in full real-space RPA calculations for finite AB-stacked (and twisted) bilayer graphene supercells.

We start our discussion with comparing EELS($\omega,q$) for different finite-size supercells of untwisted AB-stacked bilayer graphene, which is defined by:

$$\begin{equation*} \operatorname{EELS}(\omega, q) = -\operatorname{Im}\left[ \frac{1}{\varepsilon(\omega, q)} \right], \end{equation*}$$

where $\varepsilon(\omega, q)$ is the diagonal of the formal Fourier transform of $\varepsilon_{ij}(\omega)$. In figures 5(a) and (b) we show EELS($\omega,q$) for hexagonal supercells of diameter $d \approx 150~\unicode{x00C5}$ (11 028 atoms) and $d\approx80~\unicode{x00C5}$ (3252 atoms), respectively. In both cases the electron doping is about $n = 5.3\times10^{14}\,$cm$^{-2}\,$ (µ = 1.34 eV), which is around the maximum of what is achievable with double sided ionic-liquid gating techniques [55]. In both of these panels we can clearly identify the $\omega_\pm(q)$ branches. In the case of the larger supercell we find a nearly continuous $\omega_+(q)$ branch, which starts to deviate from the $\sqrt{q}$ form around $q\approx 0.05\,\unicode{x00C5}^{-1}$ showing the characteristic flattening. The lower $\omega_-(q)$ branch shows its approximate linear dispersion. For $q >0.2\,\unicode{x00C5}^{-1}$ (above $\omega \approx 2.0\,$eV) both branches smear out due to strong Landau damping induced by the particle-hole continuum. Also note that the approximate plasmonic energies are similar to the ones observed in highly doped graphene monolayers [56]. Hence, the overall characteristics of these plasmonic branches are fully in line with the expectations as discussed in section 3.1, underlining the appropriateness of our modeling scheme. For the larger supercell only the lower branch shows first signs of discretization, which becomes more evident for the smaller cell in figure 5(b). In this case we find clear gaps in both branches, which is also reflected in the discrete peaks in the full EELS(ω) depicted in figure 5(c).

Zoom In Zoom Out Reset image size

From the comparison between the momentum resolved EELS($\omega,q$) with the plasmonic levels in the full EELS(ω), we can identify the origin of the specific plasmonic states. Starting at the lowest level with $\omega \approx 0.32\,$eV (barely visible in figure 5(c)) we see that this mode arises from the layer polarized $\omega_-(q)$ branch (black doted line in figures 5(b) and (c)). The corresponding layer-resolved real-space representation of that mode (the eigenvector $\phi_1(r, \omega\approx 0.32$ eV) of the dielectric function) is shown in figure 6(a). As expected, we find a layer-polarized excitation with a homogeneous charge distribution within each layer. The next plasmonic level at $\omega \approx 0.63\,$eV is again resulting from the $\omega_-(q)$ branch, which accordingly shows the layer polarization, see figure 6(c). However, the charge distribution within each layer is not homogeneous anymore and this mode is best described as a layer-polarized dipole excitation. It is accompanied at the same q by a plasmonic excitation with $\omega \approx 1.2\,$eV, which originates from the upper $\omega_+(q)$ branch (blue dashed lines in figures 5(b) and (c)). Accordingly, this mode does not show a layer polarization (figure 6(c)), but is otherwise also of dipole character. The plasmonic mode at $\omega \approx 1.0\,$eV is again originating from the layer polarized branch, as is also evident from its real-space representation depicted in figure 6(e). This mode has a circular charge distribution, which is vastly reminiscent of a 1s atomic wave function. This mode is as well accompanied at the same q by a 1s-like plasmonic state at $\omega \approx 1.65\,$eV from the upper branch (orange dashed-doted line in figures 5(b) and (c)), which does not show any layer polarization. Finally, we find another layer-polarized plasmonic excitation around $\omega \approx 1.3\,$eV which hosts a p-like state (see figure 6(f)) which is again reminiscent of the atomic-like wave functions of electronic quantum dots. For this mode, we were however not able to identify the accompanying mode from the upper plasmonic branch, which we attribute to smaller energy separation of plasmonic states at higher energies, as also indicated in the sketches in figure 4.

Zoom In Zoom Out Reset image size

Thus, next to the finite-size induced discrete energy levels and tunability, these plasmonic charge distributions can also resemble atomic-like wave functions of electronic quantum dots, which further underlines their identification as plasmonic quantum dots. In the following we proceed with the twist-angle dependence, focusing on the dipole and 1s-like modes.

From the previous discussion, the dependence of the plasmonic quantum dot energies on doping level and system size are already clear. In the following we thus focus on the twist angle dependence. To this end we present in figure 7 the 'dark' and 'bright' dipole and 1s modes for $\theta = 0^{\circ}, \,10^{\circ}, \,20^{\circ}, \,30^{\circ}$. Modes are called 'dark' or 'bright' when they originate from $\omega_-(q)$ or $\omega_+(q)$ branches, respectively.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For the bright dipole modes, shown in figure 7(b), we observe that the boundary separating the differently charged areas is rotating when we adjust θ, synchronously in the lower and upper layer, however, only with $\theta/2$. The latter can be readily understood by overlaying the rotated dipole patterns and by remembering that the charge distributions in the two layers are not independent. The missing overlap at the corners of the hexagons additionally yields enhanced charge accumulations in the two opposite corners of each layer. The dark dipole mode, shown in figure 7(a), behaves similarly. The charge separation line again rotates with $\theta/2$, but we observe a smearing of it, such that the charge separation is not as sharp as in the bright dipole mode. Note that rotated dark dipole modes are additionally rotated by 90^∘ compared to $\theta = 0^{\circ}$ case. This, however, does not affect any of the drawn conclusions as these unrotated dipole modes have nearly the same excitation energies, see figure 9 for more details.

Zoom In Zoom Out Reset image size

In figures 7(c) and (d) we depict patterns of the 1s quantum dot mode. Except for a slight deformation of the initial hexagonal shape we do not see any major changes to the bright excitation upon twisting. Its dark counterpart also does not show any significant changes in its excitation pattern, except for a smearing of the clear charge separation.

In fact, the most important changes to these plasmonic states are their excitation energy shifts. In figure 8 we plot the excitation energies for all modes as a function of θ. We see that within the given accuracy, the excitation energies of the bright modes do not dependent on the rotation angle. The dark modes, however, do. We see for both the dark dipole and the dark 1s mode that the energy is significantly decreased upon 10 degree twist. Afterwards, these modes just mildly dependent on further rotation towards $\theta = 30^{\circ}$. The initial symmetry breaking between 0^∘ and 10^∘ thus leaves the strongest footprint in the plasmonic energies, while the larger rotation angles do not change it too drastically anymore.

We attribute this different impact of the rotation to the bright and dark modes to the single-particle properties. In figure 4, the original continuous 'bright' plasmonic branch $\omega_+(q)$ is clearly detached from the particle-hole continuum (shown in grey), while the 'dark' $\omega_-(q)$ branch is in its close vicinity. Thus, the dark states are much more prone to Landau damping and to changes in the particle-hole continuum due to modifications of the the single-particle electronic states.

To verify this hypothesis we calculated EELS(ω) from

$$\begin{equation*} \varepsilon(\omega, \theta, \alpha) = 1 - U(\alpha) \Pi^{p_z}(\omega, \theta), \end{equation*}$$

where θ and α are both referring to the twist angle, but now separated in the rotation in the polarizability (θ) and the Coulomb interaction (α) matrices. For all previous calculations we have set $\alpha = \theta$. From the resulting spectra in figure 11 we see that the rotation in the Coulomb interaction matrix, i.e. changes to α, barely affect the energy of the dark 1s mode, whereas the rotation in the polarizability matrix (changes to θ), strongly affect it. Thus, indeed the single-particle properties as reflected by the polarizability affect the dark states the most. If we additionally remember that the approximate morié lattice constant λ is inversely proportional to $\sin(\theta/2)$, yielding $\lambda(\theta = 10^{\circ}) \approx 14\,\unicode{x00C5}$, $\lambda(\theta = 20^{\circ}) \approx 8\,\unicode{x00C5}$, and $\lambda(\theta = 30^{\circ}) \approx 6\,\unicode{x00C5}$, we understand that the changes to the lattice structure are the largest when going from $\theta = 0^{\circ}$ to $\theta = 10^{\circ}$, while increasing θ further modifies λ only slightly. Under the reasonable assumption that changes to the lattice structure are also imprinted on the electronic structure, we can qualitatively understand why the excitation energy of the dark modes is affected the most upon the initial twist by $\theta = 10^{\circ}$.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size