First-principles calculations and model analysis of plasmon excitations in graphene

TL;DR

This study uses first-principles calculations to analyze plasmon excitations in graphene and graphene/hBN heterostructures, revealing how substrates and anisotropy influence plasmon behavior and proposing a new dispersion model.

Contribution

It introduces a unified theoretical framework for analyzing plasmon dispersion and lifetime in graphene systems, including substrate effects and a novel dispersion model for $oldsymbol{ ext{pi}}$ plasmons.

Findings

The dispersion of $ ext{pi}$ plasmons follows $oldsymbol{ ext{ω}}_ ext{π}(q) = ext{√}(E_g^2 + eta q)$ at small $q$.

Substrate and anisotropic effects significantly influence plasmon properties.

The proposed dispersion model aligns well with first-principles calculations.

Abstract

Plasmon excitations in free-standing graphene and graphene/hexagonal boron nitride (hBN) heterostructure are studied using linear-response time-dependent density functional theory within the random phase approximation. Within a single theoretical framework, we examine both the plasmon dispersion behavior and lifetime (line width) of Dirac and $\pi$ plasmons on an equal footing. Particular attention is paid to the influence of the hBN substrate and the anisotropic effect. Furthermore, a model-based analysis indicates that the correct dispersion behavior of $\pi$ plasmons should be $\omega_\pi(q) = \sqrt{E_g^2 + \beta q}$ for small $q$'s, where $E_g$ is the band gap at the $M$ point in the Brillouin zone, and $\beta$ is a fitting parameter. This model is radically different from previous proposals, but in good agreement with our calculated results from first principles.

Equations

$$-\text{Im}\{{\epsilon }^{-1}(\text{q},\omega )\}=-\frac{4\pi }{q^2}\text{Im}\{{\chi }_{\text{G}=0,{\text{G}}^{\prime }=0}(\text{q},\omega )\}.$$

$${\chi }_{\text{G},{\text{G}}^{\prime }}(\text{q},\omega )={\chi }_{\text{G},{\text{G}}^{\prime }}^0(\text{q},\omega )+\sum _{{\text{G}}_1,{\text{G}}_2}{\chi }_{\text{G},{\text{G}}_1}^0(\text{q},\omega )K_{{\text{G}}_1,{\text{G}}_2}(\text{q},\omega ){\chi }_{{\text{G}}_2,{\text{G}}^{\prime }}(\text{q},\omega ),$$

$${\chi }_{\text{G},{\text{G}}^{\prime }}^0(\text{q},\omega )=\frac{1}{\Omega }\sum _{\text{k}}^{\text{1BZ}}\sum _{n,n^{\prime }}\frac{f_{n,\text{k}}-f_{n^{\prime },\text{k}+\text{q}}}{\omega +{ϵ}_{n,\text{k}}-{ϵ}_{n^{\prime },\text{k}+\text{q}}+i\eta }⟨n,\text{k}|e^{-i(\text{q}+\text{G})\text{r}}|n^{\prime },\text{k}+\text{q}⟩⟨n^{\prime },\text{k}+\text{q}|e^{i(\text{q}+{\text{G}}^{\prime }){\text{r}}^{\prime }}|n,\text{k}⟩,$$

$$K_{{\text{G}}_1,{\text{G}}_2}^C(\text{q})=\frac{4\pi }{|\text{q}+{\text{G}}_1{|}^2}{\delta }_{{\text{G}}_1,{\text{G}}_2}.$$

$${\tilde{K}}_{{\text{G}}_1,{\text{G}}_2}^C(\text{q})=\frac{4\pi {\delta }_{{\mathbf{G}}_1,{\mathbf{G}}_2}}{|\mathbf{q}+{\mathbf{G}}_1{|}^2}\left[1+e^{-|\mathbf{q}+{\bar{\mathbf{G}}}_1|R}\left(\frac{G_{1,z}}{|\mathbf{q}+{\bar{\mathbf{G}}}_1|}\text{sin}(G_{1,z}R)-\text{cos}(G_{1,z}R)\right)\right]$$

$${\tilde{K}}_{{\text{G}}_1,{\text{G}}_2}^C(\bar{\mathbf{q}})=\frac{4\pi {\delta }_{{\text{G}}_1,{\text{G}}_2}}{|\bar{\mathbf{q}}+{\text{G}}_1{|}^2}\left[1-(-1{)}^{n_z}{e}^{-|\bar{\mathbf{q}}+{\bar{\text{G}}}_1|\frac{L_z}{2}}\right]$$

$$-\text{Im}\{{ϵ}_{\text{2D}}^{-1}(\bar{\mathbf{q}},\omega )\}=-\frac{4\pi (1-e^{-\bar{q}{L}_z/2})}{{\bar{q}}^2}\text{Im}\{{\tilde{\chi }}_{\text{G}=0,{\text{G}}^{\prime }=0}(\bar{\mathbf{q}},\omega )\},$$

$${\chi }_{\bar{\mathbf{G}},{\bar{\mathbf{G}}}^{\prime }}(z,z^{\prime },\bar{\mathbf{q}},\omega )={\chi }_{\bar{\mathbf{G}},{\bar{\mathbf{G}}}^{\prime }}^0(z,z^{\prime },\bar{\mathbf{q}},\omega )+\sum _{{\bar{\mathbf{G}}}_1,{\bar{\mathbf{G}}}_2}\iint d{z}_{1}d{z}_2{\chi }_{\bar{\mathbf{G}},{\bar{\mathbf{G}}}_1}^0(z,z_1,\bar{\mathbf{q}},\omega )K_{{\bar{\mathbf{G}}}_1,{\bar{\mathbf{G}}}_2}(z_1,z_2,\bar{\mathbf{q}}){\chi }_{{\bar{\mathbf{G}}}_2,{\bar{\mathbf{G}}}^{\prime }}(z_2,z^{\prime },\bar{\mathbf{q}},\omega ),$$

$$K_{{\bar{\mathbf{G}}}_1,{\bar{\mathbf{G}}}_2}^C(z_1,z_2,\bar{\mathbf{q}})=\frac{2\pi {\delta }_{{\bar{\mathbf{G}}}_1,{\bar{\mathbf{G}}}_2}}{|\bar{\mathbf{q}}+{\bar{\mathbf{G}}}_1|}{e}^{-|\bar{\mathbf{q}}+{\bar{\mathbf{G}}}_1||z_1-z_2|}.$$

$${\chi }_{\bar{\mathbf{G}},{\bar{\mathbf{G}}}^{\prime }}^0(z,z^{\prime },\bar{\mathbf{q}},\omega )=\sum _{G_z,G_z^{\prime }}{e}^{-i{G}_{z}z}{\chi }_{\mathbf{G},{\mathbf{G}}^{\prime }}^0(\bar{\mathbf{q}})e^{i{G}_z^{\prime }{z}^{\prime }}.$$

$${\chi }_{\text{G},{\text{G}}^{\prime }}^S(\text{q},\omega )=\frac{1}{\Omega }\sum _{\text{k}}^{\text{1BZ}}\sum _{n,n^{\prime }}(f_{n,\text{k}}-f_{n^{\prime },\text{k}+\text{q}})\delta (\omega +{ϵ}_{n,\text{k}}-{ϵ}_{n^{\prime },\text{k}+\text{q}})⟨n,\text{k}|e^{-i(\text{q}+\text{G})\text{r}}|n^{\prime },\text{k}+\text{q}⟩⟨n^{\prime },\text{k}+\text{q}|e^{i(\text{q}+{\text{G}}^{\prime }){\text{r}}^{\prime }}|n,\text{k}⟩,$$

$${\chi }_{\text{G},{\text{G}}^{\prime }}^0(\text{q},\omega )={\int }_0^{\infty }d{\omega }^{\prime }[\frac{1}{\omega -{\omega }^{\prime }+i\eta }-\frac{1}{\omega +{\omega }^{\prime }+i\eta }]{\chi }_{\text{G},{\text{G}}^{\prime }}^S(\text{q},{\omega }^{\prime }).$$

$$\delta(\omega_{i}-\Delta)\approx\left\{\begin{array}[]{ll}\frac{\Delta-\omega_{i-1}}{\omega_{i}-\omega_{i-1}},&\text{for}~{}~{}\omega_{i-1}<\Delta<\omega_{i}\\ &\\ \frac{\omega_{i+1}-\Delta}{\omega_{i+1}-\omega_{i}},&\text{for}~{}~{}\omega_{i}<\Delta<\omega_{i+1}\\ \end{array}\right.$$

$$H=t\sum _{⟨i,j⟩}{\hat{c}}_i^{\ast }{\hat{c}}_j$$

$${ϵ}^{c,v}(k_x,k_y)=\pm t{\left[1+4\text{cos}\left(\frac{3a{k}_y}{2}\right)\text{cos}\left(\frac{\sqrt{3}a{k}_x}{2}\right)+4{\text{cos}}^2\left(\frac{\sqrt{3}a{k}_x}{2}\right)\right]}^{1/2},$$

$$Re\{\epsilon (\mathbf{q},\omega )\}=1-\frac{2\pi }{q}\sum _{\mathbf{k}}\frac{|⟨v,\mathbf{k}|e^{-i\mathbf{q}\cdot \mathbf{r}}|c,\mathbf{k}+\mathbf{q}⟩{|}^{2}2({ϵ}_{c,\mathbf{k}+\mathbf{q}}-{ϵ}_{v,\mathbf{k}})}{{\omega }^2-({ϵ}_{c,\mathbf{k}+\mathbf{q}}-{ϵ}_{v,\mathbf{k}}{)}^2},$$

$$E_g({\mathbf{k}}^{\prime },\mathbf{q})\approx$$

$$\approx$$

$$=$$

$${\text{lim}}_{q\to 0}⟨v,\mathbf{k}|e^{-i\mathbf{q}\cdot \mathbf{r}}|c,\mathbf{k}+\mathbf{q}⟩=$$

$$=$$

$$Re\{\epsilon (\mathbf{q},\omega )\}$$

$$\approx \frac{q}{{\omega }^2-E_{g,M}^2}\frac{2\pi {\hslash }^4}{m_e^2}\sum _{|{\mathbf{k}}^{\prime }|<=k_{cut}}\frac{p^2({\mathbf{k}}^{\prime })2E_g({\mathbf{k}}^{\prime },\mathbf{q})}{\left[1-\frac{2E_{g,M}{E}_g({\mathbf{k}}^{\prime },\mathbf{q})}{{\omega }^2-E_{g,M}^2}\right]E_g^2({\mathbf{k}}^{\prime },0)}$$

$$=1-\frac{q\beta (\mathbf{q},\omega )}{{\omega }^2-E_{g,M}^2}.$$

$$1-\frac{\beta q}{{\omega }^2-E_{g,M}^2}=0$$

$$\omega (q)=\sqrt{E_{g,M}^2+\beta q}.$$

Full text

I Introduction

In recent years the plasmonic excitations in graphene and graphene-related materials have attracted considerable attention both theoretically Huang et al. (1997); Shyu and Lin (2000); Wunsch et al. (2006); Hwang and Das Sarma (2007); Jablan et al. (2009); Yan et al. (2011) and experimentally Eberlein et al. (2008); Kramberger et al. (2008); Liu et al. (2008); Lu et al. (2009); Koppens et al. (2011); Kinyanjui et al. (2012); Grigorenko et al. (2012); Woessner et al. (2015); Liou et al. (2015) due to their importance for basic physics and technological applications. Derived from the unique electronic band structure, the plasmon excitation spectra of graphene span a wide energy range, and fall into three distinct regimes. At low energies (0-2 eV) with finite electron doping, graphene can sustain the so-called Dirac plasmons, originating from the intraband transitions of Dirac fermions in the vicinity of $K$ points of the Brillouin zone (BZ) Wunsch et al. (2006); Hwang and Das Sarma (2007). At higher energies (4-15 eV), there exist the intrinsic $\pi$ plasmons in graphene due to the collective excitations of electrons from $\pi$ to $\pi^{\ast}$ bandsMarinopoulos et al. (2004). The $\pi$ plasmons are very dispersive, starting at $\sim 4$ eV for small momentum transfer $q$ all the way to $\sim 15$ eV for the approaching to the boundary of the BZ. At even higher energies, $\sigma$ bands start to contribute, the mixture of the $\pi\rightarrow\pi^{\ast}$ and $\sigma$ transitions leads to yet another set of distinct plasmon peaks, usually denoted as $\pi$+$\sigma$ plasmons. As such, the rich plasmon physics in graphene renders the system distinguished from normal metals, conventional 2-dimensional electron gas, and doped semiconductors.

The low-energy Dirac plasmon excitations are only present under finite electron or hole dopings Wunsch et al. (2006); Hwang and Das Sarma (2007), which can be readily achieved by chemical means or by electric gating.This type of plasmons attract enormous interest for potential technological applications due to the consideration given to both strong field localization and low energy loss, simultaneouslyJablan et al. (2009); Low and Avouris (2014). Furthermore, they are highly flexible in the sense that their frequencies can be tuned from Terahertz to middle infrared by varying the doping level, and can be “engineered” by encapsulating graphene in between other 2D layered materials Woessner et al. (2015). Because of this unique property, graphene has been considered as a promising material for fabricating nano-plasmonic devices. A review on the recent progress and future perspective of this field can be found in Ref. [Low and Avouris, 2014]. The $\pi$ and $\pi$+$\sigma$ plasmons at higher energies are also of significant scientific interest and have been under intensive theoretical Huang et al. (1997); Yan et al. (2011); Kramberger et al. (2008) and experimental Lu et al. (2009); Koppens et al. (2011) investigations. It should be noted that these latter types of plasmons are also present in the parent material of graphene – graphite. However,it shows that the dispersion behavior in graphene is much different from the graphite’s, especially in small $q$. Consequently, the interest here is often on the dispersion behavior of the plasmon peaks Kramberger et al. (2008); Lu et al. (2009); Liou et al. (2015) as a function of the momentum transfer $q$, as well as the influence of the substrates Yan et al. (2011) and/or the interlayer interactions within graphene-based heterostructures or multi-layer graphenes. Koppens et al. (2011).

Historically, researches on the above-mentioned different plasmon types have been largely conducted within different communities. Theoretically, the studies of Dirac plasmons have been predominantly carried out with the linear band model (also known as “Dirac-cone approximation”) Wunsch et al. (2006); Hwang and Das Sarma (2007); Jablan et al. (2009), which is valid only in the vicinity of the $K$ points in the BZ. Nevertheless, recent ab initio studies on Dirac plasmons based on time-dependent density-functional theory (TDDFT) revealed the appreciable non-isotropic effect (the variation of the plasmon dispersions along different directions in the Brillouin zone) Gao and Yuan (2011); Pisarra et al. (2014) and the existence of acoustic Dirac plasmon mode in graphene Pisarra et al. (2014). Physics of this type such cannot be captured by the linear band model. The studies of $\pi$ plasmons, on the other hand, have been mainly done with ab-initio TDDFT within the random-phase approximation (RPA) Kramberger et al. (2008); Yan et al. (2011); Liou et al. (2015), except for some earlier ones based on tight-binding model involving only $\pi$ bands Huang et al. (1997); Shyu and Lin (2000). Compared to other two types of plasmons, $\pi$+$\sigma$ plasmons have received less attention, and the study of these excitations would require including high-lying energy bands. Hence for $\pi$+$\sigma$ plasmons, the ab-initio approach would be more appropriate. All together, it appears that TDDFT-RPA can offer a unified description of all types of plasmon excitations in graphene and graphene-related materials in a fully ab initio manner, thus eliminating possible ambiguities arising from the choice of model parameters and/or approximations that might miss important physical effects.

Despite the intensive studies of the plasmon excitations in graphene, some issues remain unclear and await for further investigations. For instance, the dispersion behavior of $\pi$ plasmons has long been thought to be linear Kramberger et al. (2008); Yan et al. (2011); Lu et al. (2009); Kinyanjui et al. (2012), but this view was challenged by Liou et al. Liou et al. (2015) recently. These authors claimed that $\pi$ plasmons follow a $\sqrt{q}$ behavior for small $q$’s, similar to the case of Dirac plasmons. In this work, we present a comprehensive study of Dirac and $\pi$ plasmon excitations in graphene and graphene/monolayer hexagonal boron nitride (hBN) heterostructure using ab initio TDDFT-RPA approach. We examine both the dispersion behavior (plasmon peak positions) and the line-width (plasmon lifetimes) as a function of the momentum transfer ${\bf q}$. Our ab initio calculations are further complemented by model analysis, which allows us to gain new insight into the dispersion behavior of $\pi$ plasmons. For the lifetimes, ab initio results from TDDFT-RPA have not been reported in the literature. In this connection, we note that only the Landau damping channel is adequately treated within the standard RPA. Effects from impurities, disorders, and electron-phonon couplings are usually interpreted via the framework of RPA combined with the number conserving relaxation time (RT) approximationMermin (1970); Jablan et al. (2009). In this approach, the relaxation time due to the scattering from impurities or phonons is estimated empirically and put into the RPA-RT framework by hand. It appears a full ab initio treatment of the plasmon lifetime problem taking into all important physical channels is still not practical. In this work, our ab initio lifetime calculations are restricted to the level of electron-electron scattering within RPA. A full account of other contributions from first principles will be pursued in future work.

The rest of the paper is organized as follows. In Sec. II the basic equations of the TDDFT-RPA approach behind our implementation are presented. This is followed by a description of the implementation details and computational setups in Sec. III. The calculated results for graphene and graphene/hBN are then presented in Sec. IV, complemented by an in-depth model analysis and discussions. Finally we summarize our work in Sec. V.

II Methods

The plasmon excitations can be measured experimentally by the momentum-resolved electron energy loss spectroscopy (EELS) technique. Theoretically and computationally, TDDFT-RPA represents a powerful ab initio approach to describe EELS, as thoroughly discussed in Ref. [Onida et al., 2002; Silkin et al., 2004; Yuan and Gao, 2009; Mowbray, 2014]. A most recent analysis of the role of theoretical dielectric function models in describing EELS can be found in Ref. [Azzolini et al., 2017]. Here we present the key equations of the formalism as a basis for the pertinent technical details behind our implementations. In brief, EELS can be obtained from the inverse of the dielectric function $\varepsilon$, or equivalently the linear charge density response function $\chi$,

[TABLE]

where ${\bf G}$ and ${\bf G^{\prime}}$ are the 3-dimensional (3D) reciprocal lattice vectors, ${\bf q}$ is a wave vector within the first BZ (1BZ), and $\omega$ is the frequency. Here only the “head” term (${\bf G}={\bf G^{\prime}}=0$) of $\varepsilon^{-1}$ is of our concern, which corresponds to the response function of the system at the macroscopic scale. Within TDDFT, the system’s interacting response function $\chi$ is linked to its non-interacting counterpart $\chi_{0}$ via the Dyson equation:

[TABLE]

where the non-interacting response function $\chi_{0}$ is given explicitly by the well-known Adler-Wiser formula Adler (1962); Wiser (1963)

[TABLE]

In Eq. (3), $\Omega$ stands for the volume of the Born-von-Karmen supercell, $f_{n,k}$, $\epsilon_{n,k}$, $|n,{\bf k}\rangle$ are the Fermi occupation numbers, Kohn-Sham (KS) Kohn and Sham (1965) eigenvalues, and KS eigenvectors, respectively. The computation of $\chi^{0}_{\textbf{G},\textbf{G}^{\prime}}(\textbf{q},\omega)$ will be discussed in the next section. Within RPA, the full kernel $K_{\textbf{G}_{1},\textbf{G}_{2}}(\textbf{q},\omega)$ in Eq (2) is reduced to the static Coulomb kernel,

[TABLE]

The above formalism (Eqs. 2-4) is perfectly suitable for 3D periodic bulk materials, but needs modifications for 2-dimensional (2D) materials as explained below. In the 2D case, the system is infinite and periodic only in the basal ($x$-$y$) plane, but confined in the third ($z$) direction. A so-called supercell approach is commonly used to treat 2D systems, in which the system is modeled by repeated 2D slabs, separated by a large vacuum region in the $z$ direction. The advantage of the supercell approach is that the 3D formalism presented above can still be used, but care must be taken to remove spurious interactions between the periodic replicas. The issue is already well-known in semi-local DFT calculations for polar surfaces. The situation is more severe in the present case, because of the explicit presence of long-range Coulomb interaction in Eq. (2).

In the literature, two schemes have been employed to deal with the above-mentioned problem. The first scheme, introduced by Rozzi et al. Rozzi et al. (2006), sticks to the reciprocal-space formalism, but employs the Fourier transform of a truncated Coulomb potential, instead of the bare one (Eq. 4). Specifically, by restricting the Coulomb interaction within a window $[-R,R]$ in the $z$ direction, the Coulomb kernel in Eq. (4) becomes

[TABLE]

where ${\bf G}_{1}=(\bar{{\bf G}}_{1},G_{1,z})$ and ${\bf q}=(\bar{{\bf q}},0)$, with $\bar{{\bf G}}_{1}$ and $\bar{{\bf q}}$ being respectively the two-dimensional reciprocal lattice vector and Bloch wavevector in the basal plane. It has been suggested in Ref. [Rozzi et al., 2006] to choose $R=L_{z}/2$ where $L_{z}$ is the length of the lattice vector in the $z$ direction, and hence $G_{z}=2\pi n_{z}/L_{z}$ with $n_{z}$ being an integer number. With this choice, the truncated Coulomb kernel simplifies to,

[TABLE]

Now, by replacing the full Coulomb kernel by the truncated one $\tilde{K}^{C}_{\textbf{G}_{1},\textbf{G}_{2}}(\bar{{\bf q}})$ in Eqs. (1), one obtains the modified 3D interacting response function $\tilde{\chi}$. Finally the loss function of a genuine 2D material is computed as,

[TABLE]

where $\bar{q}=|\bar{{\bf q}}|$. It is easy to see that the Fourier transform of the truncated Coulomb interaction approaches its 2D form $2\pi/\bar{q}$ for $\bar{q}<<1/L_{z}$. For plasmon excitations in graphene, it has been shown recently by Mowbray Mowbray (2014) that the artificial interactions between periodic replicas can affect the plasmon peak positions and intensities quite significantly. Such effects can be efficiently removed by employing the truncated Coulomb kernel (6). Please note that, for consistency, in Eq. (7) we also used the truncated Coulomb interaction form $4\pi(1-e^{-\bar{q}L_{z}/2}/\bar{q}^{2}$ for the prefactor. Different from our choice, however, in literature the full Coulomb potential $\frac{4\pi}{{\bar{q}}^{2}}$ is often used for the prefactor in Eq. (7). Our numerical tests show that, using the truncated Coulomb potential here instead of the full one for the prefactor will not produce noticeable difference in the loss spectrum except for extremely small $\bar{q}$’s. Even in the small $\bar{q}$ regime, only the height of the plasmon peaks is slightly affected, but not the peak positions.

An alternative approach to deal with the artificial interaction issue, as employed by Silkin and coauthors Silkin et al. (2004); Bergara et al. (2003) and also by Yuan and Gao Yuan and Gao (2009); Gao and Yuan (2011), is a mixed representation of the TDDFT-RPA equations. In this approach, a reciprocal-space representation for the in-plane ($x$-$y$) dimensions, and a real-space representation for the out-of-plane ($z$) dimension are employed. In such a mixed representation, the Dyson equation for the response function (Eq. 2) becomes,

[TABLE]

where the Coulomb kernel is given by

[TABLE]

In this way, one is essentially dealing with an isolated 2D system, and the non-physical interactions between periodic images in the 3D reciprocal space formalism are naturally absent. In practice, however, the mixed representation of the non-interacting response function is often obtained from its 3D reciprocal-space form,

[TABLE]

Consequently, $\chi^{0}$ becomes de facto periodic in the $z$ direction of periodicity length $L_{z}$. Because of this feature, the integration over $z$ and $z^{\prime}$ in Eq. 8 has to be restricted within $(-L_{z}/2,L_{z}/2)$. The obtained $\chi_{\bar{{\bf G}},\bar{{\bf G}}^{\prime}}(z,z^{\prime})$ can then be Fourier transformed to its 3D reciprocal-space form $\chi_{{\bf G},{\bf G^{\prime}}}$, and the loss spectrum is again computed using Eq. (7). As such, the above two schemes to remove artificial interactions between periodic images become essentially equivalent Pisarra et al. (2016). In the present work, we directly follow the first scheme (Eqs. 2,3,6,7) in our implementation.

III Implementation Details

One efficient way to compute $\chi^{0}$ is to first calculate its imaginary part $\chi^{S}_{\textbf{G},\textbf{G}^{\prime}}(\textbf{q},\omega)$,

[TABLE]

and then obtain the full $\chi^{0}$ via the Hilbert transform,

[TABLE]

Following Ref. [Shishkin and Kresse, 2006], the $\delta$-function in Eq. (11) is approximated by a triangular function,

[TABLE]

where $\Delta=\epsilon_{n^{\prime},\textbf{k}+\textbf{q}}-\epsilon_{n,\textbf{k}}$, and $\omega_{i}$ is the frequency grid point. Furthermore, To obtain converged values for the peak positions and width of the loss spectrum, technical parameters such as the ${\bf k}$-point mesh, the frequency grid points, and the actual values for the positive infinitesimal parameter $\eta$ must be chosen carefully.

The linear-response TDDFT-RPA equations presented above have been implemented in a recently released first-principles code package Atomic-orbital Based Ab-initio Computations at UStc (ABACUS) Chen et al. (2010); Li et al. (2016); aba . The Troullier-Martins Troullier and Martins (1991) norm-conserving pseudopotential in its fully separable form Kleinman and Bylander (1982) is used to describe the interactions between nuclear ions and valence electrons. A “dual basis set” strategy is adopted in ABACUS which allows to use both plane waves and numerical atomic orbitals (NAOs) bases. In this work, the plane-wave basis sets are employed to expand the valence-electron wave functions in Eqs. (3) and (11). The Perdew-Zunger local-density approximation (LDA) Perdew and Zunger (1981) is used for the preceding Kohn-Sham (KS)-DFT Kohn and Sham (1965) calculations to obtain the KS orbitals and orbital energies. Unless otherwise stated, the production calculations in this work are done with the plane-wave basis with a cutoff energy of 50 Rydberg (Ry).

As mentioned in Sec. II, the supercell approach is used here to model the 2D systems. The distance ($L_{z}$ in Eq. 5) between periodic layers (slabs) is chosen to be 20 Å. This choice, together with the truncated Coulomb potential technique discussed in Sec. II, ensures a clean removal of the artificial interactions between periodic images in the $z$ direction. In this work, we study two types of systems: (a) a single freestanding graphene and a hBN layer, and (b) a combined graphene/hBN double-layer heterostructure in AA stacking. For the former, 20 bands (4 occupied plus 16 unoccupied) are included in the calculations of $\chi^{S}$ in Eq. (11); for the latter, 30 bands (8 occupied plus 22 unoccupied) are used. These choices guarantees a convergence of the plasmon spectra up to 30 eV, thus covering the entire energy range of $\pi+\sigma$ plasmons.

In the plasmon dispersion calculations, the positive infinitesimal parameter $\eta$ in Eq. (11) is set to be 0.01Ry, and the $\chi$ matrix is expanded in terms of 100 ${\bf G}$-vectors. The Brillouin zone is sampled with a 192$\times$192$\times$1 Monkhorst-Pack ${\bf k}$-point grid. Furthermore, we use a uniform frequency grid with energy spacing of 0.001 Ry up to 4 Ry, containing 4000 grid points. Extensive tests show that such a combination of parameters are sufficient to obtain converged results of the dispersion behavior of the plasmon excitations. For lifetimes, highly accurate numbers can only be obtained by extrapolating $\eta\rightarrow 0$. Benchmark tests on the convergence behavior of lifetime calculations will be shown in Appendix C.

In this work, we study the plasmon behavior of extrinsic graphene systems with finite electron doping. To this end, the Fermi level is shifted upwards by 0.05Ry(0.68eV) above the Dirac point, and this corresponds to a doping level of 0.027 electrons per unit cell (or free charge carrier concentration of $5.1\times 10^{13}$ cm*-2*) and the resultant Fermi vector is $k_{F}$ = 0.127 Å*-1* as determined from the LDA band structure. In doing so, the tiny effect of doping on the graphene electronic structure is neglected. Such a procedure is also followed by Pisarra et al. in Refs. [Pisarra et al., 2014, 2016]. Please note that the Fermi velocity is slightly underestimated within LDA, and correspondingly the the Fermi vector is slightly overestimated for a given doping. But this does not affect the discussions in the present work.

IV Results

In this section, the computed results of the full plasmon spectra of graphene and graphene/hBN will be presented. We will discuss both the dispersion relations and lifetimes. The influence of the hBN substrate on plasmons in graphene will be also examined.

IV.1 The plasmon dispersion behavior

IV.1.1 Overall features

In Fig 1 the loss spectra computed using the linear-response TDDFT-RPA approach are presented both for free-standing graphene (upper panel, Fig. 1(a)) and the mixed graphene/hBN bilayer (in AA stacking) system (lower panel, Fig. 1(c)). The result for a free-standing hBN layer is also shown (middle panel, Fig. 1(b)) for comparison. The different spectral lines aligned vertically correspond to different momentum transfer ${\bf q}$ along the $\Gamma-M$ direction in the 1BZ. For calculations of graphene and graphene/hBN, the Fermi level is shifted upwards by 0.05 Ry above the Dirac point to mimic the effect of finite electron doping. The results for graphene/hBN presented in this section are for AA stacking. However, as shown in Appendix A, other types of stacking yield essentially the same results.

The obtained loss functions of doped graphene and graphene/hBN clearly show three distinct plasmon modes, corresponding respectively to the Dirac plasmons, $\pi$ plasmons, and $\pi$+$\sigma$ plasmons from low to high excitation energies. The Dirac plasmon peaks are very sharp and pronounced, while $\pi$ plasmon peaks are much broader and carry more spectral weights. The peaks of $\pi$+$\sigma$ plasmons are even broader, and multiple sub-peak structures within this regime are clearly visible. As mentioned in the Introduction, the low-energy Dirac plasmon excitations are only present for extrinsic graphene with finite doping, but the $\pi$ and $\pi$+$\sigma$ are intrinsic and can be activated at both finite and zero dopings.

The comparison between Fig. 1(a) and (c) clearly reveals the major effect of a single-layer hBN substrate on the plasmon dispersion behavior of graphene. While the Dirac plasmons are not much affected, $\pi$ plasmons are significantly changed. An overall binodal peak structure, running parallel to each other as a function of ${\bf q}$, can be recognized. Comparing the graphene/hBN spectra (Fig. 1(c)) to the individual graphene (Fig. 1(a)) and hBN (Fig. 1(b)) ones, further suggests that the low-energy peaks stem from graphene, though substantially suppressed, while the high-energy peaks originate from hBN, staying essentially the same as their counterparts in the isolated hBN layer. In this context, we note that hBN, a wide-gap insulator, sustains plasmon excitations which behave similarly to the $\pi$ plasmons in graphene, though blue-shifted to higher energies, as can be seen from Fig. 1(b). The suppression of $\pi$ plasmons of graphene upon including a substrate was also observed in graphene/SiC, as reported in Ref. [Yan et al., 2011].

One may notice that, at small ${\bf q}$’s, there are finer plasmon peak structures within the $\pi$-region of the graphene/hBN heterostructure. To unravel the origin of these sub-peaks, we plot in Fig. 2 the band structure of graphene/hBN, obtained with DFT-LDA. Also shown is the zoom-in of the plasmon peaks at $q=0.015$Å*-1* within an energy window of $3$ - $7$ eV. It can be seen that the band structure of hBN is actually very similar to graphene around the $M$ point in BZ, with the direct gap (at $M$) of hBN about 2 eV larger than that of graphene. A close inspection of the two graphs further reveals that each small individual peak can be associated with an inter-band transition between the two occupied bands and two unoccupied bands of the hybrid graphene/hBN system around the $M$ point, as labeled in Fig. 2. For example, the first peak (cf. the inset of Fig. 2) with an energy of 3.71 eV corresponds well to the transition from the 8-th to 9-th band, while the third peak of 5.50 eV corresponds to the $8\rightarrow 10$ interband transition, and so on.

The similarity between the band structures of graphene and hBN around the $M$ point, and the similarity of their $\pi$ plasmons are by no means a coincidence. Actually, similar to what happens in graphite Marinopoulos et al. (2004), the $\pi$ plasmons stem from the collective electronic interband excitations around the $M$ point, which is a saddle point in the BZ, and contributes most strongly to the density of states (DOS). Furthermore, in contrast to 3D materials, the excitation energy approaches to the value of single-particle energy gap as $q\rightarrow 0$ for plasmon excitations in 2D materials originating from collective interband transition with a finite gap. This explains the blue shift of the $\pi$ plasmons in hBN compared to that in graphene. Below we will examine separately the dispersion relations of Dirac plasmons and $\pi$ plasmons in more details.

IV.1.2 Dirac plasmons

In Fig. 3(a), the spectra of Dirac plasmons along the $\Gamma-K$ direction in doped graphene (upper panel) and graphene/hBN (lower panel) are presented. Compared to the spectra along the $\Gamma-M$ direction, there is an additional shoulder peak along $\Gamma-K$, besides the “conventional” 2D Dirac plasmons, appearing at even lower energies. This additional peak structure has been interpreted by Pisarra et al. Pisarra et al. (2014) as an acoustic plasmon mode. Specifically, plasmon excitations along $\Gamma-K$ arise from oscillations of charge carriers with two different Fermi velocities. These two types of charge carrier oscillations can interact with each other, leading to an acoustic mode and an optical mode Pines (1956), with the latter being the conventional Dirac plasmons. Such an effect is not present for plasmons along the $\Gamma-M$ direction. Comparing the spectra of graphene and graphene/hBN, we see that little has changed for both acoustic and conventional Dirac plasmons, upon including the hBN substrate.

The plasmon peak positions as a function of the momentum transfer $q$ (the dispersion relations) are plotted in Fig.  3(b) for Dirac plasmons along both the $\Gamma-M$ and $\Gamma-K$ directions. For the acoustic plasmon mode along the $\Gamma-K$ direction, the spectral weight gets vanishingly small as $q\rightarrow 0$, making a precise determination of peak positions difficult. Hence the dispersion relation of the acoustic plasmons was not shown for very small $q$’s in Fig. 3(b). For “conventional” Dirac plasmons with $q<k_{\text{F}}$($\sim$ 0.12 Å*-1*), the dispersion roughly follows a $\sqrt{q}$ behavior, in agreement with what was originally found in model studies Hwang and Das Sarma (2007). For large $q$’s, where the Landau damping comes into play, the conventional Dirac plasmons actually display a quasi-linear dispersion behavior, as can be seen from Fig. 3(b). Furthermore, both the anisotropic and substrate effects can be seen in Fig 3(b). Namely, for large $q$’s, the plasmon energies are larger along the $\Gamma-K$ direction compared to the $\Gamma-M$ direction for the same $q$, i.e., the $\Gamma-M$ dispersion curve is below that of $\Gamma-K$. Also, a slight red shift can be observed for large $q$’s, when the hBN substrate is included.

IV.1.3 $\pi$ plasmons

Next we examine the dispersion behavior of the $\pi$ plasmons. In Fig 4, dispersion relations of $\pi$ plasmons in graphene and graphene/hBN are presented, along both $\Gamma-K$ and $\Gamma-M$ directions. For graphene/hBN, as discussed above, the plasmons in the $\pi$ region are split into several sub-peaks, and here only the majority sub-peak, arising mainly from hBN, is plotted. The dispersion behavior of graphene/hBN $\pi$ plasmons gets a bit complex for momentum transfer $q>=0.55$ Å*-1*, where kinks can be seen in the dispersion curve. This is due to the fact that the hybridization effect gets strong for large $q$’s, associated with substantial spectral weight transfer among the sub-peaks. Consequently, it is not possible any more to identity the original dominating (hBN-originated) peak. Also, for $q>0.5$ Å*-1*, the dispersion along $\Gamma-K$ and $\Gamma-M$ start to deviate from each other, and an anisotropy effect becomes visible.

In this connection, we would like to stress that the dispersion relation of $\pi$ plasmons in graphene is still a debated issue. Except perhaps for one work 111The dispersion behavior of $\pi$ plasmons for small ${\bf q}$’s was not explicitly discussed in Ref. [Huang et al., 1997], but the plot in Fig 4 of this paper suggests a parabolic behavior, as seen by eyes., most earlier theoretical and experimental studies Kramberger et al. (2008); Lu et al. (2009); Yan et al. (2011); Kinyanjui et al. (2012) reported a quasi-linear dispersion behavior for $\pi$ plasmons. This generally accepted view was recently challenged by Liou et al. Liou et al. (2015), who proposed a $\sqrt{q}$ dispersion relation based on their high-momentum-resolution EELS experiment. From Fig. 4, however, one can see that our first-principles TDDFT-RPA result cannot be fitted for by a simple linear or $\sqrt{q}$ model. This inconsistency motivated us to analyze this issue more deeply in terms of a tight-binding model involving only the $\pi$ and $\pi*$ bands of graphene. Our model analysis suggest that the theoretically more appropriate dispersion behavior of $\pi$ plasmons for small $q$’s should be $\sqrt{E_{g,M}^{2}+\beta q}$, where $E_{g,M}$ is the energy gap of graphene at the $M$ point. Further details of this derivation are given in Appendix B. The essence behind our derivation is to recognize that the $\pi$ plasmons in graphene stem from the collective excitations from the $\pi$ to $\pi^{\ast}$ band around the $M$ point. The $M$ point in the BZ represents a saddle point of the $\pi$/$\pi^{\ast}$ bands, yielding a van Hove singularity in the DOS at energy $E=E_{g,M}$. Thus the collective interband transition in the vicinity of the $M$ point dominates the contributions to $\pi$ plasmons at small $q$’s. For a 2D system like graphene, the plasmon excitation energy of this type approaches the single-particle energy gap at the $M$ point as $q$ goes to 0.

In Fig 5 we plotted the model dispersion curve of $\omega_{\pi}(q)=\sqrt{E_{g,M}^{2}+\beta q}$ with $\beta=90$  eV${}^{2}\cdot$Å, and the first-principles TDDFT-RPA dispersion curve along $\Gamma-M$, which was already shown in Fig. 4. The value of the parameter was chosen to fit well the slope of the first-principles results in the small ${\bf q}$’s region. Also presented for comparison are the experimental data by Liou et al. Liou et al. (2015) and Lu et al. Lu et al. (2009). It can be seen that, with one adjustable parameter, the model curve agrees very well with the first-principles one for $q<0.3$ Å*-1*, which is the regime that the model is expected to be valid. The overall agreement with the experimental data is satisfactory, though some remaining discrepancy is noticeable. One may notice that the experimental data from the two groups also show some scattering, and especially the data for 0.01 $<q<$ 0.1 Å*-1* are missing, which are important for an unambiguous determination for the $q\rightarrow 0$ behavior of the plasmon dispersion. We note that the argument of Liou et al. Liou et al. (2015) for the $\sqrt{q}$-dependence behavior of $\pi$ plasmons was based on a 2-dimensionally (2D) free-electron gas model for graphene. Such a model will certainly be valid for describing the Dirac plasmons in doped graphene. However, the $\pi$ plasmons stem from the collective excitations from the $\pi$ band to the $\pi^{\ast}$ around the $M$ point, which has an energy gap of $\sim$4 eV. From this standpoint, we argue that the dispersion behavior of $\pi$ plasmons should follow that of a 2D insulator, instead of a 2D metal. In this context, we call for more careful experiments to test our model.

IV.2 Lifetimes

From the perspective of technological applications, the lifetime of plasmons is an important quantity to care about. Theoretically the lifetime of plasmon excitations is inversely proportional to the width of the spectral peak, and hence the calculation of the lifetime is equivalent to determining the line width of the plasmon spectral peaks. Figure 1 clearly shows that Dirac plasmons have very sharp peaks for small $q$’s; these peaks however quickly broadens and diminishes as $q$ increases. This is consistent with the common understanding, as gained from model studies Hwang and Das Sarma (2007); Wunsch et al. (2006), that Dirac plasmons have infinite lifetimes within RPA for momentum transfers below certain cutoff $q_{c}$. (For the doping level considered in this work ($5.1\times 10^{13}$ cm*-2*), $q_{c}\sim 0.8k_{\text{F}}=0.1$ Å*-1*.) Above $q_{c}$, the Landau damping kicks in, due to the merging of plasmon excitations and individual particle-hole excitations, and consequently the Dirac plasmons start to gain a finite lifetime.

However, an accurate determination of the lifetimes numerically in an ab initio calculation is not entirely trivial. This is because the actual width of the peaks depends on the choice of the small positive parameter $\eta$ in Eqs. (3) and (12). By systematically reducing the value of $\eta$, one can in principle get more and more accurate peak widths, but care must be taken to to employ denser ${\bf k}$ and frequency grids to get smooth and converged spectral peaks. In this work, we adopt the following procedure to calculate the line width: first, the full widths at half maximum (FWHM) of the plasmon peaks are determined for several different $\eta$ values ranging from 0.01 Ry to 0.001 Ry, and for each $\eta$, the spectrum is converged with respect to the number of ${\bf k}$ and frequency points. Second, the FWHM’s at finite $\eta$’s are extrapolated to the limit of $\eta=0$. We consider the extrapolated FWHM as the final converged line-width at RPA level. Further details about our FWHM determination procedure is given in Appendix C. We successfully applied this procedure to accurately determine the FWHM’s of Dirac and $\pi$ plasmons of graphene. However, we admit that this procedure does not go without limitations. Especially, when there are several peaks very close to each other, an unambiguous determination of the width of each individual peak gets very difficult. This is the case, e.g., for the $\pi$ plasmons in graphene/hBN, and $\pi$+$\sigma$ plasmons.

In Fig 6 the extrapolated FWHM values for both Dirac and $\pi$ plasmons are presented as a function of the momentum transfer $q$. For Dirac plasmons (Fig 6(a)), the extrapolated FWHM value is strictly zero for $q<q_{c}$, indicating an infinite lifetime in this regime. For $q>q_{c}$, however, the FWHM of Dirac plasmons increases quadratically as $(q-q_{c})^{2}$. This behavior for FWHM’s of Dirac plasmons is in close agreement with the results of Wunsch et al. Wunsch et al. (2006), obtained from model studies. Interestingly, in this case, the hBN substrate seems to have the tendency of reducing the line-width, i.e., increasing the lifetimes of Dirac plasmons $q>q_{c}$. However, the effect here is purely electronic and different from the phonon and impurity scattering effects as discussed in Ref. [Woessner et al., 2015].

Distinct from Dirac plasmons, the $\pi$ and $\pi$+$\sigma$ plasmons have finite lifetimes for all momenta ${\bf q}$, because in these cases the Landau damping mechanism is always active. Here only the FWHM values for $\pi$ plasmons in freestanding graphene are presented, since the above described numerical procedure for determining FWHM does not yield reliable results for $\pi$ plasmons in graphene/hBN and for $\pi+\sigma$ plasmons. From Fig. 6(b), one can see that the FWHM value of $\pi$ plasmons first steadily increase as a function of $q$, reaching its maximum for $q\sim 0.12$ Å*-1*, and then, surprisingly, the line width starts to slowly decrease for even larger $q$’s. Such a non-monotonic behavior for $\pi$ plasmon lifetimes in pure graphene is not yet well understood.

Finally, we emphasize again that the lifetime studies reported in this work is purely electronic at the RPA level, and hence the conclusion might not be directly applicable to realistic situations. Yet this study is meaningful from a numerical point of view, since a reliable determination of plasmons line-widths from first-principles calculations appears to be a challenging task. Our procedure could provide reference numbers for simple situations. Plasmon lifetime studies incorporating effects from phonons, impurities, and disorders Jablan et al. (2009); Principi et al. (2014), as well as electronic effects beyond RPA have been reported in literature at model levelPrincipi et al. (2013). Including these effects in a first-principles way is in principle also feasible, but goes beyond the scope of this work.

V Summary

Plasmon excitations in graphene and graphene-related materials have intriguing properties that hold great promises for technological applications. In this work, we developed a first-principles TDDFT-RPA module within the ABACUS software package aba , which allows for accurately simulating the plasmon excitations in graphene and graphene/hBN heterostructure for a large energy window (from 0 to $\sim 30$ eV). Regarding the controversial nature of the dispersion behavior of $\pi$ plasmons, our first-principles results and model analysis indicated that a theoretically more sound dispersion relation should be $\omega_{\pi}(q)=\sqrt{E_{g,M}^{2}+\beta q}$ at small $q$’s, in stark distinction from previous proposals. The essential physics behind this is that the $\pi$ plasmons in graphene arise from collective interband excitations that has a finite gap, and hence the asymptotic behavior at ${\bf q}\rightarrow 0$ is different from that of Dirac plasmons, which come from collective gapless intraband transitions. Finally, we demonstrated that, to extract accurate lifetime from the computed spectra, care must be taken to extrapolate the results to the limit of $\eta\rightarrow 0$, where $\eta$ is the technical broadening parameter used in dynamical response function calculations.

Appendix A Comparison of graphene/h-BN heterostructure with different stackings.

In the main text, for simplicity, we only presented the computational results of the graphene/h-BN heterostructure with AA stacking. For completeness, here we also present in Fig. 7 the results of dispersion relations and lifetimes (FWHM’s of the peaks) for the Dirac plasmons with other types of stackings. Figure. 7 shows that the different stackings have little impact on both the dispersion behavior and lifetime. We also did a similar comparison study for $\pi$ plasmons, and arrived at essentially the same conclusion.

Appendix B On the dispersion of the $\pi$ plasmons in graphene

We start with a tight-binding model description of graphene

[TABLE]

where $t$ is the hopping integral between neighboring sites on a honeycomb lattice. The resultant dispersion of the $\pi$ (valence) and $\pi^{\ast}$ (conduction) bands reads

[TABLE]

with $a$ being the distance between neighboring C atoms, and ${\bf k}=(k_{x},k_{y})$ are the two-dimensional wavevectors in the BZ. The real part of the dielectric function of graphene given by such a tight-binding model, within the random-phase approximation, can be obtained as

[TABLE]

where $\epsilon_{v,{\bf k}}$/$\epsilon_{c,{\bf k}}$ are the valence/conduction band energies as given by Eq. (15), and $|v,{\bf k}\rangle$/$|c,{\bf k}\rangle$ are the corresponding eigenvectors. We note that that $(k_{x}=0,k_{y}=2\pi/3a)$ corresponds to one $M$ point in the BZ. The excitation energy gap in the vicinity of this $M$ point can be obtained by a Taylor expansion,

[TABLE]

where ${\bf k}^{\prime}=(k^{\prime}_{x},k^{\prime}_{y})=(k_{x},k_{y}-2\pi/3a)$ is the coordinate of a ${\bf k}$ point with reference to the $M$ point, and $E_{g,M}$ is the direct band gap at the $M$ point. From Eq. (17) one may realize that the $M$ point is a saddle point in the band dispersion of graphene, and $E_{g,M}$ represents a van Hove singularity in the single-particle density of states (DOS) of graphene.

For $q\rightarrow 0$, the oscillator strength in the numerator of Eq. (16) becomes

[TABLE]

where $p({\bf k})=\langle v,{\bf k}|\nabla_{\bf r}|c,{\bf k}\rangle$, and $m_{e}$ is the mass of the electron. Furthermore, we are interested in the $\pi$ plasmons which originate from the collective $\pi\rightarrow\pi^{\ast}$ interband excitations around the $M$ point. Therefore the $k$-integration in Eq (16) can be restricted to the vicinity of $M$ (i.e. by setting $|{\bf k}^{\prime}|<=k_{cut}$), and consequently the dielectric function at $q\rightarrow 0$ is simplified to

[TABLE]

The precise form of the parameter $\beta$ introduced in Eq. (19) in principle depends on ${\bf q}$ and $\omega$, but the dependence is of higher order. To the zeroth order approximation, $\beta$ can be taken to be a constant, which is valid for small $q$’s. Also, since a proper value for the cutting wavevector $k_{cut}$ is not known a priori, we cannot determine $\beta$ reliably within this model. Therefore, in this work $\beta$ is treated as a fitting parameter.

The plasmon dispersion relation can be determined by searching for the zeros of the dielectric function ${Re}\{\varepsilon({\bf q},\omega)\}=0$. Hence we end up with,

[TABLE]

or

[TABLE]

Appendix C Procedure to determine FWHM

In this section, we describe in further details how the plasmons lifetimes (or more precisely the FWHM’s) are determined in our work. First, the loss spectra were calculated for several different $\eta$ parameters. For a given $\eta$, the loss spectra calculations are converged with other technical parameters such as the number of ${\bf k}$-point sampling in BZ, and the number of frequency grid points. The position and height ($h$) of the spectral peak center are then determined, together with the positions of half height ($h/2$) on the left- and right-hand sides of the peak center, denoted as $\omega_{l}$ and $\omega_{r}$ respectively. The FWHM for a given $\eta$ and ${\bf q}$ is then given by $\omega_{r}-\omega_{l}$. The final FWHM value is then obtained by extrapolating the FWHM values for finite $\eta$’s to $\eta=0$.

Figure 8 shows the FWHM values of the conventional Dirac plasmon peaks in graphene for different $\eta$ values. One can clearly see that the FWHM value depends appreciably on the parameter $\eta$. For $\eta=0.01$ Rydberg, a typical default choice for first-principle TDDFT-RPA calculations, the FWHM value is overestimated by more than 0.25 eV for small $q$’s. As $\eta$ is reduced, the FWHM value decreases steadily to zero. From Fig. 8 one can also see their is a critical momentum value $q_{c}$, below which the FWHM value approaches all the way to zero as $\eta$ goes to zero.

Figure 9 illustrates how the FWHM value is extrapolated to $\eta=0$ limit. Below the critical $q$ value $q_{c}=0.107\AA^{-1}$, the actual FWHM value is linear proportional to $\eta$ (actually roughly $2\eta$). Thus the FWHM naturally goes to zero at the $\eta=0$ limit. For $q>q_{c}$, the $\eta$-dependence of FWHM follows nicely a quadratic behavior, saturating at a finite value at $\eta=0$. The thus obtained $\eta=0$ limit of FWHM for Dirac plasmons were presented in Fig. 6. Also presented are the FWHM results for $\pi$ plasmons, obtained using a similar procedure. However, for $\pi$ plasmons the FWHM is finite for all $q$ values.

Acknowledgments

The work was supported by the National Key Research and Development Program of China (Grants No. 2016YFB0201202). XR acknowledges the support from Chinese National Science Foundation (Grant number 11374276, 11574283). LH acknowledges the support from Chinese National Science Foundation (Grant number 11374275). The numerical calculations have been done on the USTC HPC facilities.