# Excitation spectrum and high energy plasmons in single- and multi-layer graphene

###### Abstract

In this paper we study the excitation spectrum of single- and multi-layer graphene beyond the Dirac cone approximation. The dynamical polarizability of graphene is computed using a full -band tight-binding model, considering the possibility of inter-layer hopping in the calculation. The effect of electron-electron interaction is considered within the random phase approximation. We further discuss the effect of disorder in the spectrum, which leads to a smearing of the absorption peaks. Our results show a redshift of the -plasmon dispersion of single-layer graphene with respect to graphite, in agreement with experimental results. The inclusion of inter-layer hopping in the kinetic Hamiltonian of multi-layer graphene is found to be very important to properly capture the low energy region of the excitation spectrum.

###### pacs:

73.21.Ac,81.05.ue,79.20.Uv## I Introduction

One of the main issues in the understanding of the physics of graphene is the role played by electron-electron interaction.For a recent review on electron-electron interaction in graphene, see et al. (2010) Several collective modes as low and high energy plasmons, as well as plasmarons, are a consequence of electronic correlations and have been measured in this material. The high energy -plasmons have been observed in electron energy-loss spectroscopy (EELS),Kramberger et al. (2008); Eberlein et al. (2008); Gass et al. (2008) inelastic X-ray scattering (IXS)Reed et al. (2010) or optical conductivity.Mak et al. (2011) Recently, a plasmaron mode (which is a result of coupling between electrons and plasmons) has been measured in angle-resoved photoemission spectroscopy (ARPES).Bostwick et al. (2010)

At low energies, long range Coulomb interaction leads, in doped graphene, to a gapless plasmon mode which disperses as ,Liu et al. (2008) and which can be described theoretically within the random phase approximation (RPA).Shung (1986); Wunsch et al. (2006); Hwang and Sarma (2007); Polini et al. (2008); Pyatkovskiy (2009); Gamayun (2011) The low energy linear dispersion relation of graphene is at the origin of a new series of collective modes predicted for this material and which do not exist for other two-dimensional electron gases (2DEG), as inter-valley plasmonsTudorovskiy and Mikhailov (2010) or linear magneto-plasmons,Roldán et al. (2009) which can be described within the RPA as well. For undoped graphene, the inclusion of ladder diagrams in the polarization can lead to a new class of collective modesGangadharaiah et al. (2008) as well as to an excitonic instability.Khveshchenko (2001); Aleiner et al. (2007); Gamayun et al. (2009); Wang et al. (2010); González (2010)

However, much less is known about the high energy -plasmon which, in the long wavelength limit, has an energy of the order of 5-6 eV, and which is due to the presence of Van Hove singularities in the band dispersion. For single-layer graphene (SLG), this mode has been studied by Stauber et al.Stauber et al. (2010) and by Hill et al.Hill et al. (2009) in the RPA. Yang et al. have included excitonic effects and found a redshift of the absorption peak,Yang et al. (2009) leading to a better agreement with the experimental results.Eberlein et al. (2008) Here we extend those previous works and study the excitation spectrum of SLG and multi-layer graphene (MLG) from a tight-binding model on a honeycomb lattice. By means of the Kubo formula, the non-interacting polarization function is obtained from the numerical solution of the time-dependent Schrödinger equation. Coulomb interactions are considered in the RPA, the validity of which is discussed. We also consider the effect of disorder in the system, which lead to a considerable smearing of the Van Hove singularities in the spectrum. Our results show a redshift of the -plasmon mode in graphene with respect to graphite, as it has been observed in the experiments.Kramberger et al. (2008); Eberlein et al. (2008) Furthermore, the inclusion of inter-layer hopping is found to be very important to capture the low energy region of the spectrum in MLG.

The paper is organized as follows. In Sec. II we describe the method that we use to compute the dynamical polarization function of SLG and MLG. In Sec. III we give results for the excitation spectrum of SLG, considering the effect of disorder and electron-electron interaction. The spectrum of MLG is described in Sec. IV. In Sec. V we compare our results to recent experimental data. Our main conclusions are summarized in Sec. VI.

## Ii Description of the method

The tight-binding Hamiltonian of a MLG is given by

(1) |

where is the Hamiltonian of the ’th layer of graphene,

(2) |

where () creates (annihilates) an electron on sublattice A (B) of the ’th layer, and is the nearest neighbor hopping parameter. The second term of accounts for the effect of an on-site potential , where is the occupation number operator. In the second term of the Hamiltonian Eq. (1), describes the hopping of electrons between layers and . In nature there are two known forms of stable stacking sequence in bulk graphite, namely ABA (Bernal) and ABC (rhombohedral) stacking, and they are schematically shown in Fig. 1. For a MLG with an ABA stacking, is given by

(3) |

where the inter-layer hopping terms and are shown in Fig. 1. Thus, all the even layers () are rotated with respect to the odd layers () by . The difference between ABA and ABC stacking is that, the third layer(s) is rotated with respect to the second layer by (then it will be exactly under the first layer) in ABA stacking, but by in ABC stacking.Koshino (2010) In this paper we use the hopping amplitudes eV, eV and eV.Castro-Neto et al. (2009) The spin degree of freedom contributes only through a degeneracy factor and is omitted for simplicity in Eq. (1). In our numerical calculations, we use periodic boundary conditions in the plane () of graphene layers, and open boundary conditions in the stacking direction ().

The dynamical polarization can be obtained from the Kubo formula Kubo (1957) as

(4) |

where denotes the volume (or area in 2D) of the unit cell, is the density operator given by

(5) |

and the average is taken over the canonical ensemble. For the case of the single-particle Hamiltonian, Eq. (4) can be written asYuan et al. (2010a)

where is the Fermi-Dirac distribution operator, where is the temperature and is the Boltzmann constant, and is the chemical potential. In the numerical simulations, we use units such that , and the average in Eq. (II) is performed over a random phase superposition of all the basis states in the real space, i.e.,Hams and De Raedt (2000); Yuan et al. (2010a)

(7) |

where are random complex numbers normalized as . By introducing the time evolution of two wave functions

(8) | |||||

(9) |

we get the real and imaginary part of the dynamical polarization as

The Fermi-Dirac distribution operator and the time evolution operator can be obtained by the standard Chebyshev polynomial decomposition.Yuan et al. (2010a)

For the case of SLG, we will further compare the polarization function obtained from the Kubo formula Eq. (4), to the one obtained from the usual Lindhard function.Giuliani and Vignale (2005) Notice that this method can be used to calculate the optical conductivity of graphene beyond the Dirac cone approximation. Stauber et al. (2008); Yuan et al. (2010a) For pristine graphene, the dynamical polarization obtained from the Lindhard function using the full -band tight-binding model isShung (1986); Hill et al. (2009); Stauber et al. (2010)

where the integral is over the Brillouin zone, is the spin degeneracy, is the energy dispersion with respect to the chemical potential, where

(12) |

Å being the in-plane carbon-carbon distance, and the overlap between the wave-functions of the electron and the hole is given by

(13) |

In the RPA, the response function of SLG due to electron-electron interactions can be calculated as

(14) |

where is the Fourier component of the Coulomb interaction in two dimensions, in terms of the background dielectric constant , and

(15) |

is the dielectric function of the system. We will be interested on the collective modes of the system, which are defined from the zeroes of the dielectric function []. The dispersion relation of the collective modes is defined from

(16) |

which leads to poles in the response function (14). The damping of the mode is proportional to , and it is given by

(17) |

For MLG, the response function is calculated as (we use )Shung (1986)

(18) |

where Å is the inter-layer separation. Because we use open boundary conditions in the stacking direction, we define the form factor as

(19) |

The expression (18) assumes that the polarization of each layer is the same, and it is exact in two different limits: bilayer graphene and graphite. Notice that a similar effective form factor has been used to study the loss function of multiwall carbon nanotubes.Shyu and Lin (2000) Eq. (19) coincides with the commonly used form factor for a multi-layer system with an infinite number of layers:Das Sarma and Quinn (1982)

(20) |

where in this last case the periodicity ensures that is independent of layer index , with the asymptotic behavior .Das Sarma and Quinn (1982)

A crucial issue is the value of the dielectric constant for each of the cases considered, because it encodes the screening due to high energy () bands which are not explicitly considered in our calculation. A good estimation for it can be obtained from the expressionWehling et al. (2011)

(21) |

where is the dielectric constant of graphite, is the total height of the multi-layer system in terms of the number of layers and the height of a monolayer graphene Å. As expected, Eq. (21) gives for SLG at and for graphite.

We notice that the accuracy of the numerical results for the polarization function Eq. (II) is mainly determined by three factors: the time interval of the propagation, the total number of time steps, and the size of the sample. The maximum time interval of the propagation in the time evolution operator is determined by the Nyquist sampling theorem. This implies that employing a sampling interval , where are the eigenenergies, is sufficient to cover the full range of energy eigenvalues. On the other hand, the accuracy of the energy eigenvalues is determined by the total number of the propagation time steps () that is the number of the data items used in the fast Fourier transform (FFT). Eigenvalues that differ less than cannot be identified properly. However, since is proportional to we only have to double the length of the calculation to increase the accuracy by the same factor. The statistic error of our numerical method is inversely proportional to the dimension of the Hilbert space,Hams and De Raedt (2000) and in our case (the single particle representation), it is the number of sites in the sample. A sample with more sites in the real space will have more random coefficients () in the initial state , providing a better statistical representation of the superposition of all energy eigenstates.Yuan et al. (2010a)

Similar algorithm has been successfully used in the numerical calculation of the electronic structure and transport properties of single- and multi-layer graphene, such as the density of states (DOS), or dc and ac conductivities.Wehling et al. (2010); Yuan et al. (2010a, b) The main advantage of our algorithm is that different kinds of disorders and boundary conditions can be easily introduced in the Hamiltonian, and the computer memory and CPU time is linearly proportional to the size of the sample, which allows us to do the calculations on a sample containing tens of million sites.

## Iii Excitation spectrum of single-layer graphene

The particle-hole excitation spectrum is the region of the energy-momentum space which is available for particle-hole excitations. For non-interacting electrons, it is defined as the region where , as given by Eq. (4) or (II), is non-zero.Giuliani and Vignale (2005) The linear low energy dispersion relation of graphene as well as the possibility for inter-band transitions lead to a rather peculiar excitation spectrum for SLG as compared to the one of a two-dimensional electron gas (2DEG) with a parabolic band dispersion.Roldán et al. (2010) Here we focus on undoped graphene (), for which only inter-band transitions are allowed. In Fig. 3 we plot for different wave-vectors at (which is the temperature that we will use from here on in our results). The first thing one observes is the good agreement between the results obtained from the Kubo formula Eq. (4), as compared to Lindhard function Eq. (II), what proofs the validity of our numerical method. Furthermore, for the small wave-vector used in Fig. 3(a)-(b), the results are well described by the Dirac cone approximation,Wunsch et al. (2006); Hwang and Sarma (2007) but only at low energies, around , where is the Fermi velocity near the Dirac points. In particular, the continuum approximation cannot capture the peaks of around . These peaks are related to particle-hole excitations between states of the valence band with energy and states of the conduction band with energy , which contribute to the polarization with a strong spectral weight due to the enhanced density of states at the Van Hove singularities of the -bands (see Fig. 4).

Second, for larger wave-vectors [Figs. 3(c)-(f)] one observes strong differences in the spectrum depending on the orientation of , effect which has been discussed previously.Stauber et al. (2010); Stauber and Gómez-Santos (2010) If is along the -K direction, there is a splitting of the peak associated to the Van Hove singularity at . At low energies, we also observe a finite contribution to the spectral weight to the left of the peak for momenta along the -M direction [plots Figs. 3(d) and (f)]. Finally, trigonal warping effects are important as we increase the magnitude of , due to the deviation of the band dispersion with respect to the linear cone approximation. As a consequence, the constant energy contours are not any more circles around the Dirac points, but present a triangular shape. The consideration of this effect leads to a redshift of the peak with respect to the Dirac cone approximation, as seen clearly in Fig. 3(e).

Once we have discussed the clean case, we consider the effect of disorder on the excitation spectrum as explained in Sec. II. Two different kinds of disorder are considered: random local change of on-site potentials and random renormalization of the hopping, which correspond to the diagonal and off-digonal disorders in the single-layer Hamiltonian Eq. (2), respectively. The former acts as a chemical potential shift for the Dirac fermions, i.e., shifts locally the Dirac point, and the later rises from the changes of distance or angles between the orbitals. In Fig. 4 we show the DOS of SLG for different kinds and magnitudes of disorder. The DOS for clean graphene has been plotted by using the analytical expression given in Ref. Hobson and Nierenberg, 1953. The DOS of the disordered systems are calculated by Fourier transform of the time-dependent correlation functions Yuan et al. (2010a)

(22) |

with the same initial state defined in Eq. (7). As shown in Ref. Yuan et al., 2010a, the result calculated from a SLG with lattice sites matches very well with the analytical expression, and here we use the same sample size in the disordered systems. We consider that the on-site potential is random and uniformly distributed (independently on each site ) between and . Similarly, the in-plane nearest-neighbor hopping is random and uniformly distributed (independently on sites ) between and . The main effect is a smearing of the Van Hove singularity at , as observed in the right hand side inset of Fig. 4.

The effect of disorder is also appreciable in the non-interacting excitation spectrum of the system, as shown by Fig. 5. A broadening of the and peaks is observed in all the cases. Furthermore, disorder leads to a slight but appreciable redshift of the peaks with respect to the clean limit. This effect is more important for higher wave-vectors, as it can be seen in Fig. 5(c)-(d). Finally, the disorder broadening of the peaks leads in all the cases to a transfer of spectral weight to low energies (below ), as it is appreciable in Fig. 5(a)-(d).

The next step is to consider both, disorder and electron-electron interaction in the system. In the RPA, the response function is calculated as in Eq. (14). The results are shown in Fig. 6, where is plotted for the same wave-vectors and disorder used in Fig. 5. We observe that the Dirac cone approximation (black line) captures well the low energy region of the spectrum. However, the large peak at cannot be captured by the continuum approximation. They are due to a plasmon mode associated to transitions between electrons in the saddle points of the -bands. Strictly speaking, those modes cannot be considered as fully coherent collective modes, as for example, the low energy -plasmon which is present in doped graphene.Shung (1986) For doped graphene, the acoustic -plasmon is undamped above the threshold until it enters the inter-band particle-hole continuum, when it starts to be damped and decays into electron-hole pairs. However, the -plasmon, although it corresponds to a zero of the dielectric function as it can be seen in Fig. 7, it is a mode which lies inside the continuum of particle-hole excitations: at the -plasmon energy , and the mode will be damped even at . In any case, it is a well defined mode which has been measured experimentally for SLG and MLG.Kramberger et al. (2008); Eberlein et al. (2008); Gass et al. (2008); Reed et al. (2010); Mak et al. (2011) Coming back to our results, notice that the height of the peaks is reduced when the effect of disorder is considered, although the position is unaffected by it. For small wave-vectors, this mode is highly damped due to the strong spectral weight of the particle-hole excitation spectrum at this energy, as seen e.g. by the peak of at in Fig. 3(a)-(b). The position of the collective modes can be alternatively seen by the zeroes of the dielectric function Eq. (16), which is shown in Fig. 7. Notice that the Dirac cone approximation (solid black lines in Fig. 7) is completely insufficient to capture this high energy -plasmon, which predicts always a finite . As for the polarization, we see that disorder lead to an important smearing of the singularities of the dielectric function, as seen in Fig. 7. Finally, we mention that the application of our method to even higher wave-vectors and energies as the ones considered in the present work, should be accompanied by the inclusion of local field effects (LFE) in the dielectric function, which are related to the periodicity of the crystalline lattice.Adler (1962) In fact, for SLG and for wave-vectors along the zone boundary between the M and the K points of the Brillouin zone (see Fig. 2), the inclusion of LFEs leads to a new optical plasmon mode at an energy of 20-25 eV.Pellegrino et al. (2010)

## Iv Excitation spectrum for multi-layer graphene

In the following, we study the excitation spectrum and collective modes of MLG. For this, we consider not only the Coulomb interaction between electrons on different layers, but also the possibility for the carriers to tunnel between neighboring layers, as described in Sec. II. The importance of considering inter-layer hopping has been already shown in the study of screening properties of MLG.Guinea (2007) First, we see that the results are sensitive to the relative orientation between layers. In Fig. 8 we show the density of states for ABA- and ABC- stacked MLG (see Fig. 1 for details on the difference between those two orientations). As seen in Fig. 8(c)-(d), all the MLGs present a finite DOS at , contrary to SLG which has a vanishing DOS at the Dirac point. The main difference between the two kinds of stacking is that for ABC there is a central peak together with a series of satellite peaks around [Fig. 8(d)], whereas for ABA the DOS follows closer the behavior of the SLG [Fig.8(c)]. The different structure in the DOS can be understood by looking at Fig. 9, where we show the low energy band structure of a trilayer graphene with ABA [Fig. 9(a)] and ABC [Fig. 9(b)] orientations. The different jumps and peaks in the DOS of Fig. 8(c)-(d) are associated to the regions of the band dispersion marked by the horizontal red lines of Fig. 9, the energy of which depends on the values of the tight-binding parameters associated to inter-layer tunneling ( and in our case). In the two cases, we observe a splitting of the Van Hove peak, as seen in Fig. 8(e)-(f). Notice that when we have a high number of layers (e.g. above 10 layers), there is a weak effect on adding a new graphene sheet to the system, as it can be seen from the similar DOS between the 10- and the 50-layers cases of Fig. 8.

In Fig. 10 we show the non-interacting (left panels) and the RPA (right panels) polarization function of MLG, for systems made of 3, 5 and 20 layers, and for ABA- and ABC-stacking. For the spectrum in the absence of electron-electron interaction, as shown in Fig. 10(a), (c) and (e), one does not observe any specific difference in the two spectra apart from different intensities depending on the kind of stacking and on the number of layers considered for the calculation. On the other hand, the energy of the -plasmon of ABC samples is redshifted with respect to the ABA-stacking. This can be see from the relative position of the peaks of in Fig. 10(b), (d) and (e). Also, notice that the separation between the two peaks grows with the number of layers, and for a 20-layers system, the difference can be of the order of 1 eV, as it can be seen in Fig. 10(f). In the following and unless we say the opposite, all the results will be calculated for the more commonly found ABA-stacking.

For a more clear understanding about the evolution of the particle-hole excitation spectrum with the number of layers, we plot in Fig. 11(a) the imaginary part of for SLG and MLG of several number of layers, and compare the results to the polarization obtained using the Dirac cone approximation. It is very important to notice that multi-layer graphene presents some spectral weight at low energies as compared to graphene, which can be seen from the finite contribution of that appears to the left of the big peak of the graphene polarizability at (eV for the used parameters), in terms of the Fermi velocity near the Dirac point, . This is due to the low energy parabolic-like dispersion of bilayer and multilayer graphene, as compared to the linear dispersion of single layer graphene, and it can only be captured by considering the inter-layer hopping contribution to the kinetic Hamiltonian Eq. (3). Furthermore, the spectrum presents a series of peaks for , the number of which depends on the number of layers. This is due to the fact that as we increase the number of coupled graphene planes, the number of bands available for particle-hole excitations also grows leading to peaks at different energies for a given wave-vector.Zhang et al. (2008)

The difference between SLG and MLG is also relevant in the low energy region of the dielectric function, as it can be seen in Fig. 11(b). In fact, the limit of calculated within the RPA grows with the number of layers. Moreover, as we have discussed above, the zeroes of signal the position of collective excitations in the system (plasmons). In Fig. 11(b) we see that , for the small wave-vector used, crosses 0 for MLG, revealing the existence of a solution of Eq. (16), but not so for SLG, as it was pointed out in Ref. Stauber et al., 2010. However, we emphasize that the very existence of solutions for the equation for MLG does not imply the existence of long-lived plasmon modes. In fact, as we have already discussed in Sec. III, these modes disperse within the continuum of particle-hole excitations [, where is the solution of Eq. (16)], so they will be Landau damped and will decay into electron-hole pairs with a damping given by Eq. (17). Furthermore, we remember that for a given wave-vector, the energy of the mode is controlled by the background dielectric constant , as given by Eq. (21). For the systems under consideration, changes between 1 (for SLG) and 2.4 (for graphite). The value of , together with the form factor Eq. (19) that takes into account inter-layer Coulomb interaction, fix the position of the modes in each case.

The effect of disorder in MLG is considered in Fig. 12, where we show the polarization function of a 20-layer graphene system for different kinds of disorder. As in the SLG, we find that disorder leads to a slight redshift of the peaks of the non-interacting spectrum [Fig. 12(a) and (b)], together with a smearing of the peaks at and . On the other hand, the interacting polarization function presents a reduction of the intensity of the plasmon peak due to disorder, as seen in Fig. 12(c) and (d), also in analogy with the SLG case.

## V Discussion and comparison to experimental results

In this section we compute quantities which are directly comparable to recent experimental results on SLG and MLG. We start by calculating the loss function , which is proportional to the spectrum measured by EELS. Our results, shown in Fig. 13, are in good agreement with the experimental data of Ref. Eberlein et al., 2008: as in the experiments, we observe a redshift of the plasmon peak as one decrease the number of layers, as well as an increase of the intensity with the number of layers. Notice that, due to finite size effects, there is an infra-red cutoff for the wave-vectors used in our calculations which prevents to reach the long wavelength limit. In Fig. 13 we show the results for the smallest wave-vector available, and we emphasize that the peaks will be further shifted to the left for smaller values of . A further redshift of the peaks would be obtained beyond RPA, as it has been reported for single- and bilayer graphene, where excitonic effects have been included.Yang et al. (2009)

We have also used our method to study the IXS experiments of Reed et al.Reed et al. (2010) In Fig. 14(a)-(b) we plot the imaginary part of the non-interacting polarization function for SLG and MLG, for two values of similar to the ones used in Ref. Reed et al., 2010. As we have discussed in Sec. IV, inter-layer hopping leads to a finite contribution to the spectral weight in the low energy region of MLG as compared to the SLG spectrum. Notice that the number of peaks at this energy scales with the number of accessible bands and therefore, with the number of layers. We emphasize that this effect is not included by the usually employed approximation of considering MLG as a series of single-layers of graphene, only coupled via direct Coulomb interaction.Reed et al. (2010) Without the possibility of inter-layer hopping, the polarization function of graphene and graphite are, apart from some multiplicative factor, the same. As we have seen in Sec. IV, this simplification does not capture the low energy part of the spectrum, with some finite spectral weight due to low energy inter-band transitions between parabolic-like bands.

At an energy of the order of one observes the peak due to transitions between electrons from the Van Hove singularity of the occupied band to the singularity of the empty band. For SLG, the peak is split into two peaks if the wave-vector points in the -K direction (as it is the case here), the separation of which increases with the modulus . However the amplitude of these peaks is highly suppressed from SLG to MLG. Finally, one observes in Fig. 14(b) that for higher values of , as the one used here, there is a redshift of the peak of at the energy with respect to the Dirac cone approximation. This is due to trigonal warping effects, which are beyond the continuum approximation. Summarizing, we find two effects that lead to a global contribution to the polarization at low energies: one is the contribution to the spectral weight due to inter-layer hopping in MLG, and the other is the redshift of the peaks at due to trigonal warping effects. Notice that, because we are studying the non-interacting polarization function, no excitonic effects are present in the results of Fig. 14(a)-(b).

Once the polarization function is known, we compute the response function at the RPA level, as shown in Fig. 14(c)-(d). Again, we find a redshift of the position of the peaks as we decrease the number of layers. The different position of the peaks is due to the different contribution of inter-layer electron-electron interaction for each case, as well as to the different value of as given by Eq. (21). Our results agree reasonably well with those of Ref. Reed et al., 2010.

Finally, we calculate the renormalization of the fine structure constant due to dynamic screening associated to the inter-band transitions from the valence band. For this, and in analogy with Ref. Reed et al., 2010, we define

(23) |

The results for the modulus for SLG and MLG are shown in Fig. 15. In this plot we have used the value of the Fermi velocity valid near the Dirac point, i.e., . Therefore, we emphasize that these results should be reliable only at low energies. At and for the smallest wave-vector we can access (), RPA predicts for SLG, which is considerably higher than the value estimated in Ref. Reed et al., 2010: . However, the results that we obtain for MLG are much closer to this value: we find that for graphite, only slightly higher (a factor of 2) than the experimental results of Ref. Reed et al., 2010, which are actually obtained from graphite.

## Vi Conclusions

In conclusion, we have studied the excitation spectrum of single- and multi-layer graphene using a full -band tight-binding model in the random phase approximation. We have found that, for MLG, the consideration of inter-layer hopping is very important to properly capture the low energy region () of the spectrum. This, together with trigonal warping effects, lead to a finite contribution to the spectral weight at low energies as well as a redshift of the peaks with respect to the Dirac cone approximation. We have also studied the high energy plasmons which are present in the spectrum of SLG and MLG at an energy of the order of eV and which are associated to the enhanced DOS at the Van Hove singularities of the -bands. The energy of the -plasmon depend also on the orientation between adjacent layers, and we find that, for a given wave-vector, the energy of the mode for ABC-stacked MLG is redshifted with respect to the corresponding energy of ABA ordering. This difference is higher as we increase the number of graphene layers of the system.

The effect of disorder has been considered by the inclusion of a random on-site potential and by a renormalization of the nearest neighbor hopping. Both kinds of disorder lead to a redshift of the and peaks of the non-interacting excitation spectrum and to a smearing of the Van Hove singularities. The position of the -plasmons is unaffected by disorder, although the height of the absorption peaks is reduced as compared to the clean limit.

Finally, we have compared our results to some recent experiments. Our calculations for the loss function show a redshift of the SLG mode with respect to graphite, and compare reasonably well with experimental EELS data.Eberlein et al. (2008); Gass et al. (2008) Furthermore, we also obtain good agreement with the IXS results for the response function obtained in Ref. Reed et al., 2010. We obtain a static dielectric function which grows with the number of layers of the system. In the long wavelength and limit, the dynamically screened fine structure constant is found to be highly reduced from graphene to graphite. The value that we find for a MLG in the RPA, without considering any excitonic effects, is about two times larger than the one estimated in Ref. Reed et al., 2010 for graphene. More accurate results could be obtained going beyond single-band RPA,van Schilfgaarde and Katsnelson (2011) which is beyond the scope of this work.

## Vii Acknowledgement

The support by the Stichting Fundamenteel Onderzoek der Materie (FOM) and the Netherlands National Computing Facilities foundation (NCF) are acknowledged. We thank the EU-India FP-7 collaboration under MONAMI.

## References

- For a recent review on electron-electron interaction in graphene, see et al. (2010) For a recent review on electron-electron interaction in graphene, see, V. Kotov, B. Uchoa, V. M. Pereira, A. H. Castro-Neto, and F. Guinea, (2010), arXiv:1012.3484 .
- Kramberger et al. (2008) C. Kramberger, R. Hambach, C. Giorgetti, M. H. Rümmeli, M. Knupfer, J. Fink, B. Büchner, L. Reining, E. Einarsson, S. Maruyama, F. Sottile, K. Hannewald, V. Olevano, A. G. Marinopoulos, and T. Pichler, Phys. Rev. Lett. 100, 196803 (2008).
- Eberlein et al. (2008) T. Eberlein, U. Bangert, R. R. Nair, R. Jones, M. Gass, A. L. Bleloch, K. S. Novoselov, A. Geim, and P. R. Briddon, Phys. Rev. B 77, 233406 (2008).
- Gass et al. (2008) M. H. Gass, U. Bangert, A. L. Bleloch, P. Wang, R. R. Nair, and A. Geim, Nature Nanotechnology 3, 676 (2008).
- Reed et al. (2010) J. P. Reed, B. Uchoa, Y. I. Joe, Y. Gan, D. Casa, E. Fradkin, and P. Abbamonte, Science 330, 805 (2010).
- Mak et al. (2011) K. F. Mak, J. Shan, and T. F. Heinz, Phys. Rev. Lett. 106, 046401 (2011).
- Bostwick et al. (2010) A. Bostwick, F. Speck, T. Seyller, K. Horn, M. Polini, A. H. MacDonald, and E. Rothenberg, Science 328, 999 (2010).
- Liu et al. (2008) Y. Liu, R. F. Willis, K. V. Emtsev, and T. Seyller, Phys. Rev. B 78, 201403 (2008).
- Shung (1986) K. W. K. Shung, Phys. Rev. B 34, 979 (1986).
- Wunsch et al. (2006) B. Wunsch, T. Stauber, F. Sols, and F. Guinea, New Journal of Physics 8, 318 (2006).
- Hwang and Sarma (2007) E. H. Hwang and S. D. Sarma, Phys. Rev. B 75, 205418 (2007).
- Polini et al. (2008) M. Polini, R. Asgari, G. Borghi, Y. Barlas, T. Pereg-Barnea, and A. H. MacDonald, Phys. Rev. B 77, 081411 (2008).
- Pyatkovskiy (2009) P. K. Pyatkovskiy, Journal of Physics: Condensed Matter 21, 025506 (2009).
- Gamayun (2011) O. V. Gamayun, (2011), arXiv:1103.4597 .
- Tudorovskiy and Mikhailov (2010) T. Tudorovskiy and S. A. Mikhailov, Phys. Rev. B 82, 073411 (2010).
- Roldán et al. (2009) R. Roldán, J.-N. Fuchs, and M. O. Goerbig, Phys. Rev. B 80, 085408 (2009).
- Gangadharaiah et al. (2008) S. Gangadharaiah, A. M. Farid, and E. G. Mishchenko, Phys. Rev. Lett. 100, 166802 (2008).
- Khveshchenko (2001) D. V. Khveshchenko, Phys. Rev. Lett. 87, 246802 (2001).
- Aleiner et al. (2007) I. L. Aleiner, D. E. Kharzeev, and A. M. Tsvelik, Phys. Rev. B 76, 195415 (2007).
- Gamayun et al. (2009) O. V. Gamayun, E. V. Gorbar, and V. P. Gusynin, Phys. Rev. B 80, 165429 (2009).
- Wang et al. (2010) J. Wang, H. A. Fertig, and G. Murthy, Phys. Rev. Lett. 104, 186401 (2010).
- González (2010) J. González, Phys. Rev. B 82, 155404 (2010).
- Stauber et al. (2010) T. Stauber, J. Schliemann, and N. M. R. Peres, Phys. Rev. B 81, 085409 (2010).
- Hill et al. (2009) A. Hill, S. A. Mikhailov, and K. Ziegler, Europhysics Letters 87, 27005 (2009).
- Yang et al. (2009) L. Yang, J. Deslippe, C.-H. Park, M. L. Cohen, and S. G. Louie, Phys. Rev. Lett. 103, 186802 (2009).
- Koshino (2010) M. Koshino, Phys. Rev. B 81, 125304 (2010).
- Castro-Neto et al. (2009) A. H. Castro-Neto, F. Guinea, N. M. R. Peres, K. Novoselov, and A. K. Geim, Rev. Mod. Phys. 81, 109 (2009).
- Kubo (1957) R. Kubo, Journal of the Physical Society of Japan 12, 570 (1957).
- Yuan et al. (2010a) S. Yuan, H. De Raedt, and M. I. Katsnelson, Phys. Rev. B 82, 115448 (2010a).
- Hams and De Raedt (2000) A. Hams and H. De Raedt, Phys. Rev. E 62, 4365 (2000).
- Giuliani and Vignale (2005) G. F. Giuliani and G. Vignale, Quatum Theory of the Electron Liquid (CUP, Cambridge, 2005).
- Stauber et al. (2008) T. Stauber, N. M. R. Peres, and A. K. Geim, Phys. Rev. B 78, 085432 (2008).
- Shyu and Lin (2000) F. L. Shyu and M. F. Lin, Phys. Rev. B 62, 8508 (2000).
- Das Sarma and Quinn (1982) S. Das Sarma and J. J. Quinn, Phys. Rev. B 25, 7603 (1982).
- Wehling et al. (2011) T. O. Wehling, E. Şaşıoğlu, C. Friedrich, A. I. Lichtenstein, M. I. Katsnelson, and S. Blügel, Phys. Rev. Lett. 106, 236805 (2011).
- Wehling et al. (2010) T. O. Wehling, S. Yuan, A. I. Lichtenstein, A. K. Geim, and M. I. Katsnelson, Phys. Rev. Lett. 105, 056802 (2010).
- Yuan et al. (2010b) S. Yuan, H. De Raedt, and M. I. Katsnelson, Phys. Rev. B 82, 235409 (2010b).
- Roldán et al. (2010) R. Roldán, M. O. Goerbig, and J.-N. Fuchs, Semicond. Sci. Technol. 25, 034005 (2010).
- Stauber and Gómez-Santos (2010) T. Stauber and G. Gómez-Santos, Phys. Rev. B 82, 155412 (2010).
- Hobson and Nierenberg (1953) J. P. Hobson and W. A. Nierenberg, Phys. Rev. 89, 662 (1953).
- Adler (1962) S. L. Adler, Phys. Rev. 126, 413 (1962).
- Pellegrino et al. (2010) F. M. D. Pellegrino, G. G. N. Angilella, and R. Pucci, Phys. Rev. B 82, 115434 (2010).
- Guinea (2007) F. Guinea, Phys. Rev. B 75, 235433 (2007).
- Zhang et al. (2008) L. M. Zhang, Z. Q. Li, D. N. Basov, M. M. Fogler, Z. Hao, and M. C. Martin, Phys. Rev. B 78, 235408 (2008).
- van Schilfgaarde and Katsnelson (2011) M. van Schilfgaarde and M. I. Katsnelson, Phys. Rev. B 83, 081409 (2011).