Ab initio study of the electron energy loss function in a graphene-sapphire-graphene composite system

Vito Despoja,1,2,3,* Tijana Djordjević,4 Lazar Karbunar,5 Ivan Radović,4 and Zoran L. Mišković6 1Department of Physics, University of Zagreb, Bijenička 32, HR-10000 Zagreb, Croatia 2Donostia International Physics Center (DIPC), P. Manuel de Lardizabal, 20018 San Sebastian, Basque Country, Spain 3Universidad del Pais Vasco, Centro de Fisica de Materiales CSIC-UPV/EHUMPC, Av. Tolosa 72, E-20018 San Sebastian, Spain 4Vinča Institute of Nuclear Sciences, University of Belgrade, P.O. Box 522, 11001 Belgrade, Serbia 5School of Electrical Engineering, University of Belgrade, Bulevar Kralja Aleksandra 73, 11120 Belgrade, Serbia 6Department of Applied Mathematics, and Waterloo Institute for Nanotechnology, University of Waterloo, Waterloo, Ontario, Canada N2L 3G1 (Received 16 June 2017; published 23 August 2017)


I. INTRODUCTION
Even though graphene is just a one-atom-thick layer of carbon, it supports a variety of electronic excitations, ranging from the high-energy π and π + σ transitions, which lie in the ultraviolet (UV) to the far-ultraviolet (FUV) frequency range, down to the low-energy excitations within the π electronic bands in doped graphene, lying in the range of terahertz (THz) to infrared (IR) frequencies [1][2][3][4].It was recently shown that the π and π + σ transitions do not behave as plasmons in the long-wavelength (or optical) limit [5], whereas for smaller wavelengths they become strongly damped by electron-hole excitations [6], and thus are not interesting from a practical point of view.On the other hand, the intraband π electron excitations in doped graphene, which are often called Dirac plasmons, are weakly damped and hence exhibit relatively long propagation lengths, while being efficiently localized in two dimensions (2D) [7][8][9].Moreover, the one property that makes graphene particularly attractive is the possibility to tune its Dirac plasmon dispersion via changing the doping density of graphene by simple application of an electrostatic gate potential [3,4,[6][7][8][9][10][11][12].Another means of tuning the collective modes in devices containing multiple graphene layers is based on hybridization between Dirac plasmons in individual graphene layers, which can be controlled by suitably * vito@phy.hrengineered interlayer spacer materials of varying thickness and dielectric properties [13,14].Therefore Dirac plasmons in graphene show great potential for applications in photonics, optoelectronics, transformation optics, and plasmonics in the THz to mid-IR range [8][9][10][11].
In the nanoscale devices, graphene layers typically appear in stacks [12,15] or sandwichlike structures [16] with insulating spacers between them made of polar materials, which often support strong Fuchs-Kliewer (FK), or surface optical phonon modes in the THz to mid-IR frequency range [17].Those surface phonons strongly interact with longitudinal plasmon modes, and may therefore completely change the dispersion and damping of the Dirac plasmons in graphene [18][19][20][21], thereby seriously affecting its tunability for optoelectronic [12,17,[20][21][22] and plasmonic [23][24][25] applications in the range of frequencies of technological interest.While this provides strong motivation for a theoretical study of plasmon-phonon hybridization between graphene and the nearby polar insulator(s), a question arises as to how one can most effectively probe the hybridized plasmon-phonon modes.Using light as direct probe is not effective because the longitudinal plasmon-phonon branches will always lie below the (transverse) light line, ω = Qc, with Q being an in-plane wave number in graphene [4], so that purely optical excitation of plasmon-phonon modes will never occur.Such excitation could occur indirectly, for example, when light is used first to excite a Mie resonance in a metallic tip placed in the vicinity of a graphene/insulator interface.Such Mie resonances are localized dipole oscillations, which (depending on the polarization and frequency) are able to excite plasmon-phonon resonances at a graphene/insulator interface.This methodology is characterized as the near-field imaging [26].
However, the most efficient way to probe surface plasmonphonon resonances is by means of an externally moving charge.The most commonly used technique is the reflection electron energy loss spectroscopy (REELS), which uses electron beams with typical energies on the order of some 10 eV [27].This technique was recently used to study hybridization of a Dirac plasmon and a surface optical phonon in a graphene layer grown epitaxially on a SiC substrate, both experimentally [28][29][30] and theoretically [31,32].Another promising experimental technique, which could be used to study plasmon-phonon hybridization on a surface is lowenergy ion grazing scattering (LEIGS) where the ions with energies in the keV range are grazingly scattered with an angle of incidence of about 1 • (almost parallel) from a surface.Although this technique is mainly used for probing the FK phonons on polar insulating surfaces [33][34][35][36], it should be noted that both REELS and LEIGS use charged particles that move at speeds comparable to the Fermi speed v F ≈ c/300 (where c is the speed of light in vacuum) of graphene's π electron bands.
As a prototype structure of interest for designing photoelectronic devices that contain multiple layers of graphene, we focus in this work on a sandwich structure with infinitely large area placed in vacuum, which consists of two layers of graphene separated by a 5-nm-thick insulating layer [16], assumed to be sapphire (aluminum oxide, Al 2 O 3 ), an experimentally often used dielectric spacer [25,37,38].We designate the sandwich as sy 1 -Al 2 O 3 -sy 2 , where the system on either side of the sandwich, sy i , may be vacuum or graphene layer gr(E F ) with Fermi energy E F .When there are two graphene layers, we allow them to be doped with different charge carrier densities, giving rise to different E F values.In particular, we consider various combinations of graphene layers with E F = 1 eV, 200 meV, or 0 (neutral, or pristine graphene).We note that the gr-Al 2 O 3 -gr composite has already found an application as double-layer graphene optical modulator [37,38], where an exactly 5-nm-thick Al 2 O 3 film was fabricated [37].In this paper, we compute the so-called loss function, or the electron energy loss spectrum (EELS) [2,39], of interest for a REELS applied to one of the surfaces on our sy 1 -Al 2 O 3 -sy 2 sandwich.For example, the EELS of a composite structure gr 1 -Al 2 O 3 -gr 2 with two doped graphene layers is expected to exhibit two types of hybridization: one between the Dirac plasmons in each graphene layer, and the other being their hybridization with the surface phonon modes in the Al 2 O 3 spacer.Similar theoretical investigations were already conducted in Refs.[18,19], where the authors investigated hybridized plasmon-phonon modes in single-layer graphene deposited on an (experimentally often used) SiO 2 semi-infinite substrate.
The propagator of dynamically screened Coulomb interaction W (Q,ω,z,z ) for a composite structure is derived in this paper in terms of single-layer graphene response functions χ 1,2 (Q,ω) and the bulk Al 2 O 3 local dielectric function S (ω).The response functions χ 1,2 (Q,ω) are obtained using two complementary methods: an ab initio random phase approximation (RPA) method for all valence electrons in carbon, and the RPA method for graphene's π electrons described as massless Dirac fermions at zero temperature [40,41], which we call the MDF-RPA method.The propagator W is used to derive an effective 2D dielectric function (Q,ω) of the composite, giving the EELS as ∼ − [1/ (Q,ω)].Special attention is paid to comparing the ab initio RPA results for EELS with the computationally less demanding MDF-RPA results.It should be noted that, while the latter method is readily implemented in the optical limit (Q ≈ 0), the ab initio RPA method has to be treated with great care to reach this limit, which is of primary interest for optoelectronic applications.In that respect, it is thoroughly demonstrated in this work how to solve the problem of ab initio calculation of the EELS in quasi-2D crystals for a very fine Q-point mesh (i.e., when Q is smaller than the minimal value allowed by the k-point sampling) in the optical limit.It is shown that, for ω < E F and for the wave-number range Q < 0.8 nm −1 studied in this work, the EELS obtained using both methods agree quite well.The only disagreement arises in the Dirac plasmon dispersion for Q > 0.3 nm −1 and ω > E F , which establish upper bounds for the range of confidence in using the MDF-RPA method to describe nonlocal effects in the optical response of graphene.
Moreover, dispersion relations of the coupled plasmon/phonon modes for different compositions of our sandwich structure are calculated using both the optical limit of the MDF-RPA method [42], which we call the "optical MDF-RPA" method, and the (analytical) Drude model, and are compared with the EELS data from both the ab initio RPA and the full "MDF-RPA" calculations.It is shown that for Q < 0.15 nm −1 the "optical MDF-RPA" dispersion relations agree extremely well with the EELS intensity patterns from both types of the RPA calculations.It is further demonstrated that even the simple Drude model agrees very well with the EELS intensity patterns for Q < 0.08 nm −1 , which enters the THz frequency range of practical interest, thereby reaffirming the utility of this popular modeling approach in graphene-based photonics.
In Sec.II, we present the methodology used to calculate dynamically screened Coulomb interaction W (Q,ω,z,z ) and the effective 2D dielectric function (Q,ω) of a gr 1 -Al 2 O 3 -gr 2 composite.In Sec.III, we analyze the EELS and the plasmon/phonon dispersion relations in the reduced (Q < 0.08 nm −1 ) and extended (Q < 0.8 nm −1 ) wave-number intervals for several composites: single-layer graphene, doublelayer graphene, as well as symmetrical (sy 1 = sy 2 ) and asymmetrical (sy 1 = sy 2 ) composites sy 1 -Al 2 O 3 -sy 2 with one and two graphene layers.In Sec.IV, we present our concluding remarks.

A. Modeling the Al 2 O 3 slab
A possible composite structure of the gr-Al 2 O 3 -gr sandwich is shown in Fig. 1(a).The Al 2 O 3 slab is cut along (1,1,1) planes, so that its thickness is a − 2h nm, while the graphene layers occupy the z = ∓a/2 planes, and are thus located a distance h = 0.32 nm away from the nearest aluminum planes.Taking a = 5 nm, the unit cell for such a huge nanostructure consists of 2878 atoms, so it is impossible to provide full ab initio ground state and structural optimization calculation.Moreover, the ab initio calculation of the response function would be even more tedious, so we seek some approximations, which would still enable the response function calculations.The easiest (and probably the best) approximation is to treat the Al 2 O 3 slab as homogeneous dielectric described by some local dielectric function S (ω), while the two graphene layers can be considered as purely 2D systems described by the response functions χ 1,2 (Q,ω), as sketched in Fig. 1(b).For simplicity, we also take h = 0, so that the dielectric layer completely fills the space between graphene layers; it was found that this approximation has negligible effects on the coupled modes.
The dielectric properties (or the dynamic response) of the bulk ionic crystals in the long-wavelength limit can be determined in terms of their optical phonons in point.More complex polar crystals such as Al 2 O 3 posses a multitude of different optical phonons with different symmetries and polarizations.However, here we suppose that the polarization of the bulk crystal mainly comes from excitations of the two optical modes, which have the largest oscillator strengths, as proposed in Refs.[17,43].Therefore we assume that the bulk dielectric function of Al 2 O 3 is [44] where ∞ ox , i ox , and 0 ox are the optical, intermediate, and static permittivities, ω TO1 and ω TO2 are the first and second transverse optical (TO) angular frequencies (with ω TO1 < ω TO2 ), and γ TO1 and γ TO2 are the damping rates of the corresponding TO phonons.

B. Dynamic response functions of graphene layers
If the graphene layers are considered as fully 2D systems, their nonlocal independent-electron response functions may be written as [3] where Q is the momentum transfer vector parallel to the x-y plane, and z 1 = −a/2 and z 2 = a/2.Here, the χ 0 1,2 (Q,ω) functions are calculated from first principles as where the 3D Fourier transform of the independent-electron response function is given by with f 1,2 nK = [e (E nK −E 1,2 F )/kT + 1] −1 being the Fermi-Dirac distributions at the temperature T .The charge vertices in (4) have the form where G = (G ,G z ) are the 3D reciprocal lattice vectors and r = (ρ,z) is a 3D position vector.Integration in ( 5) is performed over the normalization volume = S × L, where S is the normalization surface and L is superlattice constant in the z direction (separation between graphene layers in a superlattice arrangement).The plane-wave expansion of the wave function has the form where the coefficients C nK are obtained by solving the local density approximation Kohn-Sham (LDA-KS) equations selfconsistently, as will be discussed in Sec.III A. Once the 2D response function χ 0 i (Q,ω) of the ith graphene layer is determined, a screened response function of that layer in free space may be calculated in the RPA simply as where v Q = 2π Q represents a 2D Fourier transform of the bare, "intrasystem" Coulomb interaction.

C. Optical limit
The above ab initio calculation of graphene response functions χ 1,2 is straightforward, but it is not sufficient if we want to investigate hybridization between the Dirac plasmon and the FK phonons in a dielectric surface.Namely, due to the very low energy of the FK phonons (∼50 meV), the crossing with the Dirac plasmon dispersion relation occurs for very small wave numbers, on the order of Q = 0.01 nm −1 .On the other hand, even for very dense K-point mesh sampling, such as, for example, 601 × 601 × 1 used in this calculation, the minimum transfer wave number Q that can be reached (e.g., Q = 0.049 nm −1 in this calculation) is still larger than the FK phonon-Dirac plasmon crossing wave number.Therefore we have to find a way to calculate χ 1,2 (Q,ω) for a finer Q-point mesh in the optical limit, i.e., when Q ≈ 0. One possibility is that, instead of calculating the response functions χ 0 i (Q,ω), we calculate the optical conductivities σ i (ω), which can be obtained in the strict Q = 0 limit.
The optical conductivity in graphene may be written as [4] where is the intraband, or Drude conductivity, and where represents the effective number of charge carriers.The interband conductivity is where the current vertices are given by and In the optical limit, Q ≈ 0, the independent-electron response functions can be written in terms of optical conductivities (7) as [45] Needless to say, this approximation for graphene's response function should work extremely well in the technologically interesting range of frequencies from the mid-IR (Q ∼ 0.001 nm −1 ) down to THz (Q ∼ 10 −5 nm −1 ).When expressions for the conductivity given in Eqs. ( 7)- (12) are used in Eqs. ( 13) and ( 6) constitutes what we call "optical ab initio RPA" method.

D. Screened Coulomb interaction
Once we have determined the screened response functions of graphene layers, χ 1,2 , the next step is to determine the propagator of the screened Coulomb interaction throughout our sandwich structure, W (Q,ω,z,z ), where z and z are the field point and the source point of a unit charge, respectively.This function will be obtained by two equivalent methods: a Feynman diagrammatic technique and the Dyson-Schwinger equation for the Green's function of the Poisson equation.

Diagrammatic technique
We first consider the screened Coulomb interaction in the system consisting of just one graphene layer in vacuum, placed at z = a/2, and described by the response function χ 2 (Q,ω), as shown in Fig. 1(b).Suppose that a unit point charge is placed somewhere in the region z > a/2 and that it fluctuates as cos(ωt).Then, the 2D Fourier transform of the screened Coulomb potential at the point z > a/2 may be written as (14) where represents the propagator of the surface induced Coulomb interaction.
Let us now examine the situation where a second graphene layer, described by the response function χ 1 (Q,ω), is introduced at z = −a/2, as shown in Fig. 1(b), with the space between the two layers being vacuum.If the external point charge is placed in the region z > a/2, then from the many-body perturbation theory point of view, it can induce charge density fluctuations (in two graphene sheets) in four different ways [46].The first class of processes corresponds to a situation where the charge density fluctuation is created and annihilated in the first graphene layer.This class of processes gives the surface induced Coulomb interaction D 11 , as shown by the Feynman diagrams in Fig. 2(a).To the lowest order, the external charge creates a charge density fluctuation in the first sheet, which annihilates in the same sheet and produces the Coulomb field.In the next, higher order, the external charge creates charge fluctuation in the first sheet, which annihilates and creates a charge fluctuation in the second sheet, which annihilates and creates a charge in the first sheet, which finally annihilates in the same sheet and produces the field.Following the Feynman diagrams in Fig. 2(a), the interaction D 11 may be written as where v 12 (Q) = v Q e −Qa represents a 2D Fourier transform of the bare "intersystem" Coulomb interaction.The second class of processes corresponds to a situation where the charge density fluctuation is created and annihilated in the second sheet.In analogy with the term D 11 , this class of processes gives the surface induced Coulomb interaction D 22 , which may be written following the Feynman diagrams in Fig. 2(a) as where v 11 (Q) = v Q represents a 2D Fourier transform of the bare "intrasystem" Coulomb interaction in the first graphene sheet.
FIG. 2. The Feynman diagrams for interaction between charge density fluctuations in two graphene layers described by response functions χ 1 (Q,ω) and χ 2 (Q,ω) in two cases, when the region between the layers is (a) vacuum and (b) dielectric slab described by the local dielectric function S (ω).The black dashed lines represent the bare "intrasystem" v 11 and "intersystem" v 12 Coulomb interactions, whereas the blue wavy lines represent the screened intrasystem ṽ11 and intersystem ṽ12 Coulomb interactions in the presence of dielectric slab.
The third and fourth classes of processes correspond to situations where the charge density fluctuation is created in the first sheet and is annihilated in the second sheet and vice versa.These classes of processes give the surface induced Coulomb interactions D 12 and D 21 , which may be written following the Feynman diagrams in Fig. 2(a) as Finally, the total surface induced Coulomb interaction is the sum of four terms: The expression for the Fourier transform of the screened Coulomb interaction (for z,z > a/2) remains unchanged, i.e., it has the form given in Eq. ( 14) with D(Q,ω) derived in Eq. (19).The introduction of a dielectric slab with S between graphene layers only slightly modifies the previous formulation.The only change to be made is that the bare "intrasystem" and "intersystem" Coulomb interactions v 11 and v 12 should be replaced by the corresponding interactions "screened" by the dielectric slab, which have the following analytical form: where From the Feynman diagrams' point of view, this causes the black dashed lines in Fig. 2(a) to be replaced by the blue wavy lines, as shown in Fig. 2(b).This also causes that the RPA screened response functions become functions of the screened intrasystem Coulomb interaction ṽ11 , Therefore, because all quantities that enter in the total surface induced Coulomb interactions, Eqs. ( 16)-( 19), become renormalized, D(Q,ω) itself becomes renormalized and should be rewritten as where D11 , D22 , D12 , and D21 have the same form as those given in Eqs. ( 16), (17), and (18), except that v 11 , v 12 , and χ 1,2 should be replaced by ṽ11 , ṽ12 , and χ1,2 .Finally, the only expression that changes is the Fourier transform of the screened Coulomb interaction, which now has the form (for z,z > a/2) where

Dyson-Schwinger equation
The screened Coulomb interaction W (Q,ω,z,z ) is equal (in atomic units) to the Green's function (GF) of the Poisson equation, which may be easily derived for a layered structure exhibiting translational invariance in directions parallel to graphene layers by using the standard electrostatic matching conditions at the interfaces between different dielectric regions [47,48].Accordingly, this approach is quite efficient for rather complex layered structures, where inclusion of graphene layers at arbitrary positions may be described by the Dyson-Schwinger (DS) equation for W (Q,ω,z,z ) of the form (25) where we have dropped Q and ω to simplify the notation.Here, W 0 (z,z ) ≡ W 0 (Q,ω,z,z ) denotes the GF for the same layered structure without graphene, whereas the effects of N graphene sheets, which are placed in the planes z = z i , are described by the "perturbation" V(z) ≡ V(Q,ω,z) given by where ) is the independent-electron response function of the ith graphene sheet.
Considering our sandwich structure with N = 2 graphene sheets placed at z 1 = −a/2 and z 2 = a/2, we obtain from Eqs. ( 25) and ( 26) We next set z = z 1 and z = z 2 in Eq. ( 27) to obtain a system of algebraic equations for the unknown values of W (z 1 ,z ) and 29) where we have defined the "coefficients" ṽij ≡ W 0 (z i ,z j ).For an oxide film described by a (relative) dielectric function S (ω), which occupies region −a/2 z a/2 and is surrounded by vacuum, one obtains expressions for ṽ11 = ṽ22 and ṽ12 = ṽ21 given in Eqs. ( 20) and ( 21), respectively.Using those expressions in Eqs. ( 28) and ( 29), one can obtain W (z 1 ,z ) and W (z 2 ,z ), which need to be substituted in the right-hand side of Eq. ( 27) to yield a final expression for the full GF, W (z,z ), given in Eq. ( 24) for z,z a/2.

E. Effective 2D dielectric function in gr-Al 2 O 3 -gr composite
Having formulated screened Coulomb interaction enables us to define an effective 2D dielectric function and the loss function for our sandwich structure.Suppose that the structure is probed by REELS so that the incident electron trajectory is localized in the region z > a 2 with the reflection point at z = a 2 .Referring to Eq. ( 24), the 2D Fourier transform of the induced Coulomb interaction in the z = a 2 plane may be written as where we introduce the effective 2D dielectric function (Q,ω).Then, it can be shown [39] that the probability density for the parallel momentum transfer Q and the energy loss ω of the reflected electron in REELS is proportional to the imaginary part of the induced potential in Eq. (30), which is called the energy loss function, or the EELS.The loss function in Eq. ( 31) may be calculated from expressions for the effective 2D dielectric function in different systems.For example, when the system consists of just one (second) graphene layer in vacuum at z = a 2 , we have whereas for the system consisting of two graphene layers placed in vacuum at z = ∓ a 2 , we have where D(Q,ω) is given in Eqs. ( 16)- (19).Finally, when the considered system is a gr-Al 2 O 3 -gr sandwich, the effective 2D dielectric function is with ṽ11 and D(Q,ω) given in Eqs. ( 20) and ( 23), respectively.After some algebra, this expression may be written in a more explicit form as

F. Dispersion relations
One can obtain dispersion relations for the collective modes in a sandwich structure with two graphene layers by setting to zero all damping rates and by solving the equation (Q,ω) = 0 in regions of the (Q,ω) plane where the imaginary parts of χ 0 i (Q,ω) are (at least approximately) zero.In the special case of graphene layers with equal polarization functions, χ 0 1 = χ 0 2 = χ 0 , Eq. ( 35) gives rise to two decoupled equations for dispersion relations for symmetric and antisymmetric coupling, respectively, between their Dirac plasmons, We finally note that, in order to cast Eqs. ( 36) and (37) in the optical limit, one may use χ 0 = −iQ 2 σ (ω)/ω along with analytical expressions for the conductivity of doped graphene, σ (ω) = σ intra (ω) + σ inter (ω), which were obtained in Ref. [42] by treating the π electron bands in the Dirac cone approximation.In the limits of zero damping and zero temperature, those expressions yield the intraband contribution in a Drude form, σ intra (ω) = iE F /(πω), while the interband contribution is given by [42] where (ω) is a Heaviside unit step function.Using the above result for graphene conductivity in the optical limit [42] along with the bulk dielectric function S (ω) for Al 2 O 3 given by Eq. ( 1) with zero damping rates, one may solve Eqs.(36) and (37) to obtain the dispersion relations for hybridized plasmon-phonon modes in the limit of long wavelengths.Interestingly, when only the Drude contribution σ intra (ω) ∝ i/ω is retained for graphene, which is common practice in photonic applications of graphene, Eqs. ( 36) and ( 37) result in a pair of cubic equations in the square of frequency, which give the dispersion relations of all modes in analytical form.

A. Computational details
The first part of our ab initio calculations consists of determining the KS ground state of single-layer graphene and the corresponding wave functions φ nK (ρ,z) and energies E n (K).For graphene's unit-cell constant we use the experimental value of a g = 0.245 nm [49], while for the superlattice unit-cell constant (separation between the periodic replicas of graphene layers) we take L = 5a g .For calculating the KS wave functions and energies we use a plane-wave, self-consistent field DFT code (PWSCF) within the QUANTUM ESPRESSO (QE) package [50].The core-electron interaction is approximated by the norm-conserving pseudopotentials [51], and the exchange correlation (XC) potential by the Perdew-Zunger local density approximation (LDA) [52].To calculate the ground-state electronic density, we use a 21 × 21 × 1 Monkhorst-Pack K-point mesh [53] of the first Brillouin zone (BZ), and for the plane-wave cut-off energy we choose 50 Ry.The second part of our ab initio calculations consists of determining the independent-electron response function ( 4) and the optical conductivity ( 7)- (10).In order to achieve better resolution in the long-wavelength limit (Q ≈ 0) and the low energy (ω ≈ 0) limit, the response function [( 4) and ( 5)] and the conductivity ( 7)-( 12) are evaluated from the wave functions φ nK (r) and energies E n (K) calculated for the 601 × 601 × 1 Monkhorst-Pack K-point mesh, which corresponds to 361801 K points in the first Brillouin zone (1BZ).The band summations (n,m) in ( 4), ( 9) and ( 10 The results for EELS obtained from the above ab initio RPA method for structures including one or two graphene layers are compared with the results based on analytical expressions for the independent-electron response function χ 0 i (Q,ω) involving just the π electron bands in the Dirac cone approximation at zero temperature [40,41], with the effects of damping included via the Mermin procedure [39].We call such a method an 'MDF-RPA' approximation.While the ab initio RPA method is expected to describe EELS of graphene in a broad range of frequencies up to FUV, the MDF-RPA method should provide good account of nonlocal effects in doped graphene at frequencies in the THz to IR range, which is of primary technological interest.
Nonlocal effects in dynamic response of graphene are often deemed unimportant in photonic applications, where optical response of graphene is well described by a frequencydependent sheet conductivity.In that respect, besides the optical limit of the ab initio RPA method given in Eqs. ( 7)-( 12), we also use the analytical expressions for graphene's conductivity obtained in Ref. [42] for the π electron bands in the Dirac cone approximation at zero temperature, which we call the "optical MDF-RPA" approximation.In particular, those expressions are used to obtain the dispersion relations of hybridized modes in our sandwich structure in the limit of zero damping, see Eqs. ( 36)- (38).

B. EELS for single-and double-layer graphene in vacuum
In this section we shall first present the spectra of electronic excitations, or the EELS in a free-standing, doped (E F = 1 eV) single-layer graphene.Then, we shall briefly analyze the EELS in a double-layer graphene, which consists of two graphene layers in vacuum, separated by a distance a = 5 nm, and doped with equal Fermi energies, E F = 200 meV.Even though the EELS in pristine and doped graphene are widely investigated experimentally and theoretically [1][2][3][4]6], here we provide high accuracy calculations of the long-wavelength and the lowenergy EELS.Special attention is paid to exploring the means to achieve continuous matching between the spectra in the long wavelength, Q ≈ 0, limit calculated using the optical ab initio RPA method [with χ 0 i calculated using Eq. ( 13)] and the spectra calculated using the full ab initio RPA method [with χ 0 i calculated using Eq. ( 3)]. Figure 3 shows the energy loss spectrum (31) in a singlelayer doped graphene, where E F = 1 eV.The low-energy part of the spectrum is dominated by the Dirac plasmon, with dispersion that starts as square root and then merges the upper  8) and (b) the "full" optical conductivity σ = σ intra + σ inter (8,10).Vertical black lines denote the boundary wave number Q b = 0.34 nm −1 .The thin white lines represent the electron-hole excitations edges in gr(1 eV), as described in Fig. 3. edge of the π * ↔ π * intraband continuum, ω = v F Q, denoted by a white line.The higher-energy region is dominated by a broad π plasmon, starting at ω = 4 eV for Q = 0.The black vertical line denotes the boundary wave number Q b , which separates the energy loss spectrum calculated using the optical ab initio RPA approximation (left region) and using the full ab initio RPA method (right region).It can be noticed that, for the boundary wave number chosen as Q b = 0.34 nm −1 , matching between the optical ab initio RPA and the full ab initio RPA results is perfectly smooth.However, we have found that, as Q b increases the "matching" becomes worse and for, e.g., Q b ≈ 1.0 nm −1 , the two methods yield quite mismatched spectra.These conclusions are very important, suggesting that the extremely demanding ab initio calculations of the EELS for optically small wave numbers can be obtained by a "shortcut" using the optical ab initio RPA method, which only requires conductivity calculations with vanishing wave number, Q = 0.
Figures 4(a separates regions where the optical approximation and the full ab initio RPA methods are used.In order to illustrate how the accuracy of the optical conductivity calculation influences the Dirac plasmon dispersion relation, in Fig. 4(a), the EELS at Q < Q b is calculated using only the Drude optical conductivity σ = σ intra from the optical ab initio RPA, Eq. ( 8), whereas in Fig. 4(b) we use the complete optical conductivity σ = σ intra + σ inter from the optical ab initio RPA, Eqs. ( 8) and (10).It may be seen that, if the optical conductivity includes only the intra term, the Dirac plasmon does not match the full ab initio RPA result at Q b .However, if both the intra and inter terms are included in the optical conductivity, the Dirac plasmon matches the full ab initio RPA spectra quite smoothly.
Figures 5(a FIG. 5.The EELS for two graphene layers calculated using: (a) the optical ab initio RPA method in the whole range of wave numbers, (b) the optical ab initio RPA method for Q < Q b , and the full ab initio RPA method for Q > Q b .The solid magenta and white lines show the dispersion relations of odd (ω − ) and even (ω + ) Dirac plasmons, respectively, calculated using the optical MDF-RPA model.The dotted magenta and white lines show the dispersion relations of odd (ω 0 − ) and even (ω 0 + ) Dirac plasmons, respectively, calculated using only the Drude model (8).The boundary wave number is Q b = 0.098 nm −1 , the separation between graphene layers is a = 5 nm and both graphene layers are doped such that E F = 200 meV.The thin white lines represent the electron-hole excitations edges in gr(200 meV), as described in Fig. 3.
of the optical ab initio RPA method, the EELS in Fig. 5(a) is calculated by using this approximation in the full range of the shown wave numbers, whereas in Fig. 5(b) the EELS is calculated using the optical ab initio RPA method for Q < Q b and the full ab initio RPA method for Q > Q b .Here, the boundary wave number is chosen as Q b = 0.098 nm −1 .
We see two plasmon branches, which are a consequence of hybridization between the Dirac plasmons in graphene layers.Because the potential, which produces the 2D plasmons decays as ∼e −Qa , a "cutoff" wave number below which two plasmons still interact may be estimated as Q C ∼ 1/a.Therefore, for Q > Q C , the plasmons interact weakly, and the two branches become almost degenerate.As can be seen in Fig. 5, this estimate works quite well, as the Dirac plasmon splitting occurs below Q C ∼ 0.2 nm −1 , while above that value the two plasmons tend to merge.Because the system is symmetric (along the z axis), it supports modes with well-defined parity, i.e., the potential produced by such modes is either symmetric (even) or antisymmetric (odd) function of z with respect to the origin at z = 0. Consequently, this means that the symmetric or antisymmetric modes can only be excited by a symmetrically or antisymmetrically designed external perturbation.For example, the system may be imagined as being excited by two symmetrically placed point charges, one at z = −a/2 − 0 + that oscillates as Q 1 cos ωt, and the other at z = a/2 + 0 + that oscillates as Q 2 cos ωt (here, 0 + is infinitesimally small positive number).The lower plasmon branch, denoted by ω − , is excited by external charges of the opposite sign, Q 1 = −Q 2 = Q, which implies that those plasmons produce antisymmetric (odd) electrical potential, accompanied by the out-of-phase oscillations of the induced charges in graphene layers.On the other hand, the upper plasmon branch, denoted by ω + , is excited by two equal external charges Q 1 = Q 2 = Q, so those plasmons produce symmetric (even) electrical potential, accompanied by the in-phase oscillations of the induced charges in graphene layers [13].Accordingly, the corresponding hybridized plasmon modes in double-layer graphene with eigenfrequencies ω ± will be designated as the symmetric (or even) and the antisymmetric (or odd) modes, respectively.
Figure 5 shows that the optical ab initio RPA EELS nicely matches the full ab initio RPA EELS at Q = Q b .Moreover, Fig. 5 shows that the two methods qualitatively agree, even up to Q = Q C = 0.2 nm −1 , but for Q > 0.2 nm −1 , they become drastically different in two important aspects.First, the full ab initio RPA Dirac plasmons in Fig. 5(b) continue increasing following the upper edge of the π * ↔ π * intraband electron-hole continuum, while the optical ab initio RPA dispersions in Fig. 5(a) are deflected into that continuum.This is reasonable to expect because, in the optical limit, there is no linearly dispersing (ω = v F Q) π * ↔ π * intraband electron-hole continuum, which would push the plasmons towards higher energies.Second, the optical ab initio RPA Dirac plasmons are seen to produce well-defined resonances in the whole range of the shown Q values, while the ab initio RPA plasmons become strongly damped already at Q ≈ 0.4 nm −1 due to Landau damping in the π ↔ π * interband electron-hole continuum taking place at frequencies ω > 2E F − v F Q. On the other hand, in the optical ab initio RPA limit, Eq. ( 13), the π ↔ π * interband electron-hole continuum is dispersionless, i.e., its lower edge is at ω = 2E F , so that the Dirac plasmons are pushed down, away from the value ω = 2E F .
The solid magenta and white lines in Fig. 5 show the dispersion relations of the odd (ω − ) and even (ω + ) Dirac plasmons, respectively.Those relations were calculated by using the optical MDF-RPA model with zero damping and setting S = 1 in Eqs. ( 36) and (37), respectively.The agreement with the optical ab initio RPA EELS intensity patterns in Fig. 5(a) is excellent in the whole Q interval, while the agreement with the full ab initio RPA EELS intensity patterns in Fig. 5(b) becomes, as excepted, worse for larger wave numbers (Q 0.2 nm −1 ).The dotted magenta and white lines show the dispersion relations of the odd (labeled by ω 0 − ) and even (labeled by ω 0 + ) Dirac plasmons, respectively, calculated using only the Drude model of graphene conductivity, Eq. ( 8).The agreement of this simple model with the optical ab initio RPA and the full ab initio RPA Dirac plasmons is satisfactory up to the boundary wave number (Q < Q b ), but for Q > Q b , the Drude model drastically overestimates the energy of Dirac plasmons.

C. EELS for the gr-Al 2 O 3 -gr composite
In this section, we shall first study the EELS for a gr-Al 2 O 3 -gr sandwich structure in a reduced (Q,ω) space in order to emphasize and explore hybridization between four FK phonons in sapphire and two Dirac plasmons in graphene.The dispersion relations of such modes, calculated using only the Drude model and the optical MDF-RPA model, will also be presented.We shall then study the EELS in the same gr-Al 2 O 3 -gr sandwich using the ab initio RPA and MDF-RPA methods in an extended (Q,ω) space in order to investigate the influence of realistic band structure on the six hybridized modes and their interaction with the electron-hole continua.The separation between graphene layers (or Al 2 O 3 slab thickness) is a = 5 nm, while both graphene layers are doped such that E F = 200 meV.
Figure 6(a) shows the EELS in the reduced (Q,ω) space for a gr-Al 2 O 3 -gr sandwich structure calculated using the optical ab initio RPA method, while Fig. 6(b) shows the same spectrum obtained using the optical MDF-RPA method.It can be clearly seen that those two spectra coincide, which justifies the use of the computationally much simpler (semi-analytic) optical-MDF-RPA method in the reduced (Q,ω) space.Both EELS show intense patterns, which represent six hybridized plasmon-phonon modes.Because the system is symmetric, the six modes can be divided into a symmetric group and an antisymmetric group, each containing three modes.The black solid lines show the dispersion relations for those six modes calculated using Eq.(1) with zero damping for S and the optical MDF-RPA model with zero damping for graphene layers, whereby Eq. ( 36) yields three antisymmetric (odd) modes (ω − 1 , ω − 2 , and ω − 3 ) and Eq. ( 37) yields three symmetric (even) modes (ω + 1 , ω + 2 , and ω + 3 ).The six dispersion curves agree perfectly well with the EELS intensity patterns, which will help us determine the symmetry and the phononic/plasmonic character of the particular modes in the EELS.
The parity of the modes may be, as in the case of doublelayer graphene, probed by two fluctuating point charges, symmetrically placed at z = ±( a 2 + 0 + ), which have equal or opposite signs.The first linear, the third horizontal, and the more intense fifth parabolic branches, denoted by ω − 1 , ω − 2 , and ω − 3 , respectively, represent odd modes, i.e., the modes which produce antisymmetric electrical potential oscillations.On the other hand, the second and the fourth horizontal, as well as the most intense sixth, square-root-like branches, denoted by ω + 1 , ω + 2 , and ω + 3 , respectively, represent even modes, i.e., the modes which produce symmetric electrical potential oscillations.It is also appropriate to determine the phononic/plasmonic character of these modes, which strongly depends on the wave number Q.For larger Q values in Fig. 6, it is clear that the branches ω ± 1 and ω ± 2 are purely phononlike, i.e., they represent the polarization oscillations, which are localized in the two dielectric surfaces.On the other hand, modes ω ± 3 are plasmonlike, i.e., they represent the charge density oscillations which are localized in the two graphene layers.However, in the Q → 0 limit, strong hybridization between some of these modes takes place, giving rise to avoided crossings, which depend on the symmetry matching.For example, when the symmetric phonon branches ω + 1 and ω + 2 become close to the symmetric plasmon branch ω + 3 they strongly interact with the branch ω + 3 and they take the plasmonic character.On the other hand, the antisymmetric branches ω − 1 and ω − 2 cross the symmetric plasmon branch ω + 3 without interaction and retain their phononic behavior.Similarly, also the plasmonic branches ω + 3 and ω − 3 cross each other without interaction because of the symmetry mismatch.On the other hand, the linear branch ω − 1 does not interact, or interacts very weakly with any other modes.
The two red dotted lines in Fig. 6, denoted by ω 0± 3 represent dispersion relations of the corresponding hybridized modes calculated using the Drude model for conductivity of graphene layers, σ = σ intra , from Eq. ( 8) with zero damping.We do not show in this figure the corresponding four branches calculated with the Drude model for ω 0±  1 and ω 0± 2 because they fully overlap with the dispersion relations for ω ± 1 and ω ± 2 obtained using "complete" conductivity, σ = σ intra + σ inter , from the optical MDF-RPA.There is only a very small disagreement between the asymmetric modes, where ω 0− 3 negligibly overestimates the branch ω − 3 for the highest wave numbers, Q ≈ 0.08 nm −1 .The largest deviation is shown for the symmetric branch ω 0+  3 , which starts overestimating the branch ω + 3 already at Q ≈ 0.01 nm −1 , and continues to overestimate it with increasing wave number, giving even up to about 20% larger frequency for Q ≈ 0.08 nm −1 .This disagreement between the ω 0+ 3 and ω + 3 curves results from the lack of edge of the interband π ↔ π * electron-hole transitions in the Drude model, which would push the highest plasmon branch ω 0+ 3 towards lower frequencies, as it happens to the branch ω + 3 in the complete optical MDF-RPA model.However, the simple Drude approximation gives good dispersions of practically all other hybridized modes in a gr-Al 2 O 3 -gr sandwich, thereby covering the long-wavelength and the THz limits of interest in plasmonic applications.
Figure 7 shows the EELS in the gr-Al 2 O 3 -gr sandwich in an extended (Q,ω) region, with the EELS in Fig. 7(a) calculated using the MDF-RPA method and in Fig. 7(b) using the ab initio RPA method.The dotted lines show dispersion relations, which were calculated in the same way as in Fig. 6.In both panels, the EELS are dominated by the highest plasmon branches, with the peak positions that agree well with the corresponding dispersion relations ω ± 3 for Q < 0.15 nm −1 .For larger Q values, the EELS peak positions from the MDF-RPA method start to increase much faster than the corresponding dispersion relations for ω ± 3 .On the other hand, it is surprising to see that those dispersion relations agree quite well with the peak 2 and ω + 3 ) calculated using the optical MDF-RPA method.The thin white lines represent the electron-hole excitations edges in gr(200 meV), as described in Fig. 3. positions in the EELS from the ab initio RPA method, even up to the edge of the interband π ↔ π * electron-hole transitions, ω = 2E F − v F Q. One also sees that the EELS peak positions for ω ± 3 from both RPA methods merge for Q ≈ 0.3 nm −1 , while the corresponding dispersions curves become degenerate only for quite large wave numbers, Q ≈ 0.6 nm −1 .
One can further see in Fig. 7 that the EELS from both RPA methods agree quite well with each other for Q < 0.3 nm −1 .However, for larger wave numbers Q, the peak positions in the EELS from the MDF-RPA method corresponding to the ω ± 3 modes continue to increase in a direction parallel to the upper edge of the intraband π * ↔ π * continuum (ω = v F Q, resulting from the Dirac cone approximation), while in the ab initio RPA method those peak positions move to lower energies, and they cross the upper edge of the intraband π * ↔ π * continuum at Q ≈ 0.45 nm −1 .The reason for this disagreement is probably due to the use of the Dirac cone approximation in the MDF-RPA method, as well as due to the lack of the π plasmon in the MDF-RPA response function, which even though energetically is positioned very high, polarizes (screens) the system and pushes the Dirac plasmons towards lower energies in the ab initio RPA method.Moreover, besides the transitions within and between the π and π * bands, the ab initio RPA response function in Eq. (4) also includes transitions between higher occupied and unoccupied bands (e.g., π ↔ σ * and σ ↔ σ * transitions, [6]) which represent additional polarizing mechanisms that may affect the ω ± 3 modes.In addition, the ab initio RPA method allows a more realistic calculation of the charge vertices in Eq. ( 5), with an accuracy that becomes especially important for shorter wave numbers.
Finally, considering the dispersion relations for increasing wave numbers Q, one can see in Fig. 7 that all three even/odd pairs of curves ω ± i with i = 1,2,3 begin to merge due to increasing degeneracy.While merging of the plasmonlike curves ω ± 3 occurs at quite large Q values, one sees that the dispersion curves ω ± 1 and ω ± 2 become degenerate already at Q > Q C = 0.2 nm −1 , and they settle at the bulk values of the first and second TO phonon frequencies, ω ± 1 ≈ ω TO1 = 48 meV and ω ± 2 ≈ ω TO2 = 71 meV, respectively.However, the EELS from both RPA methods exhibit in Fig. 7 two faint horizontal lines at energies around 52 and 84 meV for Q > 0.4 − 0.5 nm −1 , which are slightly higher in comparison with the corresponding values reached by the phonon dispersion curves, ω ± 1,2 ≈ ω TO1,2 .Since those faint lines occur at large wave numbers, one may estimate their energies by setting Qa 1 in Eqs. ( 36) and ( 37), either of which then gives a dispersion relation 1 + S (ω) − 2v Q χ 0 (Q,ω) = 0, which describes single graphene layer on the surface of a semi-infinite region filled with Al 2 O 3 .One is further tempted, at such "low" frequencies, to invoke static approximation for the response function of graphene layers, χ 0 (Q,ω) ≈ χ stat (Q), and moreover to use its large-Q form, χ stat = −Q/(4v F ), which is strictly valid for Q k F [40,41].As a result, we obtain a very simple dispersion relation, 1 + S (ω) + π/v F = 0 (recall that v F ≈ 0.46 in atomic units), giving two phonon eigenfrequencies of about 54 and 86 meV, which are in close agreement with the positions of the two faint horizontal lines in the EELS in Fig. 7.

D. EELS for the sy 1 -Al 2 O 3 -sy 2 composite
In this section, we shall concentrate on the study of the EELS in asymmetric sy 1 -Al 2 O 3 -sy 2 composite systems, where sy i represents either an independently doped single-layer graphene, sy i =gr(E F,i ), or simply vacuum, sy i =vacuum.As in Sec.III C, we shall first study the EELS and dispersion relations of hybridized modes in a reduced (optical) range, followed by a study in an extended (Q,ω) range.
Figure 8 shows the EELS in the reduced (Q,ω) space for the following structures: (a) vacuum-Al 2 O 3 -gr(200 meV), (b) gr(200 meV)-Al 2 O 3 -vacuum, (c) gr(0 meV)-Al 2 O 3 -gr(200 meV), and (d) gr(200 meV)-Al 2 O 3 -gr(0 meV), calculated using the optical ab initio RPA method.Because now one side of our sandwich structure is either vacuum or pristine graphene (E F = 0), which does not support Dirac plasmons, there are five prominent intensity patterns in all panels of Fig. 8, which correspond to five branches, resulting from hybridization between four FK phonons in the dielectric surfaces and the Dirac plasmon in one doped graphene.The black lines in all four panels show the dispersion relations of the five modes (ω 1 , ω 2 , ω 3 , ω 4 , and ω 5 ) obtained from the equation (Q,ω) = 0 by using Eq. ( 35) where we set one of the response functions of graphene layers, χ 0 i , to zero and calculate the other response function using the optical MDF-RPA method for doped graphene with zero damping [42].Because now the symmetry of the system is broken, the parity is not a good quantum number, and the modes can no longer be classified as symmetric or antisymmetric.For shorter wavelengths, the horizontal branches ω 1 , ω 2 , ω 3 , and ω 4 are phononlike, while the "square-root" branch ω 5 is plasmonlike.In the longwavelength limit (Q → 0), a strong hybridization (indicated by avoided crossings) between those modes occurs, such that the modes ω 1 , ω 3 , and ω 4 take plasmonic character and the branch ω 5 takes phononic character.The phonon branch ω 2 has very low intensity, so just from its dispersion relation (which does not exhibit avoided crossing) we conclude that it does not hybridize with the plasmon.The spectrum shown in Fig. 8(a) is dominated by a very intense plasmonlike branch.This is because the external probing electron, which pumps energy into the system, is placed at the point z = a/2 + 0 + , right next to the graphene layer gr(200 meV) in the z = a/2 plane, which supports a strong absorber-Dirac plasmon. .This is a consequence of the fact that the (Q,ω) phase space for the allowed interband π ↔ π * electron-hole excitations in gr(0 eV) covers the region ω > v F Q, which causes strong Landau damping of all hybridized plasmon/phonon modes shown in Figs.8(a) and 8(b).Moreover, the EELS intensity patterns in Figs.8(c) and 8(d) no longer agree so well with the dispersion relations ω 1 , ω 3 , ω 4 , and ω 5 .This is because the interband π ↔ π * electron-hole excitations in gr(0 eV) represent a weak but finite polarization mechanism, which "screens" the existing plasmon/phonon modes within the optical ab initio RPA method, while that mechanism is not included in the optical MDF-RPA model, which treats gr(0 eV) the same way as a vacuum layer.
The red dotted lines in Fig. 8 show the dispersion relations of the two highest modes calculated using only the Drude part of the conductivity in the optical MDF-RPA model, Eq. ( 8), with zero damping, which are labeled by ω 0 4 and ω 0 5 .The corresponding dispersion relations ω 0 1 , ω 0 2 , and ω 0 3 are not shown because they perfectly overlap with the dispersion relations ω 1 , ω 2 , and ω 3 , obtained from the optical MDF-RPA model with complete conductivity.The agreement between ω 0 4 and ω 4 is almost perfect, whereas ω 0 5 slightly overestimates ω 5 , but only for large wave numbers (Q > 0.02 nm −1 ), which is due to the lack of interband π ↔ π * transitions when only the Drude model is used to describe gr(200 meV).
Figure 9 shows EELS in the extended (Q,ω) space for (a) vacuum-Al 2 O 3 -gr(200 meV), (b) gr(200 meV)-Al 2 O 3vacuum, (c) gr(0)-Al 2 O 3 -gr(200 meV), and (d) gr(200 meV)-Al 2 O 3 -gr(0 meV), calculated using the MDF-RPA method.The results obtained using the ab initio RPA method are not shown because they are in excellent agreement with the results from the MDF-RPA method.There is only small quantitative disagreement that appears for larger wave numbers (Q > 0.3 nm −1 ), and it mostly affects the highest plasmon branch ω 5 (in cases when it still retains some intensity), in such a manner that its peak intensity appears at somewhat lower energies in the ab initio RPA method when compared with the MDF-RPA result.This is similar to the disagreement between the two RPA methods seen in Figs.7(a) and 7(b) for the branches ω ± 3 , where the mechanism responsible for such disagreement was discussed in some detail.
Figures 9(a) and 9(b) mostly confirm the conclusions of Figs.8(a) and 8(b).When doped graphene is placed next to the external probing electron, the spectra are dominated by the Dirac plasmon ω 5 , and when the bare dielectric surface is exposed to the external probing electron, the spectra are dominated by the topmost phonon ω 4 .However, the extended spectral range reveals some intriguing features in the spectra from the MDF-RPA method.Namely, Fig. 9(a) shows that the Dirac plasmon is, after entering the π ↔ π * interband transitions (i.e., after crossing the line ω = 2E F − v F Q) still quite intense, whereas for larger wave numbers (Q > 0.4 nm −1 ) it becomes broader, but still remains a well defined eigenmode.On the other hand, Fig. 9(b) shows that the intensity of the Dirac plasmon rapidly decreases with increasing Q, so that soon after crossing the π ↔ π * edge (ω = 2E F − v F Q) it drops to zero.There are two reasons for this behavior, one of a quantum mechanical origin due to Landau damping in the continuum of interband π ↔ π * transitions, and the other of a purely electrostatic origin.Namely, doped graphene layer in Fig. 9(a) is placed at z = a/2 and is probed by an external electron right next to it, while in Fig. 9(b), graphene layer is placed at the opposite side of the sandwich structure, a distance a from the external probing electron.In the latter case, the electric field produced by the Dirac plasmon in doped graphene at the opposite side is reduced by an additional factor e −2Qa , which acts together with Landau damping to attenuate the intensity of the peak corresponding to the ω 5 mode in Fig. 9(b).These conclusions also suggest that the polarization of the Al 2 O 3 slab does not screen enough the electric field due to the Dirac plasmon, i.e., it freely penetrates through the Al 2 O 3 slab.
The white dotted lines in Fig. 9 show dispersion relations of five modes (ω 1 , ω 2 , ω 3 , ω 4 and ω 5 ) obtained by solving the equation (Q,ω) = 0 using the optical MDF-RPA model in Eq. ( 35) with zero damping rates.It should be noted that in all panels of Fig. 9 the dispersion curves ω 1 and ω 3 settle for large Q values at about the bulk TO phonon frequencies, ω 1 ≈ ω TO1 = 48 meV and ω 3 ≈ ω TO2 = 71 meV, respectively.At the same time, it may be shown that the dispersion curves ω 2 and ω 4 settle for large Q values at about the corresponding FK phonon frequencies for the bare Al 2 O 3 surface, ω 2 ≈ ω FK1 ≈ 56 meV and ω 4 ≈ ω FK2 ≈ 110 meV, respectively, which are obtained by solving the dispersion relation 1 + S (ω) = 0 with zero damping.
It is interesting that the spectrum from the MDF-RPA method for large wave numbers Q shows in Fig. 9(a) a phononlike mode as an intense horizontal line, placed between the dispersion curves ω 3 and ω 4 , which is accompanied by a much weaker line placed between the dispersion curves ω 1 and ω 2 .On the other hand, one sees that the spectrum in Fig. 9(b) exhibits a very intense horizontal line, which fits perfectly well the upper FK phonon mode ω 4 ≈ ω FK2 ≈ 110 meV for all wave numbers Q > 0.03 nm −1 [see also Fig. 8(b)], and a much weaker horizontal line, which fits the lower FK phonon mode ω 2 ≈ ω FK1 ≈ 56 meV.
The energies of phononlike lines in the spectra in both Figs.9(a) and 9(b) may be understood by considering large wave numbers, such that Qa 1, wherefrom Eq. ( 35) gives a dispersion relation 1 + S (ω) − 2v Q χ 0 2 (Q,ω) = 0 describing a single graphene layer with the response function χ 0 2 , sitting on the surface of a semi-infinite region filled with Al 2 O 3 .For the case in Fig. 9(a), we may again approximate the graphene response function at large wave numbers and low frequencies by χ 0 2 ≈ χ stat = −Q/(4v F ), as in the discussion of Fig. 7.The resulting approximate dispersion relation, 1 + S (ω) + π/v F = 0 gives two phonon frequencies of about 54 and 86 meV, which closely match the positions of the two horizontal lines in Fig. 9(a) for Q > 0.4 nm −1 .For the case in Fig. 9(b), we simply obtain the dispersion relation 1 + S (ω) = 0 for two FK phonons on the surface of a semi-infinite Al 2 O 3 , ω FK1 ≈ 56 meV and ω FK2 ≈ 110 meV, which closely match the energies of two horizontal lines in that figure.Finally, we note that the large difference in intensity of those lines comes from the fact that the spectral weights of two narrow lines in the corresponding loss function, − {1/[1 + S (ω)]}, obtained with vanishing damping, result in the ratio ∼1 : 22 for the FK1:FK2 phonons.
Figures 9(c) and 9(d) exhibit quite strong impact of the replacement of the vacuum layer by a pristine graphene layer.Its effect is mainly manifested as additional Landau damping in the region ω > v F Q due to the interband π ↔ π * transitions in pristine graphene, so that the Dirac plasmon originating from the doped graphene layer exhibits reduced intensity in comparison to Figs. 9(a) and 9(b).It is interesting that this damping of the Dirac plasmon is not so strong in Fig. 9(c) when pristine graphene is placed at the opposite side from the probing electron.On the other hand, intrinsic graphene in Fig. 9(d) practically destroys the Dirac plasmon for Q > 0.04 nm −1 [see also Fig. 8(d)] when the doped graphene is placed at the opposite side from the probing electron.
Considering the two phononlike horizontal lines, which are present in the spectra in both Figs.9(c) and 9(d) for large wave numbers, one notices that they lie at approximately the same energies in both figures, albeit with different spectral weights.As in Fig. 9(a), the stronger line in Figs.9(c) and 9(d) is positioned between the phonon dispersions ω 3 and ω 4 , and the weaker line between ω 1 and ω 2 .Again as in Fig. 9(a), positions of the two horizontal lines in the spectra may be estimated in the limit of large wave numbers, Qa 1, giving a semi-infinite layer of Al 2 O 3 with its surface covered by either doped graphene in Fig. 9(c) or intrinsic graphene Fig. 9(d).By using the static limit, χ stat = −Q/(4v F ), for the graphene response function in both figures, we again obtain the dispersion 1 + S (ω) + π/v F = 0, giving two phonon energies of about 54 and 86 meV, which closely match the positions of the two horizontal lines in Figs.9(c) and 9(d).Noting that χ stat = −Q/(4v F ) is an exact representation of the static response function for intrinsic graphene in Fig. 9(d), but is only approximately valid for doped graphene for Q k F ≈ 0.3 nm −1 in Fig. 9(c), it is remarkable that the energies of the horizontal lines coincide in those figures.This indicates that the Dirac plasmon in doped graphene plays negligible role in screening of the substrate phonons at wave numbers greater than the Fermi wave number, where the regime of static screening by graphene prevails.
It is further interesting that the phononlike horizontal lines emerge in Figs.9(a) and 9(c) already for wave numbers Q k F ≈ 0.3 nm −1 , whereas in Fig. 9(d) they emerge right after crossing the edge ω = v F Q of the interband π ↔ π * transitions in pristine graphene.This behavior of the phononlike lines in Fig. 9(d) shows that the main role of pristine graphene layer is to provide a screening mechanism, which shifts phonon modes at large wave numbers without a significant reduction of their spectral weight.

IV. CONCLUSIONS
We have presented an ab initio calculation of the dynamic response function of single-layer graphene at the level of random phase approximation (RPA) that includes all electronic transitions within and between σ and π bands and covers a broad range of wave vectors [4].The results of such ab initio RPA approach were compared with simpler, analytical and semianalytical approaches at different levels of approximation by using the model system of a sandwich structure consisting of two graphene layers separated by a slab of Al 2 O 3 between them.Our main focus is on the range of frequencies from THz to midinfrared (MIR), where such sandwich structures represent a typical building block for devices used in nanophotonic and nanoplasmonic applications of graphene.With the Dirac plasmon in doped graphene being of the foremost interest in those applications, its hybridization with the Fuchs-Kliewer phonons in the nearby insulating layers, such as Al 2 O 3 , presents an ideal testing ground for various theoretical approaches to the dynamic response.
Using a typical thickness of 5 nm for the Al 2 O 3 layer, we have neglected the direct electron coupling between two graphene layers, and evaluated the dynamically screened Coulomb interaction of interest for probing the surface of our sandwich structure, e.g., via the reflection electron energy loss spectroscopy.Given the special interest for graphene-based devices in nanophotonic applications in the THz to MIR frequency range, we have focused on exploring the limit of long wavelengths in our theoretical approaches.Since ab initio calculations are extremely demanding for optically small wave numbers [4], we have devised a "shortcut" by using what we call optical ab initio RPA method, which only requires calculation of the conductivity of graphene with a vanishing wave number, Q = 0.Moreover, we have shown that the electron energy loss spectra obtained with the optical ab initio RPA method at small wave numbers can be seamlessly sewn together with the full ab initio RPA method at a wave number on the order of 0.1 nm −1 .Since ab initio calculations of the frequency-dependent conductivity of graphene are far less demanding than the full ab initio RPA calculations, we offer in this way a faster yet accurate extension of the ab initio RPA method in the optical limit.
Next, we have made a comparison of ab initio calculations of the spectra with those obtained by using the response function at the RPA level for graphene's π electron bands in the Dirac cone approximation, where electrons behave as massless Dirac fermions (MDF) [40,41].We have found that the "sewn-together" ab initio RPA method agrees quite well with the MDF-RPA method for wave numbers up to about 0.3 nm −1 (approximately the Fermi wave number k F when graphene layers are doped with the Fermi energy of E F = 0.2 eV).Moreover, we have explored the optical, Q = 0, limit of the MDF-RPA method as well, where simple analytical expressions are available for a frequency-dependent conductivity of graphene in the form of a sum of the intraband (i.e., Dirac) contribution and the interband contribution [42].It was found that the optical ab initio RPA and the optical MDF-RPA methods are in perfect agreement for wave numbers up to about 0.1 nm −1 .Both methods gave rather rich spectra exhibiting intricate patterns corresponding to eigenmodes due to plasmon-phonon hybridization taking place in our sandwich structure.
The exact character of those hybridized modes was further assessed by using dispersion relations based on the optical MDF-RPA expressions for graphene conductivity with zero damping, along with a model dielectric function for the bulk of Al 2 O 3 , which supports two principal transverse optical phonons.Besides the symmetric graphene-Al 2 O 3 -graphene structure, where both graphene layers were doped with E F = 0.2 eV, we have also studied in some detail hybridization taking place in asymmetric structures where one of the graphene layers is either pristine (E F = 0) or is replaced with vacuum.While the complete plasmon-phonon hybridization picture was clearly explained by the dispersion relations from the "optical MDF-RPA" method for wave numbers Q 0.1 nm −1 , we have shown that using only the Drude contribution in graphene's conductivity suffices at wave numbers Q 0.01 nm −1 .In this way, we have reaffirmed the widespread popularity of the Drude model for doped graphene in the THz applications.
On the other hand, it was also found that, while the Dirac plasmonlike modes generally subside with increasing wave numbers, Q k F = 0.3 nm −1 , due to Landau damping, i.e., upon entering the continuum of interband electron-hole transitions in graphene, there are nondispersing, phononlike modes that emerge in the same range of wave numbers.Energies of those phononlike modes could be reproduced by considering conditions for the Fuchs-Kliewer phonons in the surface of a semi-infinite Al 2 O 3 layer, where the presence of a graphene layer on that surface provides a mechanism for additional static screening of those phonon modes.
We believe that our work sheds light on the mechanism of plasmon-phonon hybridization in a sandwichlike structure with two graphene layers and, more importantly, provides a thorough assessment of a broad range of methods and approximations that can be used in calculations of the response function of graphene in the optical limit for such structures.In particular, our findings reinforce confidence in using approaches based on the frequency-dependent conductivity of graphene, where analytical expressions of the optical MDF-RPA method provide an accurate description of the dynamic response even up to wave numbers on the order of 0.1 nm −1 .

FIG. 1 .
FIG. 1.(a) Crystal structure of a gr-Al 2 O 3 -gr composite.(b) Simplified model where the sapphire layer is approximated by a homogeneous dielectric slab described by local dielectric function S (ω), whereas graphene sheets are described by 2D response functions χ 1,2 (Q,ω).

5 FIG. 3 .FIG. 4 .
FIG. 3. The EELS in single-layer doped graphene with E F = 1 eV.The black vertical line denotes the boundary wave number Q b = 0.34 nm −1 , which separates the EELS obtained using the optical ab initio RPA method (left region) and the full ab initio RPA method (right region).The thin white lines show the lower, ω = v F (Q − 2k F ), and the upper, ω = v F Q, edges of the intraband π * ↔ π * electronhole excitations, as well as the lower edge, ω = 2E F − v F Q, of the interband π ↔ π * electron-hole excitations in the Dirac cone approximation.
) and 4(b) show the low-energy EELS in singlelayer doped graphene (E F = 1 eV) where the boundary wave number (denoted by a vertical black line) Q b = 0.34 nm −1 ) and 5(b) show the EELS for double-layer graphene in vacuum, where the separation between the layers is a = 5 nm and the graphene layers are equally doped such that E F = 200 meV.In order to examine the range of validity -Im[1/ε(ω

FIG. 6 .
FIG.6.The EELS in reduced (Q,ω) space in a symmetric gr-Al 2 O 3 -gr sandwich calculated using (a) the optical ab initio RPA method and (b) the optical MDF-RPA method.The separation between graphene layers (or the Al 2 O 3 slab thickness) is a = 5 nm and both layers are doped such that E F = 200 meV.The black lines show dispersion relations of six hybridized modes (ω ± 1 , ω ± 2 and ω ± 3 ) calculated using the optical MDF-RPA method, and the red dotted lines show dispersion relations of the two highest modes (ω 0±3 ) calculated using only the Drude model.The thin white lines represent the upper edge, ω = v F Q, of the intraband π * ↔ π * electron-hole excitations in Dirac cone approximation.

FIG. 7 .
FIG. 7. The EELS in a symmetric gr-Al 2 O 3 -gr sandwich calculated using (a) the MDF-RPA method and (b) a combination of the optical ab initio RPA and the full ab initio RPA methods with the boundary wave number Q b = 0.049 nm −1 .The separation between graphene layers (or the Al 2 O 3 slab thickness) is a = 5 nm and they are both doped with equal densities, such that E F = 200 meV.The magenta dotted lines show dispersion relations of three odd modes (ω − 1 , ω − 2 and ω − 3 ), while the white dotted lines show dispersion relations of three even modes (ω + 1 , ω +2 and ω + 3 ) calculated using the optical MDF-RPA method.The thin white lines represent the electron-hole excitations edges in gr(200 meV), as described in Fig.3.

Figure 8 (
b) shows the opposite situation where the gr(200 meV) layer is placed in the z = −a/2 plane, while the external probing electron is located right next to a free Al 2 O 3 surface at z = a/2.This causes a reduction in the plasmon intensity, but the intensity of the uppermost FK phonon ω 4 becomes greatly enhanced (see the discussion of Fig. 9 below).It can be seen that the intensity patterns the EELS in Figs.8(a) and 8(b) are in excellent agreement with the dispersion relations for ω 1 , ω 3 , ω 4 , and ω 5 .The next step is to compare Figs.8(a) and 8(b) with Figs.8(c) and 8(d) in order to examine how spectra change when the layer of vacuum is replaced by a pristine graphene layer, gr(0 eV).It may be clearly seen in Figs.8(c) and 8(d) that the presence of gr(0 eV) reduces the intensities of plasmon and phonon modes in comparison with those shown in Figs.8(a) and 8(b)