Dielectric screening and plasmon resonances in bilayer graphene

M. Pisarra,1,2,* A. Sindona,1,3,† M. Gravina,1,3 V. M. Silkin,4,5,6,‡ and J. M. Pitarke7,8,§ 1Dipartimento di Fisica, Università della Calabria, Via P. Bucci, Cubo 30C, I-87036 Rende (CS), Italy 2Departamento de Quı́mica, Universidad Autónoma de Madrid, Calle Francisco Tomás y Valiente 7 (Módulo 13), 28049, Madrid, Spain 3INFN, sezione LNF, Gruppo collegato di Cosenza, Cubo 31C, 87036 Rende (CS), Italy 4Donostia International Physics Center (DIPC), Paseo Manuel de Lardizabal 4, E-20018 Donostia San Sebastian, Basque Country, Spain 5Departamento de Fı́sica de Materiales, Facultad de Ciencias Quı́micas, UPV/EHU, Apartado 1072, E-20080 Donostia San Sebastian, Basque Country, Spain 6IKERBASQUE, Basque Foundation for Science, 48011 Bilbo, Basque Country, Spain 7CIC nanoGUNE, Tolosa Hiribidea 76, E-20018 Donostia San Sebastian, Basque Country, Spain 8Materia Kondentsatuaren Fisika Saila, DIPC, and Centro Fisica Materiales CSIC-UPV/EHU, 644 Posta Kutxatila, E-48080 Bilbo, Basque Country, Spain (Received 3 September 2015; revised manuscript received 23 December 2015; published 22 January 2016)


I. INTRODUCTION
Research on collective electronic excitations (plasmons) in graphene [1,2] has received considerable interest due to the very promising technological applications of this nanomaterial and its potential to be integrated in the design of a new generation of nanodevices.Indeed, graphene plasmonics has burst on the scene, as graphene layers can be cut and/or assembled with other materials to form more complex nanostructures or heterostructures [3][4][5].These novel assemblies can be suitably designed to match the desired properties, which attracts an even wider audience [6].Apart from the outstanding electronic properties of graphene itself, the plasmon properties of graphene interfaces represent a promising field of investigation; it comes as no surprise that the first graphene-based heterostructured plasmonic device is already a reality [7].In this context, bilayer graphene (BLG) represents the simplest possible system where graphene sheets are brought together to form a more complex nanostructure whose properties show remarkable similarities and differences as compared to monolayer graphene (MLG).The great interest in BLG is due to the fact that besides displaying the electronic properties of a two-dimensional (2D) material [8] it also offers peculiar features that appear as the result of the spatial separation of the two MLG planes.It is therefore important to reach a thorough understanding of the dielectric screening properties of this system.
BLG can be considered to be the simplest possible prototype for interacting 2D systems, which are known to display two plasmonic resonances [27,28], namely, the acoustic and the optical modes, as long as the two 2D systems are close enough for the modes confined on each system to significantly influence each other.BLG can also be interpreted as a multicomponent isolated quasi-2D system, which displays two plasmon excitations as long as two bands are partially occupied [29,30].
In this work, we calculate the EL function and the plasmon dispersion of both pristine (intrinsic) and doped (extrinsic) BLG in the framework of ab initio time-dependent densityfunctional theory (TDDFT) at the level of the random-phase approximation (RPA).The method relies on a self-consistent calculation of the ground-state electronic structure of pristine BLG, where linear combinations of plane waves (PWs) are used to represent the Kohn-Sham (KS) single-particle orbitals of density-functional theory (DFT).As usual, in PW codes, a fictitious 3D geometry is chosen, which in our case is made of BLG slabs replicated in the out-of-plane direction.The mutual interaction between the replicas is known to play a significant role in many-electron features like excitation energies, the frequency-dependent electronic response, photo-absorption spectra [12,13], as well as exciton bands dispersion [14].In order to get rid of these artificial effects, we use a mixed representation where the system's susceptibility is expressed as a function of the in-plane wave vector and the out-ofplane coordinate.Then, we apply a cutoff on the Coulomb potential along the out-of-plane direction [15,16] to obtain a 2D permittivity that yields an accurate description of the energy-loss function and plasmon energy dispersions.
The advantage of this approach is that it includes the fully ab initio electronic structure of BLG, which is inherently anisotropic around the Dirac point; moreover, it offers the possibility to treat BLG as an isolated 2D system without making any assumptions on the strength of the interaction between the two constituent graphene layers, the interlayer interaction being naturally embedded in the KS single-particle orbitals.While both the anisotropic dispersion and the interlayer interaction have been treated as effective parameters in previous studies [20][21][22][23][24][25], here we report a fully ab initio RPA calculation of the electronic response of BLG.In particular, we present the energy and momentum dependent loss function of undoped (intrinsic) and doped (extrinsic) BLG in an energy range from a few meV to ∼30 eV, along the inequivalent K and M Brilloin-zone (BZ) paths.First of all, we analyze the high-energy features of the EL spectrum of undoped BLG, focusing on the collective nature of the encountered EL peaks and on the anisotropy of the EL spectrum, which reflects the anisotropic band structure of the system.Then, we look at the low-energy end of the spectrum, which we find to be very sensitive to the doping level and highly anisotropic.We demonstrate that doped BLG can be correctly visualized as a multicomponent quasi-2D system with up to four inequivalent charge carriers giving rise to a complex system of collective modes and single-particle (SP) excitations interacting with one another.This is a consequence of the peculiar anisotropic band dispersion of BLG near the Fermi level.
The rest of the paper is organized as follows.In Sec.II, we provide the input parameters of our ground-state densityfunctional calculation yielding the KS single-particle energies and orbitals.In this section, we also apply a TDDFT scheme suitable for a 3D periodic system, and then we introduce a mixed space representation that will allow us to cut off unwanted interactions between the system replicas that are present in slab calculations.In Sec.III, we present and analyze our high-energy calculations.In Sec.IV, we investigate the low-energy end of the spectrum in the presence of doping.In Sec.V, we draw our conclusions.

A. DFT calculation
We begin with a density-functional calculation of BLG, as obtained in the local-density approximation (LDA) by using (i) the Perdew-Zunger parametrization [31] of the uniform-gas correlation energy, (ii) a PW basis set (normalized to the space volume and truncated with a cutoff energy of 25 Hartrees), and (iii) a norm-conserving pseudopotential of the Troullier-Martins type [32].The system geometry is made of periodically repeated BLG units (slabs), separated by a distance of ∼40 Bohr.The bilayer structure is obtained by stacking two MLG sheets with the alternate A-B stacking (or Bernal stacking) typical of graphite.The nearest-neighbor inplane and interplane distances are optimized, respectively, to the values 1.418 and 3.348 Å.The energy cutoff and the space geometry are associated to a finite set of 3D reciprocal lattice vectors G.The first-BZ integration for the self consistent run and geometry optimization is implemented with an unshifted 60 × 60 × 1 Monkhorst-Pack grid [33], which results in a 3600 k-point sampling of the BZ.The converged electron density is then used to calculate the KS electronic structure, i.e., the single-particle energies ε ν,k and orbitals |ν,k on a denser k-point mesh (720 × 720 × 1), which includes up to 60 bands, indexed by ν.Such numbers provide a reliable sampling for the study of the electronic excitations and the energy loss up to 30 eV.

B. Density-response function of a 3D periodic system
The noninteracting density-response function of a 3D periodic many-electron system is obtained, in reciprocal space, as follows [34,35]: where ρ kq νν (G) represent matrix elements of the form ρ kq νν (G) = νk|e −i(q+G)•r |ν k + q and η is a positive infinitesimal [36].The sums over band indices and wave vectors run within the first BZ over both occupied and unoccupied levels, whose population is established by the Fermi-Dirac distribution f νk , the factor of 2 accounting for the spin degeneracy.
The exact interacting density-response function can be obtained in the framework of TDDFT, as follows [38]: where v G 1 G 2 represent the Fourier coefficients of an effective electron-electron interaction.In the framework of linear-response theory, the inelastic scattering cross section corresponding to a process in which (after the scattering of external electrons or photons) an electronic excitation of wave vector q + G and energy ω is created at the 3D periodic system is proportional to the imaginary part of one diagonal element of the so-called inverse dielectric matrix: where v 0 GG 1 represent the Fourier coefficients of the bare Coulomb interaction: v 0 GG 1 = 4πδ GG 1 /|q + G| 2 .In the RPA, the effective electron-electron interaction is taken to coincide with its bare Coulomb counterpart, i.e., v G 1 G 2 = v 0 G 1 G 2 and, therefore, The so-called crystal local-field effects are included by taking the full matrix GG in the process of obtaining the diagonal part of the inverse dielectric matrix ( −1 ) GG .We have verified that with the use of 51 G vectors of the form {0,0,G z }, sorted in length order from smallest to largest, crystal local-field effects are properly taken into account for the momenta q and the energies ω under consideration.
Collective excitations (plasmons) are dictated by the zeros in the real part of the so-called macroscopic dielectric function (permittivity): The so-called energy-loss function is related to the imaginary part of the inverse permittivity:

C. Energy-loss function and plasmons in 2D systems
The approach outlined above is accurate for purely 3D periodic systems, repeated slabs of BLG for example.The properties of this system, however, do not necessarily concide with the properties of one single BLG, because of the longrange behavior of the Coulomb interaction.A reasonably large distance between the system replicas, such as the 40 Bohr value that we use here, ensures a negligible overlap between the KS single-particle orbitals entering χ 0 .Nonetheless, in the building up of the interacting density-response function χ (and, therefore, the energy-loss funcion E LOSS ) the electron densities that are localized at each slab may still interact with one another through the long-range Coulomb interaction.Such an unwanted effect decreases too slowly with the increase of the separation between the repeated slabs, so that a wellconverged calculation of the inverse dielectric matrix becomes computationally unaffordable even for finite momenta.In the long-wavelength limit (q → 0), the G = G = 0 component of the Coulomb interaction diverges; hence, in this limit contributions to the interacting density-response function coming from the interaction between the electron densities that are localized in the repeated slabs are not negligible, no matter how large their physical distance is.It turns out that in this long-wavelength limit, the electronic response and the collective modes of a system composed of replicated 2D slabs are artificially similar to those encountered in a truly 3D periodic system.
A possible work-around to this problem is to cut the unwanted contribution from the interaction between the replicas by using a mixed (direct-reciprocal) space representation [15,16].First of all, we separate the in-plane and outof-plane components of the reciprocal lattice vector, which we represent as G = {g,G}.We then select excitation processes where the momentum transfer q lies within the plane of the slab (i.e., q = {q x ,q y ,0}), and we look at the following Fourier transforms: where z represents the out-of-plane coordinate, which we restrict to a finite region [a,b] containing the unit cell of the 3D system and excluding, therefore, the replicas.Introducing Eq. ( 2) into Eq.( 8), one finds χgg (z,z ) = χ 0 gg (z,z ) where (see Appendix) and vg 1 g 2 (z 1 ,z 2 ) is (in the RPA) A back Fourier transform of Eq. ( 9) with respect to the out-of-plane coordinates z,z yields the Fourier coefficients χGG , which have no spurious contribution from the interaction between the BLG replicas.Such an operation restores the advantage of dealing with purely algebraic matrix operations and, at the same time, keeps unaltered the formal definitions of Eqs. ( 2)-( 6), on which we simply add a bar symbol to recall that the usual Coulomb coefficients v 0 GG need to be replaced by those of Eq. ( 10) [37].
For symmetry reasons and computational advantage, we choose a = −L z /2 and b = L z /2, in which case the matrix elements of Eq. ( 10) reduce to the simple form would bring us back to the original quantities (with no bar symbol), as defined in the previous paragraph.

III. INTRINSIC BLG
We now focus on the 2D EL function ĒLOSS of undoped (intrinsic) BLG, which we report in Fig. 1 for energies ω below 30 eV and momenta q along the high-symmetry inequivalent directions K (a) and M (b).We see that ĒLOSS is dominated by two dispersive features, characterized by a narrow and a broad peak that at long wavelengths (q ∼ 10 −2 Bohr −1 ) are located at ω ∼ 4-5 eV and ω ∼ 14-15 eV, respectively.These structures reflect coherent many-body excitations following one-electron transitions that connect the sampled k-points with high density of states (underlying the well-known Van Hove singularities) and involve the π − π * and σ − π * (π − σ * ) bands [39].The lower-energy peak lying (for q = 1.5 × 10 −2 Bohr −1 ) at 4.5 eV is the so-called π plasmon, while the higherenergy peak lying (for q = 0.9 × 10 −2 Bohr −1 ) at ∼ 14.5 eV is the σ + π plasmon.These two structures are also observed in the EL functions of both MLG [9] and graphite [40].As occurs in MLG, the π and σ + π plasmons of BLG are red-shifted with respect to their graphite counterparts and display a distinct long-wavelength energy dispersion.
In a recent study, the plasmon character of the π and σ + π peaks that are observed in MLG was analyzed carefully [19].The conclusion was that there is a continuous transition from a single-particle character to a collective character of the corresponding excitations, going from a regime (at q 0.03 Bohr −1 ) where the π and σ + π peaks simply correspond to π → π * and σ → σ * single-particle transitions (with very little screening) to a regime (at q 0.03 Bohr −1 ) where these excitations can be safely treated as π and σ + π plasmons.In Fig. 2, we plot the real and imaginary parts of the permittivity ¯ M of bilayer graphene, together with the corresponding energy-loss function ĒLOSS , as a function of ω, for two representative values of the momentum q along the K direction, namely q = 0.06 and q = 0.15 Bohr −1 .In both cases, the π and σ + π peaks in the EL spectra correspond to dips in Re(¯ M ) in an energy region where Im(¯ M ) is rather small, so we can treat them as damped plasmons.For the smallest momentum under consideration (q = 0.06 Bohr −1 ), Re(¯ M ) is indeed zero at the energy of the π peak, which represents, therefore, a well-defined plasmon; indeed, we have found that this is also the case for similar values of the momentum q along the M direction.In the case of the σ + π peak, Landau damping is larger; but one can still treat this excitation as a damped plasmon.
A deeper analysis of the curves of Fig. 1 shows that for magnitudes of the momentum that are below 0.2-0.3Bohr −1 the EL spectrum is basically independent of the direction of q, while substantial differences are observed for magnitudes of the momentum that are above 0.4-0.5 Bohr −1 .Such an anisotropic dispersion reflects the anisotropic band structure of BLG and its sixfold symmetry.A more thorough picture of the EL spectrum and its anisotropic properties is given in Fig. 3, where ĒLOSS is represented as a density plot vs q and ω along the directions K and M. In Figs.3(a) and 3(b), we also draw blue dots corresponding to the values of ω (for each q) where the plasmon condition is fulfilled, i.e., where the real part of the permittivity is zero in an energy region where the imaginary part is small.
A direct comparison of Figs.3(a) and 3(b) shows that the EL spectrum is isotropic for q < 0.2 Bohr −1 , as discussed above.The σ + π peak appears (along both directions under study) as a strongly dispersing broad peak whose width increases with q; it is not possible to determine whether this dispersion is linear or quadratic, due to the variation of the peak width and the limited energy window under study.The anisotropic effect is evident in the dispersion of the π plasmon, as the momentum increases.In particular, the π peak splits along M for q > 0.4 Bohr −1 [two peaks are clearly resolved in Fig. 3 FIG.3. EL function ĒLOSS of undoped BLG vs ω (in eV) for a momentum q (in Bohr −1 ) along the following directions: (a) K and (b) M. The blue dots represent the (q,ω) values at which Re(¯ M ) = 0 in a region where Im(¯ M ) is small.In (c), the 3D EL function E LOSS is shown, for q along K.
dispersion along both directions.These features were also observed in MLG [12,13].At longer wavelengths (q 0.1 Bohr −1 ), however, the linear dispersion is replaced by a square-root dispersion (as in MLG [18,19]), and at the longest wavelengths q 0.05 Bohr −1 (not visible in Fig. 3), we observe an overall drop of intensity in the π peak, which acquires a q 2 dispersion also predicted for MLG [19].
When we evaluate the uncorrected 3D energy-loss function E LOSS (with no bar symbol) described in Sec.II B, the dispersion of the π plasmon is found to be linear for all values of q, as reported for MLG in Refs.[12,13].The low-q linear dispersion of the π plasmon found in these calculations is simply the result of the spurious interaction between the system replicas that are present in a 3D approach.
Apart from the π and σ + π peaks, the EL spectrum of intrinsic BLG also exhibits a smooth background intensity that appears as the result of single-particle interband excitations.For q ∼ 0, the single-particle background starts at ω = 0, reflecting the semimetallic feature of BLG that is characterized, at the absolute zero, by a filled valence band and an empty conduction band touching each other at the Fermi level in correspondence of the K points; one-electron transitions are, therefore, allowed for arbitrarily small energies when q = 0. On the other hand, for q = 0 the EL spectra have a low-energy gap [white region in Figs.3(a) and 3(b)], whose width depends on q.This feature is the result of one-electron transitions not being allowed for (q,ω) within the gap, whose shape reflects the band dispersion near the Fermi level.

IV. EXTRINSIC BLG
In order to investigate the EL spectrum of extrinsic BLG, we mimic the effect of injected/ejected electrons by changing the occupation factors in Eq. ( 1) and assuming that the doping has a negligible effect on the KS electronic structure.In semimetals like MLG and BLG, even a small doping yields the formation of a quasi-2D gas of charge carriers (electrons or holes).The fact that the conduction (valence) band becomes partially occupied (unoccupied) causes the appearance of low-energy intraband SP transitions and collective excitations; the doping also causes the opening of a gap for interband transitions at low momenta.Both these effects are simply absent in pristine (intrinsic) BLG with a fully occupied valence band and a fully unoccupied conduction band.
In the following, we consider the case of positive/negative doping inducing Fermi-level shifts E F in the range of −1 to 1 eV, which corresponds to (i) a fraction of removed/added electrons per unit cell that is below 0.12 and (ii) an absolute value of the charge-carrier density that is below ∼2.3 × 10 14 cm −2 [see Table I].In order to understand how these doping levels may affect the EL spectrum, let us have a look at the π and π * band dispersions of BLG near the K point, along high-symmetry paths of the first BZ (Fig. 4).Unlike MLG, BLG presents four bands of p z symmetry, say,   extrinsic Fermi level crosses both bands, π * 1 and π * 2 (or π 1 and π 2 ), thereby yielding two Fermi surfaces per inequivalent K point, which follow the anisotropic features of the band dispersions.A rich EL spectrum is, therefore, visible when the system is probed along high-symmetry directions of the first BZ.The anisotropic Fermi surfaces lead to a peculiar situation in which inequivalent charge carriers are present, giving rise to the appearance of various collective excitations, as it happens in the general case of a multicomponent 2D electron gas.We also remark that positive and negative doping are not equivalent and must, therefore, be treated separately.This is due to the fact that the KS band dispersions are not symmetric for energies above/below the pristine Fermi level.Such differences could not be detected in previous continuum or tight-binding calculations that were based on two or four symmetric bands [20][21][22][23][24][25].
We have carried out energy-loss calculations for energies ω up to 30 eV, and we have found that at least for values of E F that are below 1 eV the π and σ + π structures are insensitive to the doping.At the low-energy end, however, the energy-loss spectrum is very sensitive to the doping, as expected.In Figs. 5  and 6, we give a density plot of the energy-loss function ĒLOSS of BLG vs q and ω, as obtained for the values of E F listed in Table I at energies ω that are below 3 eV.A case-by-case analysis is given below.

A. E F = ±1.0 eV
We begin with Fig. 5(b), which exhibits the energy-loss function of extrinsic BLG for E F = 1.0 eV and q along The plot markers (i.e., dots and diamonds) represent the (q,ω) values at which Re(¯ M ) = 0 in a region where Im(¯ M ) is zero or small.The blue dots correspond to the conventional 2D plasmon.The red diamonds correspond to an acoustic plasmon.The (q,ω) regions where SP intraband and interband transitions occur are indicated in (b).The color code is the same as in Fig. 3.
M. Some remarkable features appear here with respect to the undoped case of Fig. 3. First of all, the loss function is not zero (for q = 0) even at very low energies, which is due to the fact that doping allows for the occurrence of intraband SP excitations.Secondly, we observe that (i) there is a triangular region (at long wavelengths) with no energy loss and (ii) a very well defined energy-loss peak (that is absent in the case of intrinsic BLG) is visible outside the triangular region for q 0.05 au.This is a signature of a well-defined intraband 2D plasmon (similar to the conventional 2D plasmon of MLG), which also exists at lower momenta in a region (the triangle) where there is no damping and the energy-loss function becomes a delta function.At low momenta, we track the energy dispersion of this plasmon [blue dots in Fig. 5(b)] by looking at the zeros of the real part of the permittivity in a region where the imaginary part is either zero (undamped plasmon inside the triangle) or small (damped plasmon outside the triangle).The energy of this conventional 2D plasmon is FIG.6.The same as in Fig. 5, but now for negative doping, i.e., negative values of E F , with the black triangles corresponding to an interband 2D plasmon.
found to be (at long wavelengths) proportional to √ q.Thirdly, we observe that at lower energies there is a second sequence of zeros of the real part of the permittivity in a region where the imaginary part is either zero or small [red diamonds in Fig. 5(b)].The energy of this new acoustic plasmon (which is absent in MLG) is proportional to q.
The acoustic plasmon that is visible in Fig. 5(b) is related to the fact that for E F = 1.0 eV the extrinsic Fermi level crosses two bands, π * 1 and π * 2 , giving rise to two Fermi surfaces, as is the case of a 2D electron gas with two components [29].Hence, two different types of charge carriers are present, leading to the appearance of two collective modes [both visible in Fig. 5(b)]: one in which the two types of carriers oscillate in phase (the optical mode that in our case is the conventional 2D plasmon represented by blue dots) and another one where the two types of carriers oscillate out of phase (the acoustic mode represented by red diamonds).The presence of an acoustic plasmon in doped BLG was predicted by Das Sarma et al. [24]; they found, however, an overdamped plasmon.Our ab initio calculation shows that this acoustic plasmon happens to be (at long wavelengths) a very well defined undamped collective excitation along the M direction.
Along the K direction [Fig.5(a)], we still observe (i) a triangular (q,ω) region with no energy loss and (ii) the conventional 2D plasmon whose energy is (at long wavelengths) proportional to √ q; but now instead of the well-defined acoustic plasmon of Fig. 5(b) we see two new distinct damped dispersive energy-loss peaks at ω 1 eV, q > 0.06 Bohr −1 .These resonances occur, however, in a region where intraband SP transitions occur and do not correspond to zeros in the real part of the permittivity.In order to understand this situation, we look again at the band structure of Fig. 4. A momentum transfer q along the K direction can excite the electron gas along two inequivalent branches of the band structure at the inequivalent K points, namely, Kand K-M [see Fig. 4(a)].As a result, at E F = 1.0 eV, four different types of carriers are expected to exist along K in the presence of two anisotropic Fermi surfaces.However, along the Kbranch the slopes of the π * 1 and π * 2 bands (and, therefore, the corresponding Fermi velocities) are nearly indistinguishable.The combination of this single type of carrier with two different types of carriers along the K-M branch yields the appearance of the two distinct damped acoustic collective excitations that are visible in Fig. 5(a).Now we consider the case of negative doping ( E F = −1.0eV).First of all, we note that now one-electron transitions towards all the π and π * bands are allowed already for E F < −0.3, which implies that (i) the triangular no-loss region that was visible for a positive doping [Figs.5(a .Secondly, we note that now along both the Kbranch and the K-M branch (not only along the Kbranch), the slopes of the π 1 and π 2 bands are nearly indistinguishable, which implies the appearance of one single damped acoustic collective excitation for momenta q along the K direction [Fig.6(a)].Finally, we see that the conventional 2D plasmon remains insensitive to the sign of the doping.

B. E F = ±0.5 eV
In this case, a richer spectrum of acoustic collective excitations would be expected to occur, as the π * 1 and π * 2 bands now cross the extrinsic Fermi level with different slopes along both branches of the K and M directions.Nevertheless, for this shift of the Fermi level one-electron interband transitions occur in a wider (q,ω) region, which is reflected in the fact that all acoustic collective excitations happen to be strongly damped, both for positive and negative doping.The conventional 2D plasmon remains well defined, although considerably more damped than in the case of a higher doping.

C. E F = ±0.25 eV
This case reflects an interesting situation, since the π * 2 (π 2 ) band remains completely unoccupied (occupied); this means that only one Fermi surface is present.Furthermore, for singleparticle energies below 0.3 eV, the band structure is mainly isotropic (see Fig. 4), which means that (i) there is no room for acoustic excitations (as one single type of carrier is present) and (ii) the energy loss is insensitive to the direction of the momentum transfer q [see Figs.5(e), 5(f), 6(e), and 6(f)].The energy loss is also insensitive to whether the doping is positive or negative.
Figures 5(e), 5(f), 6(e), and 6(f) all show the existence of two dispersive features characterized by narrow peaks.The higher-energy peak (located at ω ∼ 0.35 eV at q = 0) is a signature of an interband plasmon, that is a collective excitation that follows one-electron transitions connecting points with high density of states in the π * 1 and π * 2 (π 1 and π 2 ) bands.This feature has the same nature as the so-called π plasmon, although zeros in the real part of the permittivity [black triangles in Figs.6(e) and 6(f)] are now found only in the case of negative doping.The lower-energy peak represents the conventional 2D plasmon that is considerably distorted (at the longest wavelengths) by the presence of the higherenergy peak and has a dispersive energy that is not zero in the long-wavelength limit.

V. CONCLUSIONS
We have presented a fully ab initio RPA analysis of the plasmon structure of intrinsic and extrinsic BLG along the inequivalent K and M directions.The quasi-2D nature of BLG has been treated carefully by combining a scheme of periodically repeated slabs of BLG with a 2D correction that completely eliminates the artificial long-range interaction between the replicas.
We have carried out energy-loss calculations for energies ω up to 30 eV, and we have found that at least for values of E F that are below 1 eV the π and σ + π structures are insensitive to the doping.The quasi-2D nature of BLG is reflected in the 2D nature of the π mode at long wavelengths, which we have found to have (at q 0.05 0.1 Bohr −1 ) an energy that is proportional to √ q, as shown by recent experiments on MLG.
At low energies, energy-loss spectra are found to be very sensitive to doping and even to whether doping is positive (injected electrons) or negative (injected holes).In addition to a conventional 2D plasmon whose energy is (at long wavelengths) proportional to √ q, we have found a myriad of acoustic collective excitations due to the presence of two distinct Fermi surfaces and various types of charge carriers.In particular, a well-defined undoped acoustic plasmon is predicted to exist for momenta along the M direction and a positive shift in the Fermi level of 1 eV, and various damped acoustic plasmons (of a different origin) are found to exist for momenta along the K direction and shifts in the Fermi level that are (in absolute value) above ∼0.3eV.A new low-energy interband plasmon (similar to the so-called π plasmon) has also been found to exist in the case of low doping ( E F = ±0.25 eV), as a consequence of the π * 2 (π 2 ) band being completely unoccupied (occupied).
Our results call for a more thorough experimental investigation of the dielectric screening properties and many body excitations of the BLG system.In particular, the presence of multiple distinct plasmon modes in exstrinsic BLG, if confirmed by experiments, could have a strong impact influencing the future design of BLG-based plasmonic nanodevices.

APPENDIX : INTERACTION KERNEL IN THE DIRECT-RECIPROCAL SPACE REPRESENTATION
In a mixed (direct-reciprocal) space representation where the system's susceptibility is expressed as a function of the in-plane wave vector and the out-of-plane coordinate, Eq. ( 2) becomes χgg (z,z ) = χ 0 gg (z,z ) + ( χ 0 v χ) gg (z,z ), (A1) where the second term on the RHS has the explicit form ( χ 0 v χ ) gg (z,z ) = Performing this double space integration, one finds f G 1 G 2 = e i(aG 1 −bG 2 ) e −q(b−a) (q + iG 1 )(q + iG 2 ) + e i(bG 1 −aG 2 ) e −q(b−a) (q − iG 1 )(q − iG 2 ) The choice a = −L z /2 and b = L z /2 simplifies the implementation and speeds up the calculation, but is not mandatory.A different choice of the integration region, being large enough to include the z extension of the unit-cell electron density and small enough to exclude the replicas, should yield equivalent results.
Finally, we note that the interaction kernel of Eq. ( 10) is simply )

E
FIG.1.EL function ĒLOSS vs ω (in eV) for a momentum q (in Bohr −1 ) along the following directions: (a) K and (b) M. Curves corresponding to different magnitudes (q) of the momentum are displaced vertically for clarity of the presentation.

4 ω
FIG.Real and imaginary parts of the macroscopic dielectric function: Re(¯ M ) and Im(¯ M ), and the EL function ĒLOSS , vs ω 30 eV, for q = 0.06,0.15Bohr −1 along the K direction.
(b)], while no splitting is exhibited along K.In the intermediate momentum region, 0.1-0.4Bohr −1 , the π plasmon has a linear

π 1 ,
π 2 , π * 1 , and π * 2 [Figs.4(a) and 4(b)].The dispersion of these bands is parabolic near the K point, and their occupation depends on the doping values.In pristine BLG, the π 1 and π 2 bands are fully occupied, while the π * 1 and π * 2 bands are fully unoccupied, offering markedly anisotropic dispersions already for in-plane wave-vectors larger than ∼0.1 Bohr −1 .All these features complicate the interpretation of the low-energy loss spectrum of BLG in terms of SP and collective excitations.In particular, for | E F | 0.3 eV the

FIG. 4 .
FIG. 4. (a) and (b) Anisotropic band structure of BLG at the K point along the inequivalent k paths shown in (c), with the Fermi level of pristine BLG set to zero energy.The horizontal dashed lines represent the Fermi-level shifts listed in Table I.(c) First BZ of BLG, with the inequivalent directions used for the band plots of (a) and (b).

FIG. 5 .
FIG. 5. EL function ĒLOSS of positively doped BLG vs ω (in eV)for a momentum q (in Bohr −1 ) along the following directions: [(a), (c), and (e)] K and [(b), (d), and (f)] M. The plot markers (i.e., dots and diamonds) represent the (q,ω) values at which Re(¯ M ) = 0 in a region where Im(¯ M ) is zero or small.The blue dots correspond to the conventional 2D plasmon.The red diamonds correspond to an acoustic plasmon.The (q,ω) regions where SP intraband and interband transitions occur are indicated in (b).The color code is the same as in Fig.3.
) and 5(b)] is now significantly reduced [Figs.6(a) and 6(b)] and (ii) the well-defined acoustic plasmon that was visible along the M direction [Fig.5(b)] is now nearly completely damped [Fig.6(b)] ACKNOWLEDGMENTSM.P. and M.G. acknowledge financial support by the European Commission, the European Social Fund, and the Regione Calabria, (POR) Calabria -FSE 2007/2013, and the hospitality of CIC nanoGUNE and the Donostia International Physics Center (DIPC).V.M.S. acknowledges the partial support from the Basque Departamento de Educación, UPV/EHU (Grant No. IT-756-13) and the Spanish Ministry of Economy and Competitiveness MINECO (Grant No. FIS2013-48286-C2-1-P).