*Article* **Resolving the Mechanism of Acoustic Plasmon Instability in Graphene Doped by Alkali Metals**

**Leonardo Maruši´c 1, \*, Ana Kalini´c 2 , Ivan Radovi´c 2 , Josip Jakovac 3 , Zoran L. Miškovi´c <sup>4</sup> and Vito Despoja 3,5**


**Abstract:** Graphene doped by alkali atoms (AC*x*) supports two heavily populated bands (*π* and *σ*) crossing the Fermi level, which enables the formation of two intense two-dimensional plasmons: the Dirac plasmon (DP) and the acoustic plasmon (AP). Although the mechanism of the formation of these plasmons in electrostatically biased graphene or at noble metal surfaces is well known, the mechanism of their formation in alkali-doped graphenes is still not completely understood. We shall demonstrate that two isoelectronic systems, KC<sup>8</sup> and CsC8, support substantially different plasmonic spectra: the KC<sup>8</sup> supports a sharp DP and a well-defined AP, while the CsC<sup>8</sup> supports a broad DP and does not support an AP at all. We shall demonstrate that the AP in an AC*<sup>x</sup>* is not, as previously believed, just a consequence of the interplay of the *π* and *σ* intraband transitions, but a very subtle interplay between these transitions and the background screening, caused by the out-of-plane interband *C*(*π*) → *A*(*σ*) transitions.

**Keywords:** plasmon; acoustic plasmon; graphene; graphene intercalation compounds; EELS

#### **1. Introduction**

Studying acoustic plasmons (APs) in single-layer [1–5], double-layer [3,6–9] and multilayer [10–13] graphene or in metal/dielectric/graphene superstructures [14] is a very active field of research. Recently, an AP was detected in a graphene–dielectric–metal structure using near-field scattering microscopy [2]. That experiment demonstrated a very robust character of the AP with small propagation loss, which enables its efficient coupling to the electromagnetic field at infra-red frequencies. The graphene/graphene-nanoribbon superstructure also supports a strong AP which can be excited by light and provide strong field enhancement inside the nano-gap, allowing efficient biosensing [6]. In this work, we focus on the AP in graphene doped by alkali metals.

As alkali atoms are added to graphene deposited on various substrates, they intercalate between the graphene layer and the substrate and metalize, i.e., they arrange in a twodimensional (2D) crystal lattice, which supports a partially filled parabolic band [15–21]. The alkali atoms then usually donate electrons to the graphene *π* band, lifting the Fermi level above the Dirac point and transforming the graphene from a gapless semiconductor into a metal. Therefore, the alkali adlayer results in the formation of two thin, i.e., quasi two-dimensional (q2D), plasmas able to support a variety of different plasmonic modes, which do not exist in pristine or electrostatically doped graphene. Our previous theoretical research of electronic excitations in graphene doped by alkali atoms (AC*x*) showed that this system supports two kinds of intraband or 2D plasmons: the Dirac plasmon (DP) and the

**Citation:** Maruši´c, L.; Kalini´c, A.; Radovi´c, I.; Jakovac, J.; Miškovi´c, Z.L.; Despoja, V. Resolving the Mechanism of Acoustic Plasmon Instability in Graphene Doped by Alkali Metals. *Int. J. Mol. Sci.* **2022**, *23*, 4770. https://doi.org/10.3390/ ijms23094770

Academic Editor: Ana María Díez-Pascual

Received: 15 March 2022 Accepted: 23 April 2022 Published: 26 April 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

acoustic plasmon (AP). The intensity of these modes depends on the type of the dielectric substrate, as well as on electrostatic or natural doping [22–25]. Moreover, these theoretical studies suggest the possibility of simple manipulation of the DP and AP intensities (they can be switched 'on' and 'off' in a controlled way) [25], thereby opening new possibilities for their application in various fields, such as plasmonics, photonics, transformation optics, optoelectronics, light emitters, detectors and photovoltaic devices [26–36]. Additionally, these 'tunable' 2D plasmons could be very useful in the area of chemical or biological sensing [27,37–40].

In this paper, we focus on the AP, firstly because of its linear dispersion, but also because of the still rather unclear mechanism of its occurrence in the AC*x*. The mechanism of formation of the AP on a metal surface was resolved a long time ago [41–43]. In that case, a partially filled surface state supports a 2D plasmon, which hybridizes with the ordinary surface plasmon, resulting in the transformation of the usual 2D plasmon with square-root 2D dispersion into a linear AP [41]. However, our previous studies of plasmons in AC*<sup>x</sup>* systems show that the AP is very sensitive and strongly dependent on several parameters: the type of the alkali metal atom (A), the coverage (*x*) and the electrostatic doping. In other words, it is substantially unstable, and the aim of this research is to explore the mechanism of the (in)stability of the AP in the AC*x*.

In order to do that, we investigate two model systems of doped graphene, KC<sup>8</sup> and CsC8, which are isostructural, but the first system does support the AP, while the second one does not. Moreover, the DP in the KC<sup>8</sup> is very sharp, while in the CsC<sup>8</sup> it is quite broad, even though the occupancies of the graphene *π* band and the alkali atom *σ* band in these two systems are almost identical. We shall show that the AP in the AC*<sup>x</sup>* is not, as previously believed, just a consequence of the hybridization between the intraband 2D plasmons in the *π* and *σ* bands, but rather a very subtle interplay between these plasmons and the background screening caused by the interlayer interband *C*(*π*) → *A*(*σ*) transitions. We shall also demonstrate that the electronic screening coming from the high-energy interband transitions significantly reduces the intensity of the AP in both systems.

We study electronic excitations in the AC*<sup>x</sup>* using two complementary methods: the ab initio random phase approximation (RPA) approach [24] and a reduced model. In the reduced model, the electrons in the alkali atom layer are approximated by a parabolic band with parameters (the effective mass *m<sup>σ</sup>* and the Fermi energy *EFσ*) taken from the ab initio calculations. We refer to this approximation as the "massive" 2D electron gas (m2DEG). The graphene band structure is approximated by the conical *π* bands with the occupation corresponding to the occupation of the *π* band (Fermi energy *EFπ*) in the ab initio AC<sup>8</sup> sample. This approximation is also known as the 'massless' Dirac fermion (MDF) approximation [44,45]. The polarization effects coming from the highenergy interband transitions are included through the polarizability parameters *α<sup>m</sup>* and *α<sup>g</sup>* deduced from the ab initio calculations.

The paper is organized as follows. In Section 2, we present the geometry of the system, a derivation of the ab initio ground state and the RPA spectra of the electronic excitations *S*(*Q*, *ω*) in the AC*<sup>x</sup>* deposited on a dielectric substrate. We also present a theoretical formulation of the reduced model. In Section 3, we show the band structures of the KC<sup>8</sup> and CsC8, as well as the intensities of the electronic excitations in these systems. After that, using a detailed analysis, which combines the ab initio RPA method and the reduced model, we resolve the mechanism of the AP instability. The conclusions are presented in Section 4.

#### **2. Theoretical Formulation**

In this section, we derive the spectral function of the electronic excitations in the self-standing KC<sup>8</sup> and CsC8. However, to demonstrate the robustness of these electronic excitations, we show how a dielectric substrate can be included in our expressions for the spectral function, to enable the calculation of the electronic excitation spectra for systems consisting of crystals (CsC<sup>8</sup> and KC8) deposited on Al2O<sup>3</sup> substrate described by a local dielectric function *ǫS*(*ω*). Additionally, we make use of the plausible assumption [25] that

the Al2O<sup>3</sup> substrate does not affect the ground state electronic structure (Kohn–Sham (KS) wave functions and energies) of the AC<sup>8</sup> crystal, but it can influence its dielectric properties.

Figure 1 shows the geometry of the studied system. The coordinates are oriented so that the AC<sup>8</sup> crystal is positioned in the *x* − *y* plane, the *z* direction is perpendicular to the crystal plane, the graphene layer occupies the *z* = 0 plane, the alkali atom layer occupies the *<sup>z</sup>* <sup>=</sup> <sup>−</sup>*<sup>d</sup>* plane and the dielectric substrate occupies the *<sup>z</sup>* <sup>&</sup>lt; <sup>−</sup>*<sup>h</sup>* half-space.

**Figure 1.** Schematic representation of the AC<sup>8</sup> crystal deposited on a dielectric substrate described by a local dielectric function *<sup>ǫ</sup>S*(*ω*). The substrate (blue) occupies the region *<sup>z</sup>* <sup>&</sup>lt; <sup>−</sup>*h*, the graphene layer (yellow) is in the *z* = 0 plane and the alkali atom (turquoise) layer is in the *z* = −*d* plane.

#### *2.1. Calculation of the Surface Electronic Excitations Spectra*

We shall briefly describe the method of calculation of the surface electronic excitation spectra at the arbitrary position *z* > 0, previously used in several studies of electronic excitations in 2D crystals [22,23,46–50].

We start from the 3D Fourier transform of the noninteracting (free) electron response function

$$\begin{split} \chi^{0}\_{\mathbf{G}\mathbf{G}'} (\mathbf{Q}, \omega) = \frac{2}{\Omega} \sum\_{\mathbf{K} \in \text{SBZ}} \sum\_{n,m} \frac{f\_n(\mathbf{K}) - f\_m(\mathbf{K} + \mathbf{Q})}{\omega + i\eta + E\_n(\mathbf{K}) - E\_m(\mathbf{K} + \mathbf{Q})} \\\\ \times \rho\_{n\mathbf{K}, m\mathbf{K}+\mathbf{Q}}(\mathbf{G}) \,\rho^\*\_{n\mathbf{K}, m\mathbf{K}+\mathbf{Q}}(\mathbf{G}'), \end{split} \tag{1}$$

where <sup>Ω</sup> = *<sup>S</sup>* × *<sup>L</sup>* is the normalization volume and *<sup>f</sup>n***<sup>K</sup>** = [*<sup>e</sup>* (*En***K**−*EF*)/*kT* + 1] −1 is the Fermi–Dirac distribution at temperature *T*. The matrix elements (or the charge vertices) are

$$\rho\_{n\mathbf{K},m\mathbf{K}+\mathbf{Q}}(\mathbf{G}) = \left\langle \Phi\_{n\mathbf{K}} \Big| e^{-i(\mathbf{Q}+\mathbf{G})\cdot\mathbf{r}} \Big| \Phi\_{n\mathbf{K}+\mathbf{Q}} \right\rangle\_{V'} \tag{2}$$

where **K** = (*Kx*, *Ky*) is the 2D wave vector, **Q** is the momentum transfer vector parallel to the *<sup>x</sup>* − *<sup>y</sup>* plane, **<sup>G</sup>** = (**G**<sup>k</sup> , *Gz*) are the 3*D* reciprocal lattice vectors and **r** = (*ρ*, *z*) is the 3*D* position vector. The integration is performed over the normalization volume Ω. The plane wave expansion of the wave function has the form

$$\Phi\_{n\mathbf{K}}(\rho, z) = \frac{1}{\sqrt{\Omega}} e^{i\mathbf{K} \cdot \mathbf{\mathcal{P}}} \sum\_{\mathbf{G}} \mathbb{C}\_{n\mathbf{K}}(\mathbf{G}) e^{i\mathbf{G} \cdot \mathbf{r}},$$

where the coefficients *Cn***<sup>K</sup>** are obtained by solving the KS equations within the local density approximation (LDA) self-consistently.

The Fourier expansion of the free electron response function in the *z* and *z* ′ coordinates is:

$$
\chi^{0}\_{\mathbf{G}\_{\parallel}\mathbf{G}\_{\parallel}'} (\mathbf{Q}\_{\prime}\omega, z\_{\prime}z^{\prime}) = \frac{1}{L} \sum\_{\mathbf{G}\_{z}\mathbf{G}\_{z}'} \chi^{0}\_{\mathbf{G}\mathbf{G}^{\prime}} (\mathbf{Q}\_{\prime}\omega) e^{i\mathbf{G}\_{z}z - i\mathbf{G}\_{z}^{\prime}z^{\prime}}\,\tag{3}
$$

where we assume that our system is periodical in the *z* and *z* ′ direction as well, i.e., that it repeats periodically from supercell to supercell, and the supercells are q2D crystals separated by the distance *L*. We now need to determine the screened response function *χ* 0 **<sup>G</sup>**k**G**′ k (**Q**, *ω*, *z*, *z* ′ ) of one supercell without including the polarization of the surrounding supercells. This spurious interaction with the replicas of the q2D crystal can be eliminated easily, by using the RPA Dyson equation

$$
\chi\_{\mathbf{G}\_{\parallel}\mathbf{G}\_{\parallel}'}(\mathbf{Q},\omega,z,z') = \chi\_{\mathbf{G}\_{\parallel}\mathbf{G}\_{\parallel}'}^{0}(\mathbf{Q},\omega,z,z') + \sum\_{\mathbf{G}\_{\parallel}\mathbf{G}\_{\parallel}} \int\_{-L/2}^{L/2} dz\_{1} dz\_{2} \tag{4}
$$

$$
\times \chi\_{\mathbf{G}\_{\parallel}\mathbf{G}\_{\parallel}}^{0}(\mathbf{Q},\omega,z,z\_{1}) \, v\_{\mathbf{G}\_{\parallel}\mathbf{G}\_{\parallel}}^{\mathbf{D}}(\mathbf{Q},z\_{1},z\_{2}) \chi\_{\mathbf{G}\_{\parallel}\mathbf{G}\_{\parallel}'}(\mathbf{Q},\omega,z\_{2},z'),
$$

where the matrix of the bare Coulomb interaction is

$$v^{\text{2D}}\_{\mathbf{G}\_{\parallel 1}\mathbf{G}\_{\parallel 2}}(\mathbf{Q}, z, z') = v^{\text{2D}}\_{\mathbf{G}\_{\parallel 1}}(\mathbf{Q}, z, z') \delta\_{\mathbf{G}\_{\parallel 1}\mathbf{G}\_{\parallel 2}}\tag{5}$$

and the 2D Fourier transform of the bare Coulomb interaction is

$$v\_{\mathbf{G}\_{\parallel 1}}^{\text{2D}}(\mathbf{Q}, z, z') = \frac{2\pi}{|\mathbf{Q} + \mathbf{G}\_{\parallel 1}|} e^{-|\mathbf{Q} + \mathbf{G}\_{\parallel 1}||z - z'|}. \tag{6}$$

Since the integrations in (4) are performed from −*L*/2 to *L*/2, the interaction between the density fluctuations, via the Coulomb propagator *v* 2D **<sup>G</sup>**k1**G**k<sup>2</sup> (**Q**, *z*, *z* ′ ), is limited to one supercell located at <sup>−</sup>*L*/2 <sup>&</sup>lt; *<sup>z</sup>* <sup>&</sup>lt; *<sup>L</sup>*/2. After inserting the Fourier expansion (3), and a similar one for *χ*, in RPA Dyson Equation (4), it again becomes a matrix equation

$$\chi\_{\mathbf{G}\mathbf{G}'}(\mathbf{Q},\omega) = \chi\_{\mathbf{G}\mathbf{G}'}^{0}(\mathbf{Q},\omega) + \sum\_{\mathbf{G}\_1\mathbf{G}\_2} \chi\_{\mathbf{G}\mathbf{G}\_1}^{0}(\mathbf{Q},\omega) \, v\_{\mathbf{G}\_1\mathbf{G}\_2}^{2D}(\mathbf{Q}) \chi\_{\mathbf{G}\_2\mathbf{G}'}(\mathbf{Q},\omega),\tag{7}$$

where the matrix of the bare Coulomb interaction is

$$\begin{split} \upsilon\_{\mathbf{G}\_{1}\mathbf{G}\_{2}}^{\mathcal{D}D}(\mathbf{Q}) &= \upsilon\_{\mathbf{G}\_{1}\mathbf{G}\_{2}}^{\mathcal{D}D}(\mathbf{Q}) - p\_{\mathbf{G}\_{21}}p\_{\mathbf{G}\_{22}} \frac{4\pi(1-\varepsilon^{-}|\mathbf{Q}+\mathbf{G}\_{\parallel 1}|^{L})}{\left|\mathbf{Q}+\mathbf{G}\_{\parallel 1}\right|^{L}} \\ &\times \frac{\left|\mathbf{Q}+\mathbf{G}\_{\parallel 1}\right|^{2}-\mathbf{G}\_{\perp 1}\mathbf{G}\_{\mathbf{r} 2}}{\left(\left|\mathbf{Q}+\mathbf{G}\_{\parallel 1}\right|^{2}+\mathbf{G}\_{\perp 1}^{2}\right)\left(\left|\mathbf{Q}+\mathbf{G}\_{\parallel 1}\right|^{2}+\mathbf{G}\_{\perp 2}^{2}\right)} \delta\_{\mathbf{G}\_{\parallel 1}\mathbf{G}\_{\parallel 2}\prime} \end{split} \tag{8}$$

with *v* 3*D* **G**1**G**2 (**Q**) = <sup>4</sup>*<sup>π</sup>* |**Q**+**G**1| <sup>2</sup> *δ***G**1**G**<sup>2</sup> , *G<sup>z</sup>* = <sup>2</sup>*π<sup>k</sup> L* , *pG<sup>z</sup>* = (−1) *k* , and *<sup>k</sup>* <sup>∈</sup> <sup>Z</sup>. The solution of Equation (7) has the form

$$\chi\_{\mathbf{G}\mathbf{G}'}(\mathbf{Q},\omega) = \sum\_{\mathbf{G}\_1} \mathcal{E}\_{\mathbf{G}\mathbf{G}\_1}^{-1}(\mathbf{Q},\omega) \chi\_{\mathbf{G}\_1\mathbf{G}'}^{0}(\mathbf{Q},\omega),\tag{9}$$

where the dielectric matrix is

$$\mathcal{E}\_{\mathbf{G}\mathbf{G}'} (\mathbf{Q}\_{\prime}\omega) = \delta\_{\mathbf{G}\mathbf{G}'} - \sum\_{\mathbf{G}\_1} V\_{\mathbf{G}\mathbf{G}\_1}^{2D} (\mathbf{Q}) \chi\_{\mathbf{G}\_1\mathbf{G}'}^{0} (\mathbf{Q}\_{\prime}\omega). \tag{10}$$

After solving Equation (7), the nonlocal screened response function in the *z* and *z* ′ direction becomes:

$$
\chi\_{\mathbf{G}\_{\parallel}\mathbf{G}\_{\parallel}'} (\mathbf{Q}\_{\prime}\omega, z\_{\prime}z^{\prime}) = \frac{1}{L} \sum\_{\mathbf{G}\_{z}\mathbf{G}\_{z}'} \chi\_{\mathbf{G}\mathbf{G}^{\prime}}(\mathbf{Q}\_{\prime}\omega) e^{i\mathbf{G}\_{z}z - i\mathbf{G}\_{z}'z^{\prime}}.\tag{11}$$

The propagator of the induced dynamically screened Coulomb interaction can be calculated from the response function (9) as [50]

$$\mathcal{W}\_{\mathbf{G}\_{\parallel}}^{\rm ind}(\mathbf{Q},\omega,z,z') = \int\_{-L/2}^{L/2} dz\_1 dz\_2 v\_{\mathbf{G}\_{\parallel}}^{\rm D}(\mathbf{Q},z,z\_1) \chi\_{\mathbf{G}\_{\parallel}0}(\mathbf{Q},\omega,z\_1,z\_2) \, v\_0^{\rm 2D}(\mathbf{Q},z\_2,z'),\tag{12}$$

where the index zero means that **G**′ <sup>k</sup> <sup>=</sup> 0. After using the expansion (11) and Equation (6), the integrations over *z*<sup>1</sup> and *z*<sup>2</sup> can be performed analytically, and the induced dynamically screened interaction at *z*, *z* ′ > 0 can be written as

$$\mathcal{W}\_{\mathbf{G}\_{\parallel}}^{\rm ind}(\mathbf{Q},\omega,z,z') = e^{-|\mathbf{Q}+\mathbf{G}\_{\parallel}|z-Qz'} D\_{\mathbf{G}\_{\parallel}}(\mathbf{Q},\omega) \tag{13}$$

where the propagator of the surface excitations is

$$D\_{\mathbf{G}\_{\parallel}}(\mathbf{Q},\omega) = \sum\_{\mathbf{G}\_{\text{z1}},\mathbf{G}\_{\text{z2}}} \chi\_{\mathbf{G}\_{\parallel},0,\mathbf{G}\_{\text{z1}},\mathbf{G}\_{\text{z2}}}(\mathbf{Q},\omega) F\_{\mathbf{G}\_{\text{z1}}}(\mathbf{Q}+\mathbf{G}\_{\parallel}) F\_{\mathbf{G}\_{\text{z2}}}^{\*}(\mathbf{Q}),\tag{14}$$

and the form factors *F* are

$$F\_{G\_z}(\mathbf{Q}) = \frac{4\pi p\_{G\_z}}{Q\sqrt{L}} \frac{\sinh\left(\frac{QL}{2}\right)}{Q + iG\_z}.\tag{15}$$

The spectral function, which defines the intensity of the energy loss by an external perturbation to the excitation of the (**Q**, *ω*) modes, can now be calculated as

$$S(\mathbf{Q}\_{\prime}\omega) = -ImD\_{\mathbf{G}\_{\parallel} = 0}(\mathbf{Q}\_{\prime}\omega). \tag{16}$$

Up to this point, we derived the expressions for the self-standing q2D systems, but including a dielectric substrate polarization is now straightforward. It is obtained simply by replacing the bare Coulomb interaction (6) with the Coulomb interaction screened by the substrate

$$\imath \upsilon\_{\mathbf{G}\_{\parallel}}^{\text{2D}}(\mathbf{Q}, z, z') \rightarrow \imath\_{\mathbf{G}\_{\parallel}}^{\text{2D}}(\mathbf{Q}, \omega, z, z') = \imath\_{\mathbf{G}\_{\parallel}}^{\text{2D}}(\mathbf{Q}, z, z') + \frac{2\pi}{|\mathbf{Q} + \mathbf{G}\_{\parallel}|} D\_{\mathbf{S}}(\omega) \epsilon^{-|\mathbf{Q} + \mathbf{G}\_{\parallel}|(2k + z + z')},\tag{17}$$

where the dielectric surface response function is

$$D\_S(\omega) = \frac{1 - \epsilon\_S(\omega)}{1 + \epsilon\_S(\omega)}$$

This also means that the matrix (8) used in matrix Dyson Equation (7) should be modified as

$$\begin{split} \boldsymbol{v}^{2\text{D}}\_{\mathbf{G}\_{1}\mathbf{G}\_{2}}(\mathbf{Q}) &\to \boldsymbol{\vartheta}^{2\text{D}}\_{\mathbf{G}\_{1}\mathbf{G}\_{2}}(\mathbf{Q},\omega) = \\\\ \boldsymbol{v}^{2\text{D}}\_{\mathbf{G}\_{1}\mathbf{G}\_{2}}(\mathbf{Q}) + \frac{|\mathbf{Q} + \mathbf{G}\_{\|1\|}|}{2\pi} \boldsymbol{e}^{-2|\mathbf{Q} + \mathbf{G}\_{\|1\|}|h} D\_{\mathcal{S}}(\omega) \, F\_{\mathcal{G}\_{21}}(|\mathbf{Q} + \mathbf{G}\_{\|1\|}|) F\_{\mathcal{G}\_{22}}^{\*}(|\mathbf{Q} + \mathbf{G}\_{\|1\|}|) \delta\_{\mathbf{G}\_{\|1}\mathbf{G}\_{\|2}}. \end{split} \tag{18}$$

.

#### *2.2. Calculation of the 2D Dynamical Polarizability Function α*(*ω*)

In the long-wavelength or optical limit (*Q* → 0), the in-plane dielectric function of a 2D crystal can be approximated as

$$
\epsilon(\mathbb{Q}, \omega) = 1 + 2\pi \text{Qa}(\omega). \tag{19}
$$

The 2D polarizability *α* can be divided into intraband and interband contributions

$$\mathfrak{a}(\omega) = \mathfrak{a}^{\text{intra}}(\omega) + \mathfrak{a}^{\text{inter}}(\omega), \tag{20}$$

where

*v* 2*D*

$$\alpha^{\text{intra}/\text{inter}}(\omega) = i\sigma\_{\mu}^{\text{intra}/\text{inter}}(\omega) / \omega \,, \tag{21}$$

and *µ* = *x* or *y*. The intraband optical conductivity is [51]

$$\sigma\_{\mu}^{\text{intra}}(\omega) = i \frac{e^2}{m} \frac{n\_{\mu}}{\omega + i \eta\_{\text{intra}}},\tag{22}$$

where the effective number of charge carriers is

$$m\_{\mu} = -\frac{m}{Se^2} \sum\_{n} \sum\_{\mathbf{K} \in \text{1.SBZ}} \frac{\partial f\_{n\mathbf{K}}}{\partial E\_{n\mathbf{K}}} \left| \dot{j}\_{n\mathbf{K}, n\mathbf{K}}^{\mu} \right|^2. \tag{23}$$

The interband optical conductivity is determined from the optical limit of the nonlocal interband conductivity

$$
\sigma\_{\mu}^{inter}(\omega) = L \sigma\_{\mu}^{inter}(\omega, \mathbf{Q} \to 0).
$$

The nonlocal interband conductivity is [51]

$$\sigma\_{\mu}^{\text{inter}}(\mathbf{Q},\omega) = -\frac{i\hbar}{\Omega} \sum\_{n \neq m} \sum\_{\mathbf{K} \in \mathbb{L}.\text{SEZ}} \frac{\left| \boldsymbol{\upbeta}\_{n\mathbf{K},m\mathbf{K}+\mathbf{Q}}^{\mu} \right|^{2}}{E\_{n\mathbf{K}} - E\_{m\mathbf{K}+\mathbf{Q}}} \frac{f\_{n\mathbf{K}} - f\_{m\mathbf{K}+\mathbf{Q}}}{\hbar\omega + i\eta\_{inter} + E\_{n\mathbf{K}} - E\_{m\mathbf{K}+\mathbf{Q}}},\tag{24}$$

where the current vertices are

$$j\_{n\mathbf{K},m\mathbf{K}+\mathbf{Q}}^{\mu} = \int\_{\Omega} d\mathbf{r} e^{-i\mathbf{Q}\cdot\mathbf{r}} j\_{n\mathbf{K},m\mathbf{K}+\mathbf{Q}}^{\mu}(\mathbf{r}) \tag{25}$$

and the current produced by the transitions between the Bloch states *φ* ∗ *<sup>n</sup>***<sup>K</sup>** → *φm***K**+**<sup>Q</sup>** is defined as

$$j\_{n\mathbf{K},m\mathbf{K}+\mathbf{Q}}^{\mu}(\mathbf{r}) = \frac{e\hbar}{2im} \{\phi\_{n\mathbf{K}}^{\*}(\mathbf{r})\partial\_{\mu}\phi\_{m\mathbf{K}+\mathbf{Q}}(\mathbf{r}) - [\partial\_{\mu}\phi\_{n\mathbf{K}}^{\*}(\mathbf{r})]\phi\_{m\mathbf{K}+\mathbf{Q}}(\mathbf{r})\}\dots$$

Figure 2 shows the interband contribution to the dynamical polarizability *α* inter(*ω*) in the KC<sup>8</sup> (black), CsC<sup>8</sup> (orange) and doped graphene (red dashed). The graphene is doped so that the Fermi level is 1eV above the Dirac point, which corresponds to the doping of the *π* bands in the KC<sup>8</sup> and CsC8. All three systems show qualitatively equal behavior; the peak at about *ω* = 2 eV, indicating the onset for the *π* → *π* ∗ interband transitions, and the dip at *ω* ≈ 4 eV is a consequence of the high density of the *π* → *π* ∗ interband transitions at the M point of the Brillouin zone. Even though the *π* bands in all three systems are almost equally doped, we can see a substantial difference in the KC<sup>8</sup> and CsC<sup>8</sup> statical polarizabilities which are *α* inter(*ω* = 0) = 2.2 Å and *α* inter(*ω* = 0) = 3.05 Å, respectively. This difference probably comes from the difference in the intensities of the *C*(*π*) → *A*(*σ*) interlayer (interband) excitations in the two systems. These excitations are manifested as the peak at *ω* ≈ 0.25 eV for the CsC8, which does not exist for the KC8. Finally, we can see that the agreement between the dynamical polarizabilities in the KC<sup>8</sup> and the equivalently doped graphene is almost perfect. As we shall demonstrate later, this small deviation in the low-energy part of the dynamical polarizability is responsible for the disappearance of the AP in the CsC8, but only if we take into account that these transitions represent a perpendicular polarization.

**Figure 2.** The interband contribution to the dynamical polarizability *α* inter(*ω*) for the KC<sup>8</sup> (black), CsC<sup>8</sup> (orange) and doped graphene (red dashed). The graphene is doped so that the Fermi energy is 1eV above Dirac point, which corresponds to the doping of the *π* bands in the KC<sup>8</sup> and CsC8.

#### *2.3. Calculation of the Substrate Dielectric Function*

We assume that the dielectric media is vacuum (i.e., *ǫ*<sup>0</sup> = 1), and that the substrate is aluminium oxide Al2O<sup>3</sup> described by the macroscopic dielectric function *ǫs*(*ω*). To calculate the *ǫs*(*ω*), we start from the 3D Fourier transform of the independent electron response function

$$\chi^{0}\_{\mathbf{GC}'}(\mathbf{q},\omega) = \frac{2}{\Omega} \sum\_{\mathbf{k} \in \mathbb{L}, \mathbf{BZ}} \sum\_{n, \mathbf{m}} \frac{f\_{\text{n}}(\mathbf{k}) - f\_{\text{m}}(\mathbf{k} + \mathbf{q})}{\omega + i\eta + E\_{\text{n}}(\mathbf{k}) - E\_{\text{m}}(\mathbf{k} + \mathbf{q})} \rho\_{n\mathbf{k},\text{m} + \mathbf{q}}(\mathbf{G}) \,\rho^{\*}\_{n\mathbf{k},\text{m} + \mathbf{q}}(\mathbf{G}'), \tag{26}$$

where **k** ∈ 1.BZ indicates that the summation is performed within the first Brillouin zone. The charge vertices are defined as

$$\rho\_{n\mathbf{k},m\mathbf{k}+\mathbf{q}}(\mathbf{G}) = \int\_{\Omega} d\mathbf{r} \, \phi\_{n\mathbf{k}}^{\*}(\mathbf{r}) e^{-i(\mathbf{q}+\mathbf{G}) \cdot \mathbf{r}} \phi\_{m\mathbf{k}+\mathbf{q}}(\mathbf{r}) \,. \tag{27}$$

where **k** = (*kx*, *ky*, *kz*), **q** = (*qx*, *qy*, *qz*) and **G** = (*Gx*, *Gy*, *Gz*) are the 3D wave vector, the transfer wave vector and the reciprocal lattice vector, respectively. The integration is performed over the normalization volume Ω. We use the response matrix (26) to determine the dielectric matrix as

$$\mathcal{E}\_{\mathbf{G}\mathbf{G}'} (\mathbf{q}, \omega) = \delta\_{\mathbf{G}\mathbf{G}'} - \sum\_{\mathbf{G}\_1} v\_{\mathbf{G}\mathbf{G}\_1}(\mathbf{q}) \chi^0\_{\mathbf{G}\_1\mathbf{G}'}(\mathbf{q}, \omega), \tag{28}$$

where the bare Coulomb interaction is *<sup>v</sup>***GG**′(**q**) = <sup>4</sup>*<sup>π</sup>* |**q+G**| <sup>2</sup> *δ***GG**′ . Finally, the macroscopic dielectric function is determined by inverting the dielectric matrix

$$
\epsilon\_s(\omega) = \epsilon\_1(\omega) + i\epsilon\_2(\omega) = 1/\mathcal{E}\_{\mathbf{G}=0\mathbf{G}'=0}^{-1}(\mathbf{q} \approx 0, \omega). \tag{29}
$$

#### *2.4. Reduced Model*

Analytical modeling of the energy loss spectra is achieved by representing each of the systems, CsC<sup>8</sup> and KC8, by a two-layer structure consisting of a single sheet of doped graphene and an m2DEG, placed in vacuum at a distance *d* apart, as shown in Figure 3. In the reduced model, the substrate is neglected, i.e., we assume *ǫS*(*ω*) = 1, and the energy

loss function −ℑ{1/*ǫ*(**Q**, *ω*)} is then obtained from the effective 2D dielectric permittivity for this two-layer structure, in the RPA given by [44,52]

$$\epsilon(\mathbf{Q},\omega) = \frac{1}{2} \left[ 1 + \coth(Qd) - 2v\_Q \chi\_{\mathcal{S}}^0 \right] - \frac{1}{2} \frac{\text{cosech}^2(Qd)}{1 + \coth(Qd) - 2v\_Q \chi\_m^0},\tag{30}$$

where *v<sup>Q</sup>* = 2*π*/*Q* is the Coulomb interaction, while *χ* 0 *g* (*Q*, *ω*) and *χ* 0 *<sup>m</sup>*(*Q*, *ω*) are the response functions of the noninteracting electrons in the graphene and the m2DEG layers, respectively.

**Figure 3.** Schematic representation of the reduced model. The alkali atom layer is approximated by 'massive' 2D electron gas (parabolic *σ* band), and the graphene layer is described by the 'massless' Dirac fermion (MDF) approximation (conical *π* band).

For the doped graphene, we follow the method proposed by Gjerding et al. [53], and write the response function as *χ* 0 *g* (*Q*, *<sup>ω</sup>*) = *<sup>χ</sup>*Dirac(*Q*, *<sup>ω</sup>*) <sup>−</sup> *<sup>α</sup>gQ*<sup>2</sup> , where *χ*Dirac is the response function given in Refs. [54,55], which describes both intraband and low-energy interband electron transitions within the *π* electron bands approximated by the Dirac cones with the Fermi energy *EFπ*, while *α<sup>g</sup>* is the phenomenological parameter providing the correction due to the high-energy interband transitions. For the m2DEG, we similarly write *χ* 0 *<sup>m</sup>*(*Q*, *<sup>ω</sup>*) = *<sup>χ</sup>*2DEG(*Q*, *<sup>ω</sup>*) <sup>−</sup> *<sup>α</sup>mQ*<sup>2</sup> , where *χ*2DEG is the polarization function given in Ref. [56], describing the intraband transitions in the 2DEG which occupies a single parabolic energy band with the effective mass *m<sup>σ</sup>* and the Fermi energy *EFσ*, while *α<sup>m</sup>* is a phenomenological parameter taking into account interband transitions in the m2DEG. The expressions used for both response functions, *χ* 0 *g* (*Q*, *ω*) and *χ* 0 *<sup>m</sup>*(*Q*, *ω*), are formulated for zero temperature, but they are corrected by the Mermin procedure to take into account a finite damping parameter *η* in both graphene and m2DEG layers [53]. We use *η* = 40 meV as in the ab initio calculations, as well as *η* = 65 meV to take into account the additional smearing at room temperature.

We note that the relevant parameters for the KC<sup>8</sup> (*d* = 2.92 Å, *EF<sup>π</sup>* = 1.01 eV, *EF<sup>σ</sup>* = 0.9 eV, *m<sup>σ</sup>* = 0.92) and CsC<sup>8</sup> (*d* = 3.13 Å, *EF<sup>π</sup>* = 1.03 eV, *EF<sup>σ</sup>* = 1.03 eV, *m<sup>σ</sup>* = 0.72) are obtained from the electronic band structure calculation for these two systems, shown in Figure 4a,b, respectively. The parameters *α* are obtained from the static limit (*ω* → 0) of the ab initio results for the corresponding dynamic polarizability functions *α* inter(*ω*) in the optical limit. In particular, the value *α<sup>g</sup>* ≈ 1.3 Å for high-energy interband transitions is deduced from the data for the intrinsic graphene [53], and by adding the value for the low-energy *π* → *π* ∗ interband transitions, estimated from the Dirac model in the optical limit at zero temperature [57] as <sup>4</sup>*πEFg*−<sup>1</sup> <sup>≈</sup> 1.15 Å, we obtain that the total contribution of the interband transitions in the doped graphene is around 2.45 Å. This is close to the

result *α* inter(*<sup>ω</sup>* <sup>=</sup> <sup>0</sup>) = 2.2 Å shown in Figure 2 for the KC8, which indicates that *<sup>α</sup><sup>m</sup>* <sup>≈</sup> <sup>0</sup> for that system. On the other hand, the result *α* inter(*ω* = 0) = 3.05 Å shown in Figure 2 for CsC<sup>8</sup> indicates that *α<sup>m</sup>* ≈ 0.6 Å for that system.

**Figure 4.** The electronic band structure in (**a**) KC<sup>8</sup> and in (**b**) CsC8. The magenta line denotes the lowest unoccupied band (LUCB). The blue circles represent the parabolic fit *E*(*K*) = *EF<sup>σ</sup>* + *<sup>h</sup>*¯ 2*K* 2 2*m<sup>σ</sup>* of the alkali atom *σ* band.

#### *2.5. Computational Details*

The KC8, CsC<sup>8</sup> and graphene KS wave functions *φn***<sup>K</sup>** and energies *En***<sup>K</sup>** are determined using the plane-wave self-consistent field DFT code (PWSCF) within the QUANTUM ESPRESSO (QE) package [58]. The core–electron interaction is approximated by the normconserving pseudopotentials [59,60]. For the KC<sup>8</sup> and CsC8, the exchange correlation (XC) potential is approximated by the Perdew–Burke–Ernzerhof (PBE) generalized gradient approximation (GGA) functional [61]. The ground state electronic densities are calculated using the 8 × 8 × 1 Monkhorst–Pack K-mesh [62] and the plane-wave cut-off energy is 60Ry. For both AC<sup>8</sup> crystals, we used the hexagonal Bravais lattice, where *a* = 4.922 Å and the

separation between the AC<sup>8</sup> layers is *L* = 2.5*a*. The atomic and the unit cell relaxations were performed until maximum force below 0.001 Ry/a.u. was obtained. After performing the structural optimization, the obtained separations between the K and Cs layers and the graphene layer are d = 2.92 Å and d = 3.13 Å respectively. The graphene XC potential is approximated by the Perdew–Zunger (PZ) LDA [63]. The ground state electronic density is determined using 21 × 21 × 1 K-point mesh, and for the plane-wave cut-off energy we choose 50 Ry. For the graphene unit-cell constant we use the experimental value *a<sup>g</sup>* = 2.45 Å [64], while for the superlattice unit-cell constant we take *L* = 5*ag*.

The AC<sup>8</sup> response functions (1) and conductivities (22)–(25) are evaluated from the wave functions <sup>Φ</sup>*n***K**(**r**) and energies *<sup>E</sup>n*(**K**) calculated for the 201 × 201 × 1 Monkhorst– Pack K-point mesh. The band summation (*n*, *m*) is performed over 60 bands, the damping parameters are *η* = *ηintra* = *ηinter* = 10 meV and the temperature is *T* = 25 meV. Due to the large spatial dispersivity of the dielectric response in the perpendicular (*z*) direction, the crystal local field effects are taken into account only in the *z* direction and neglected in the *<sup>x</sup>* − *<sup>y</sup>* plane, i.e., we set the **<sup>G</sup>**<sup>k</sup> to zero. To calculate the matrix *χGz*,*G*′ *<sup>z</sup>* we set the energy cut-off to 10Ry, which corresponds to a 23 × 23 matrix. Graphene conductivity (22)–(25) is calculated by performing the summation over 601 × 601 × 1 K-point mesh, the band summations (*n*, *m*) are performed over 20 bands, the damping parameters are *η*intra = 10 meV and *η*inter = 25 meV and the temperature is *T* = 25 meV. The conductivity of the doped graphene is calculated using the rigid band approximation, i.e., the occupation parameter (the Fermi energy relative to the Dirac point) is adjusted a posteriori.

The ground state electronic density of the bulk Al2O<sup>3</sup> is calculated using 9 × 9 × 3 K-mesh, the plane-wave cut-off energy is 50Ry and the Bravais lattice is hexagonal (12 Al and 18 O atoms in the unit cell) with the lattice constants *a* = 4.76 Å and *c* = 12.99 Å. The response function (26) of the Al2O<sup>3</sup> is calculated using the 21 × 21 × 7 *k*-point mesh and the band summations (*n*, *m*) are performed over 120 bands. The damping parameter is *η* = 100 meV and the temperature is *T* = 10 meV. For the optically small wave vectors **q** ≈ 0, the crystal local field effects are negligible, so the crystal local field effect cut-off energy is set to zero. Using this approach, the Al2O<sup>3</sup> dielectric function is estimated to be *<sup>ǫ</sup>S*(*ω*) <sup>≈</sup> *<sup>ǫ</sup>S*(*<sup>ω</sup>* <sup>=</sup> <sup>0</sup>) = 2.9 for *<sup>ω</sup>* <sup>&</sup>lt; 2 eV, which is in good agreement with the experimental value 3.4 [65]. However, in the calculations we used the full dynamical *ǫS*(*ω*).

#### **3. Results and Discussion**

Figure 4 shows (a) KC<sup>8</sup> and (b) CsC<sup>8</sup> band structures, respectively, along the high symmetry <sup>Γ</sup> → K → M → <sup>Γ</sup> directions. The blue circles represent the parabolic fit *E*(*K*) = *EF<sup>σ</sup>* + *<sup>h</sup>*¯ 2*K* 2 2*m<sup>σ</sup>* for the alkali atom *σ* bands. In both figures, we can see that the graphene *π* band and the alkali atom *σ* band are significantly doped by the electrons, so the Fermi level is around 1eV above the Dirac point and the *σ* band bottom. It is also evident that the *σ* bands behave almost as an ideal 2D free electron gas. They are parabolic (described by effective masses *m<sup>σ</sup>* = 0.916 and *m<sup>σ</sup>* = 0.72, respectively) up to the Fermi energy, especially in Figure 4a where the *σ* band is parabolic almost through the entire Brillouin zone. Finally, the most important feature is the obvious similarity between the two band structures, which both fulfill conditions for the occurrence of the AP (two bands of different effective masses crossing the Fermi energy) [22,41]. However, as we shall see, the spectra of the low-energy electronic excitations in these two systems are quite different. In both systems, we can identify the lowest band that remains above the Fermi lever for all wave-vectors, and we call it the lowest unoccupied band (LUCB), even though for some wave-vectors there is another unoccupied band below that one, since this one turns out to be of particular importance, as explained below.

Figure 5 shows the intensities of the electronic excitations in the (a) KC<sup>8</sup> and (b) CsC<sup>8</sup> deposited on the dielectric Al2O<sup>3</sup> surface. Since the graphene *π* bands are abundantly doped by the electrons, both systems support a strong DP. However, although the occupancies of both *π* bands are almost identical (*EF<sup>π</sup>* = 1.01 eV for the KC<sup>8</sup> and 1.03 eV for the CsC8), the intensities of the DP are quite different. We can see that the DP in the CsC<sup>8</sup>

is broader and much more efficiently damped by the interband *π* → *π* ∗ electron hole excitations than the DP in the KC8. Moreover, we can see that the DP in the KC<sup>8</sup> is very sharp and it extends deep into the interband *π* → *π* ∗ continuum while the DP in the CsC<sup>8</sup> decays immediately upon entering the interband *π* → *π* ∗ continuum. The most interesting difference between the two systems is that, even though the occupancies of their *σ* bands are very similar (*EF<sup>σ</sup>* = 0.9 eV in the KC<sup>8</sup> and 1.03 eV in the CsC8), the KC<sup>8</sup> does support the AP while the CsC<sup>8</sup> does not. Moreover, we can see that the *π* plasmon (denoted by the blue symbol *π*) in the KC<sup>8</sup> is well defined, and it appears to be a very intensive optically active plasmon (its intensity does not vanish for *Q* → 0), while the *π* plasmon in the CsC<sup>8</sup> is less intensive, more diffuse and not in optically active mode. It is very interesting that these two very similar band structures lead to very different excitation spectra.

To demonstrate that these effects are not driven by the dielectric substrate, in Figure 5 we also show the intensities of the electronic excitations in the self-standing (c) KC<sup>8</sup> and (d) CsC8. The magenta dots in Figure 5a,b denote the DP and AP dispersion relations in the self-standing samples, for comparison. It is worth noting that the insulator surface does not change the qualitative behavior of the electronic modes. The only quantitative differences are: the dielectric screening slightly reduces the DP energy and slightly reduces the intensities of the AP and *π* plasmons. However, the AP in the KC<sup>8</sup> is still clearly visible. Since the substrate does not affect the qualitative behavior of the excitation spectra, it will be omitted from further consideration. Therefore, the fact that one system does and the other does not support the AP, even though the bands crossing the Fermi energy in the two systems are almost identical, suggests that the mechanism responsible for this AP instability is probably in the screening coming from the interband excitations beyond the Fermi energy.

**Figure 5.** The ab initio spectra of the electronic excitations *S*(**Q**, *ω*) in the (**a**) KC<sup>8</sup> and (**b**) CsC<sup>8</sup> deposited on dielectric Al2O<sup>3</sup> surface with *h* = 5.92 Å and *h* = 6.13 Å respectively (i.e., the separation between the dielectric surface and the alkali atom layer is chosen to be 3 Å ). The ab initio spectra of the electronic excitations in the self-standing (**c**) KC<sup>8</sup> and (**d**) CsC8. The magenta dots in (**a**,**b**) denote the DP and AP dispersion relations in the self-standing samples, i.e., the DP and AP in (**c**,**d**). The energy loss function −ℑ{1/*ǫ*(**Q**, *ω*)} in the self-standing (**e**) KC<sup>8</sup> and (**f**) CsC<sup>8</sup> was obtained using the reduced model.

To investigate this mechanism, we take advantage of the reduced model which clearly distinguishes between the different interband contributions to the dynamical response in the two systems. As discussed in Section 2.4, the statical polarizabilities due to the highenergy interband excitations in the metallic subsystem (characterized by the parameter *αm*) are *α<sup>m</sup>* ≈ 0 in the KC<sup>8</sup> and 0.6 in the CsC8. This difference suggests that this polarization mechanism may be responsible for suppressing the AP. Figure 5e,f show the energy loss function −ℑ{1/*ǫ*(**Q**, *ω*)} in the self-standing KC<sup>8</sup> and CsC8, respectively, obtained using the reduced model. We can see that the agreement between the ab initio and the reduced model intensities for the KC<sup>8</sup> for *ω* < 3 eV and *Q* < 0.14 a.u. is almost perfect. Beyond these values, the plasmons in the reduced model appear at higher energies, which is reasonable since the reduced model neglects the dynamical effects of the high-energy interband transitions. For example, one can notice the absence of the *π* plasmon in the reduced model. In the case of the CsC8, the agreement is no longer so good, and the most important difference is that the reduced model (at this level of approximation) is obviously not able to reproduce the disappearance of the AP. The stronger metallic interband screening in the CsC<sup>8</sup> is still not sufficient to suppress the AP. Moreover, even the implementation of the ab initio dynamic polarizability *α* inter(*ω*) in the reduced model (taking into account that the low-energy interband excitations in the *χ*Dirac and *χ*2DEG need to be extracted) fails to cause the disappearance of the AP. This means that the mechanism of the AP disappearance is more complex and obviously the explanation requires more accurate treatment of the spatial dispersivity of the interband dynamical response, beyond the 2D model. Therefore, we shall again exploit the ab initio method, which incorporates the effects of the local crystal field (the spatial dispersivity of the dynamical response in the *z* direction), where the interband screening will be modified by changing the number of valence bands participating in the interband screening.

#### *Resolving the Mechanism of the AP Instability*

To explore how the screening coming from the interband excitations beyond the Fermi energy influences the AP, we omit one or more unoccupied bands from the calculation, to determine exactly how each band influences the excitation spectra. As we can see from the band structures shown in Figure 4a,b, in both systems, CsC<sup>8</sup> and KC8, there is one band crossing the Fermi level, i.e., for some wave-vectors that band is the highest occupied valence band (HOVB), while for other wave-vectors that same band is the lowest unoccupied band. Therefore, to avoid confusion, we use the expression lowest unoccupied band for the next one, i.e., the first band above the Fermi level at the Γ point (magenta line in Figure 4a,b), since that one remains above the Fermi level for all wave-vectors. That particular band is the one that has to be omitted from the calculation to achieve a qualitative difference with respect to the complete calculation (the one taking into account all bands). However, the difference becomes much more significant if we omit more bands. To keep track of the bands omitted from the calculation, we introduce the integer *n* which denotes the number of omitted bands. For example, *n* = 0 denotes that no bands are omitted, *n* = 1 denotes that only the lowest unoccupied band is omitted, *n* = 2 denotes that the first two unoccupied bands are omitted, etc.

Figure 6a,b show how the AP peak changes as we omit the unoccupied bands from the calculation, for *Q* = 0.1 a.u. and for the CsC<sup>8</sup> and KC8, respectively. The thick black lines show the actual spectra, i.e., the ones obtained without omitting any bands (*n* = 0), and we can see that for the selected wave-vector the AP in the KC<sup>8</sup> is clearly present (though not very strong), while in the CsC<sup>8</sup> it does not exist. On the other hand, the red line shows that the AP peak exists in both systems if we omit the LUCB from the calculation (*n* = 1). In the CsC8, that peak is very weak but it exists, while in the KC<sup>8</sup> it is, surprisingly, lower than the actual (*n* = 0) peak. The other lines show that the intensity of the peak increases as we omit more and more bands, in both systems, and it is the strongest when all the unoccupied bands are omitted (*n* = ∞). The same can be seen in Figure 6c, which shows how the intensity of the peak changes with *n*. It is important to point out that we also attempted

to keep the lowest unoccupied band in the calculation while omitting any number of the other unoccupied bands, above that one, and that did not lead to the occurrence of the AP in the CsC8. This means that omitting the lowest unoccupied band (*n* = 1) is crucial for the occurrence of the AP, while omitting the higher bands as well (*n* > 1) only enhances its intensity. This leads to the conclusion that the mechanism responsible for the disappearance of the AP in the CsC<sup>8</sup> is the small difference in the 'out-of-plane' polarization coming from the interband (interlayer) transitions between the graphene *π* band and the lowest unoccupied alkali atom *σ* band (denoted by red arrows in Figure 4). This difference is also manifested as the small peak at *ω* ≈ 0.25 eV in the CsC<sup>8</sup> dynamical polarizability *α* inter(*ω*) shown in Figure 2, which is missing for the KC8. However, even if we include this difference in the dynamical polarizability in the reduced model, that still fails to reproduce the disappearance of the AP. This is because the reduced model is inherently a 2D model, i.e., it allows only the 'in-plane' polarization, while the *π* → *σ* transitions induce the 'out-of-plane' polarization.

**Figure 6.** AP peak for *Q* = 0.1a.u. when various numbers of the unoccupied bands are omitted from the calculation (*n*) for the (**a**) CsC<sup>8</sup> and (**b**) KC8. (**c**) Intensity of the AP peak as a function of *n* for the CsC<sup>8</sup> (black line) and KC<sup>8</sup> (red line).

Another interesting point is that the difference in the intensities of the *π* → *σ* transitions is obviously correlated with the strength of the hybridization of the *π* and *σ* bands at the crossing point. We can see that the *π* and *σ* bands in the KC<sup>8</sup> intersect without any distortion while in the CsC<sup>8</sup> we notice a significant avoided crossing. In other words, the *π* and *σ* in the CsC<sup>8</sup> hybridize significantly while this is not the case in the KC8. This leads to one very unusual conclusion; the intensity of the AP, which was to date believed to be a consequence of various long-range screening mechanisms, can be significantly affected by the short-range electronic correlations occurring between the atomic orbitals (i.e., by the chemical bonding between the alkaline atoms and graphene).

#### **4. Conclusions**

We analyzed the origin of AP instability in graphene doped by alkali metals on two prototype systems, KC<sup>8</sup> and CsC8. Even though the band structures of these systems are almost identical, we proved that the hybridization between the C(*π*) and the K(*σ*) or Cs(*σ*) bands at the crossing point causes a significant modification of the dynamic polarizability in the perpendicular direction, which has a substantial influence on the lowenergy excitation spectra in these systems. We demonstrated that the net perpendicular screening in the CsC<sup>8</sup> causes the disappearance of the AP and significant weakening of the DP. We also demonstrated that the electronic screening coming from the high-energy interband transitions (beyond the HOVB–LUCB interval) significantly reduces the intensity of the AP in both systems, KC<sup>8</sup> and CsC8. This illustrates the importance of the nature of the chemical bonding between the alkaline atoms and graphene, as well as the importance of perpendicular dispersivity of the dynamical response in theoretical simulations of technologically interesting low-energy plasmons in chemically doped graphenes.

**Author Contributions:** Conceptualization, V.D., L.M., I.R. and Z.L.M.; methodology, V.D. and Z.L.M.; software, V.D., J.J., I.R. and A.K.; investigation, V.D., L.M., J.J., A.K., I.R. and Z.L.M., visualization, V.D., L.M., J.J., A.K., I.R. and Z.L.M., writing—original draft preparation, V.D., L.M. and I.R.; writing review and editing, V.D., L.M., Z.L.M. and I.R.; supervision V.D. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research is funded by the Croatian Science Foundation (Grant No. IP-2020-02-5556), the European Regional Development Fund for the 'QuantiXLie Centre of Excellence' (Grant No. KK.01.1.1.01.0004), the Ministry of Education, Science and Technological Development of the Republic of Serbia, the Serbia–Croatia bilateral project (Grant No. 337-00-205/2019-09/28), and the Natural Sciences and Engineering Research Council of Canada (Grant No. 2016-03689).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not Applicable.

**Acknowledgments:** A.K., I.R. and V.D. would like to acknowledge networking support from COST Action CA19118 EsSENce, supported by the COST Association (European Cooperation in Science and Technology). Computational resources were provided by the Donostia International Physics Center (DIPC) computing center.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

