Changing character of electronic transitions in graphene: From single particle excitations to plasmons

In this paper we clarify the nature of $\pi$ and $\pi+\sigma$ electron excitations in pristine graphene. We clearly demonstrate the continuous transition from single particle to collective character of such excitations and how screening modifies their dispersion relations. We prove that $\pi$ and $\pi+\sigma$ plasmons do exist in graphene, though occurring only for a particular range of wavevectors and with finite damping rate. The particular attention is paid to compare the theoretical results with available EELS measurements in optical ($\mathrm{Q\approx 0}$) and other ($\mathrm{Q\neq 0}$) limits. The conclusions, based on microscopic numerical results, are confirmed in an approximate analytical approach.


INTRODUCTION
Plasmon spectra of a pristine single layer graphene were first obtained in Ref. [1], where the authors observed two distinct structures which they attributed to the so-called π and π + σ plasmons. They observed that these two plasmonic modes were redshifted in comparison to the corresponding modes in the bulk graphite [2][3][4], due to the reduction of macroscopic screening when going from graphite to graphene [5]. The early momentumdependent theoretical and experimental measurements observed linear dispersion of this π plasmon in graphene [6][7][8][9], which differs from the Q 2 dispersion reported in graphite [2,3,5,10].
Recently a resolute claim was made [11] that the previously accepted attribution of the two strong structures in the graphene excitation spectra was wrong, and the π and π + σ plasmons are in fact strong single particle (SP) π → π * and σ → σ * excitations, respectively, with a characteristic Q 2 excitation energy dependence. Another group [12] found strong evidence for 2D plasmon character of π and σ electron excitations, based on the electron energy loss spectroscopy (EELS) experiment showing the √ Q dependent dispersion. Even taking into account possible uncertainties arising from experimental difficulties in EELS measurements for low Q values, it is obvious that this apparent controversy deserves to be analysed and resolved.
In this paper we solve this problem using both a numerical method and analytic arguments, providing a rigorous method of determining collective vs. single particle excitation character in solids, and apply it to analyse π and σ electron excitations in a self-supporting monolayer of pristine graphene. We find that the character of π → π * (and σ → σ * ) transition changes, depending on the wavevector Q. For small Q ≈ 0 these are unscreened single particle transitions, but with increasing Q they acquire collective character as the dynamical screening mechanism becomes more efficient. This explains the gradual change from the Q 2 dependence of excitation energies near Q ≈ 0 to the quasi-linear dependence at larger Q. And finally, for even larger Q the collective nature of this modes in graphene is suppressed and they again emerge as the single particle excitations. Same kind of dispersion is observed in [13][14][15], but the authors did not analyse it in detail. Although the described dispersion seems like the characteristic √ Q dependence of the 2D plasmon, we show that this cannot be true because of the complex nature of this mode.
The described analysis is quite general, and its application to graphene provides a very nice illustration how an electronic process can change its character from an interband single particle transition to a collective mode as the dynamical screening takes over with the increasing wavevector Q.
In Sec. II we describe the derivation of the electronic excitation spectra S(Q, ω) in terms of the dielectric tensor E GG (Q, ω), using the method of Ref. [20], and define the dynamical screening factor D(Q, ω). In Sec. III we calculate numerically the macroscopic dielectric functions Im E M (Q, ω), Re E M (Q, ω) and −Im 1/E M (Q, ω), as well as the dynamical screening factor D(Q, ω), and discuss the excitation spectra. In Sec. IV we summarize the results and their relation to previous experimental and theoretical work.

II. FORMULATION OF THE PROBLEM
The first part of the calculation consists of determining the Kohn-Sham (KS) ground state of graphene and the corresponding wave functions and energies. For the unit cell constant we use the experimental value of a = 4.651 a.u. [16], and we separate the graphene layers with the distance L = 5a. For calculating KS wave functions and energies we use a plane-wave selfconsistent field DFT code (PWSCF) within the QUAN-TUM ESPRESSO (QE) package [17]. The core-electron interaction was approximated by the norm-conserving pseudopotentials [18], and the exchange correlation (XC) potential by the Perdew-Zunger local density approximation (LDA) [19]. To calculate the ground state electronic density we use 30 × 30 × 1 Monkhorst-Pack K-point mesh of the first Brillouin zone (BZ) and for the plane-wave cut-off energy we choose 50 Ry.
Using the wave functions and energies obtained in the described way we perform the calculation of the electronic excitation spectra within the random phase approximation (RPA). In the quasi-2D systems such as graphene, with the electronic density in the region 0 < z < L, the spectral function of electronic excitations can be defined as [20] where the dielectric matrix in the RPA is given by The noninteracting charge-charge response function is given in matrix form as The V = S × L is the normalization volume, S is the normalization surface and f n (K) = θ[E F − E n (K)] is the Fermi-Dirac distribution at T = 0. In this summation we use 201×201×1 K-mesh sampling and up to 70 electronic bands. For the broadening parameter η we use 0.05 eV. Matrix elements in (3) have the form where Q is the momentum transfer vector parallel to the x-y plane, G = (G , G z ) are 3D reciprocal lattice vectors and r = (ρ, z) is a 3D position vector. Wave functions Φ nK (r) are KS wave functions from the ground state calculation and E n (K) are the corresponding energies. In this approach the superlattice consists of periodically repeated layers of graphene such that the charge densities of adjacent layers do not overlap. We restrict our consideration to one layer placed in the region 0 < z < L, where the interaction with the adjacent layers is avoided by allowing Coulomb interaction between charge oscillations only within this region. This is done by integrating the Dyson equation for χ within the limits of 0 < z < L [20][21][22]. The resulting Coulomb interaction matrix elements have the explicit form with Similar approach was also carried out for layered structures in Ref. [23]. In our analysis of electronic excitation spectra we will use the macroscopic dielectric function which is defined as and includes the crystal local field effects in the perpendicular, though not in the parallel direction. By this we mean that we put G || = 0, while leaving the reciprocal lattice vectors in z direction, G z , which is justified by the fact that the parallel local field effects are not so important for describing the surface plasmons [24]. To get a well converged spectra we use 71 G z vectors. With this macroscopic dielectric function the spectral function (1) can be written as where we have defined the dynamical screening factor D(Q, ω) and simplified the notation with E 1 ≡ Re E M and E 2 ≡ Im E M . For a vanishing screening in the system we have that D → 1, while for the screened D = 1. In the    first case the spectral function S is equal to E 2 and all the structures in spectra have purely SP character.
To analyse the electronic excitation spectrum of a 2D material one can also use the dielectric function within the tight-binding approximation (TBA) [25][26][27]. Instead of the KS wave functions and energies, here one uses the states and energies of the TBA hamiltonian. If we consider graphene beyond the Dirac cone approximation then the 2D charge-charge response function is given by where the band index µ = −1 represents the occupied π band and µ = 1 the unoccupied π * band in graphene. For the numerical calculation of (8) we use 600 × 600 × 1 Kpoint mesh and the broadening parameter η = 0.05 eV. The TBA band energies are E µ (K) = µγ|g(K)| with the hopping parameter γ ≈ 2.02 eV [28,29], while the hopping function is given by where a is the lattice constant of graphene. Here the spectral function is also defined as in (1) and (7), but the dielectric function does not have crystal local field effects included and can be written as III.

RESULTS AND DISCUSSION
In order to better understand the spectrum of SP electronic excitations we will first analyze the graphene band structure. Fig. 1 shows the graphene band structure along the M − M (black lines) and M − Γ (blue lines) directions of the BZ, which are relevant directions when the wavevector of external perturbation, Q, is in the ΓM direction. We see that π electrons exhibit two different kinds of interband transitions. The first is attributed to transitions between two almost dispersionless π bands along the M − M direction, as denoted by black thick arrow in Fig. 1. The second kind of π interband transitions are attributed to transitions along the M − Γ direction, as denoted by a thick blue arrow in Fig. 1. We shall see that the π plasmon can be formed from the latter transitions when they are dynamically screened. Two other transitions are between occupied and unoccupied σ bands. They can be divided into σ → σ * 1 and σ → σ * 2 transitions [30], as denoted by thin blue and black arrows in Fig.  1. Moreover, the σ 1 and σ 2 plasmons, usually treated as one π + σ plasmon, originate from these transitions. FIG . 3: a) Energies of π → π * interband transitions (blue) and the π plasmon (red) as functions of the wavevector Q along ΓM direction obtained from DFT-RPA method. Black unfilled diamonds represent the experimental data from Ref. [12]. b) Same as in a) for σ → σ * 1 and σ → σ * 2 (blue) transitions and σ1 and σ2 plasmons (red). The point Q = 0 is treated separately [31]. c) Same as in a) but showing the results obtained with the TB-RPA method. In addition the points where ε1(Q, ω) = 0 are shown by the black dashed line and the corresponding √ Q fit with the orange line. Inset: Real (green) and imaginary (dashed blue) parts of ε(Q, ω) and S(Q, ω) (black) for Q = 0.01 a.u. as obtained with TB-RPA method.
means that the mentioned transitions are not screened, i.e. the peaks in S(Q, ω) are pure SP excitations which appear at the same energies as peaks in E 2 (Q, ω) [11,25]. Blue dots in Figs. 3a and 3b show the energies of the peaks in E 2 (Q, ω), and red dots the energies of the peaks in S(Q, ω) as functions of the wavevector Q. Blue points show characteristic Q 2 dispersion of π and σ SP transitions. It can be clearly seen that for small Q peaks in E 2 (Q, ω) and peaks in S(Q, ω) coincide and follow the same Q 2 dependence which confirms their SP character. This quadratic dispersion of SP excitations is a result of the π and σ band structure around the saddle point M, as sketched in Fig. 1. As Q increases the screening factor D(Q, ω) increases, enhancing the spectral weight of the peaks and moving them to higher energies, i.e. away from the initial π → π * and σ → σ * energies (Figs. 2b,c). This is also visible in Figs. 3a and 3b which represent the gradual modification of the initial SP transitions into collective excitations as the dynamical screening becomes more efficient. In this regime (Q 0.03 a.u.) one can with confidence treat these excitations as π and π + σ plasmons, though their broad spectral shapes indicate the presence of Landau damping. Diamonds in Fig. 3a show the energies of π plasmon peaks in the measured spectra [12]. We see very nice agreement with our theoretical calculation throughout the whole energy region. Therefore, the pronounced spectral structures which appear for Q ≈ 0 (e.g. in optical absorption spectra) have purely SP character, while spectral structures which appear for finite Q (e.g. in EELS) represent collective excitations or plasmons. For Q ≈ 0.1 a.u. the screening becomes most efficient and E 1 (Q, ω) approaches its lowest value. Accordingly, in Fig. 3a we can see the largest shift of the π plasmon energy compared to the π → π * transition energy. For even higher wavevectors Q this shift becomes smaller, the plasmon energy slowly approaches the unscreened π → π * transition energies, but the shape of the plasmon peak remains almost unchanged (Figs. 2e and 2f).
These effects can also be observed in Fig. 4 which shows energy dependence of the real part of the dielectric function obtained with DFT-RPA and TB-RPA methods. In the TB-RPA approximation where only π electrons participate in the screening ε 1 (Q, ω) crosses zero for all Q above some minimum value where the resulting dispersion relation reaches the SP continuum, as shown in Fig. 3c. So, it is obvious that the π electrons in TB-RPA for higher Q's behave like a 2D electron gas showing the √ Q dispersion relation. However, for Q → 0 the excitation energy approaches the finite value at the upper boundary of the π → π * continuum (≈ 4 eV) with the Q 2 dependence (red curve in Fig. 3c), and does not follow the √ Q line to zero (orange curve in Fig. 3c). The reason for this is high interband Landau damping and general reduction of macroscopic screening for finite systems in the low Q region [5,32]. This is all in qualitative agreement with our resluts, but we can conclude that √ Q dispersion is obtained only if we neglect other electronic transitions in graphene. So in the DFT-RPA calculation, including full graphene band structure, i.e. σ electrons, a number of zeros of E 1 (Q, ω) is strongly reduced, and if we further include finite lattice effects, i.e. LFE, E 1 (Q, ω) never cross zero in the whole (Q, ω) range (Fig. 4). The spectra still shows well defined though broader peaks, but the plasmon energies are in a much better agreement with the experimental results than the TB-RPA calculation. Also, for very large Q the π plasmon energy should approach the SP π → π * transition energy, which is not achieved in TB-RPA.
It might seem that our TB-RPA results qualitatively agree with the theoretical interpretation of the measured data in Ref. [12]. There the authors demonstrated that the graphene π plasmons represent the in-plane charge density oscillations which led them to the conclusion that they behave like plasmons in a 2D electron gas with ∼ √ Q dispersion relation. Now we clearly see that this conclusion is indeed partially true, but at the same time we see its limitations. First of all, even the TB-RPA re-sults confirm ∼ √ Q dispersion only for bigger Q's, while the authors missed the fact that for small Q's it shows a Q 2 dependence. This is due to absence of dynamical screening and collective mode changes its character to single particle excitations. Possible reason for this failure is because the first measured nonzero Q point is at Q ≈ 0.05 a.u. (or 0.1Å −1 ) above the Q 2 region, as can be clearly seen by observing the authors results in Fig. 3a. So after inclusion of all other points, for Q 0.05 a.u., one gets the impression that the dispersion has a square root behavior.
Apart from these deviations for small Q's there are other arguments which violate this simple 2D plasma attribution. As we already mentioned, the inclusion of realistic crystal and band structure (σ bands and LFE) substantially modifies the TB-RPA ∼ √ Q dispersion which for bigger Q's even becomes linear. This is also confirmed by the excellent agreement of our DFT-RPA dispersion and the experimental data shown in Fig. 3a.
The dispersion of the so-called π + σ plasmon is more complicated, as visible in Figs. 2 and 3b. There are in fact not one, but two modes, one originating from σ → σ * transitions (σ 1 plasmon) and the other from σ → σ * 2 transitions (σ 2 plasmon), as shown in Fig. 1. For low Q values their dispersion is shown in Fig. 3b, together with energies of unscreened SP transitions. While they generally follow the behaviour described above for π plasmons, from SP excitations to the Landau damped plasmons and back, they give much broader structures, and furthermore there are additional features arising from their mutual interference. For small Q σ → σ * 1 excitation dominates in intensity, but around Q ≈ 0.08 a.u. it becomes modified by σ → σ * 2 transitions which gain spectral weight, so the high energy spectra of graphene are dominated by σ 2 plasmons. For large Q 0.4 a.u. values the spectra again show unscreened SP excitations.
Therefore, as for the π plasmon, in Ref. [12] the authors make a hasty conclusion about the √ Q dispersion of σ plasmon, also not taking into account that two kinds of σ → σ * excitations exist close in energy and influence each other. In experimental papers Ref. [13,14] the authors also obtained π plasmon dispersion relation which is in accordance with our theoretical results and conclusions, however they are also prone to describe it as √ Q behaviour, with a help of the hydrodinamic model [33,34].
Our above analysis partially agrees with the results recently published in Ref. [11], though not entirely, since its authors claim that the π and π + σ plasmons in graphene do not exist at all. Namely, they derived their conclusions by analysing the dielectric function in the Q → 0 (optical) limit which is indeed the region where the π and σ excitations behave as unscreened SP transitions. However, as we showed, this is not the case for higher Q's. One of the methods to demonstrate that π and π + σ excitations represent self-sustaining charge density oscillations is to induce them by some external perturbation and see how they behave. If these charge oscillations survive at least one period, after the perturbation is being switched off, then they represent collective modes. In Ref. [21] the charge density oscillations in graphene are driven by suddenly created point charge. The authors have identified two periods of oscillations whose frequencies are associated with π and π + σ plasmon frequencies. In Ref. [22] the charge density in graphene is driven by a point charge moving with constant velocity parallel to the graphene surface. There appear several rows of bow waves whose wavelength, depending on the speed of the point charge, was associated with the excitations of π or π + σ plasmons. These two observations undoubtedly confirm that π and π + σ plasmons indeed exist.
A brief analytic discussion could help to clarify the shape of numerically calculated spectra. The spectrum in the absence of screening (D = 1) will have maxima at the energies of π → π * and σ → σ * SP transitions, as in Fig. 2a, with Q 2 dispersion coming from the shape of the π and σ bands. As the screening term D(Q, ω) increases with larger Q, these SP peaks will be enhanced and shifted to the maxima Ω of this term, given by the condition or explicitly where the prime denotes the derivative with respect to ω. This systematic transition from the SP peaks to the maxima of the screening function can be seen in Figs. 2b-2f, and are the proof of the increasingly collective character of these excitations for higher Q. At the same time the dispersion deviates from the initial Q 2 dependence, as seen in Figs. 3a and 3b. From (12) we can explicitly see that the peak positions are not given by the condition as would be expected for an ideal collective mode, and indeed, in our full calculation including LFE (13) is not satisfied for any (Q, ω). Nevertheless we find broad but well defined spectral peaks corresponding to plasmons in graphene, though restricted to the region of finite wavevectors and Landau damped. We can also derive the condition (12) assuming that the spectrum (10) has a resonant form, with the maximum at the complex pole of the dielectric function: by expanding around Ω with Γ Ω. Dielectric function near the resonance energies Ω shows an interesting behaviour. By inspection of numerical results for finite Q we see which is also in agreement with the condition (12), valid in the resonance or plasmon region. So the properties (15) can be connected with the collective character of the excitations in this region.

IV. CONCLUSION
In this paper we have presented the method to determine the character of electronic excitations in solids and applied it to the π and σ electron excitations in a single layer of pristine graphene. Analysing the energy and momentum dependence of the dielectric function components, which determine the excitation spectra, and especially the dynamical screening factor D(Q, ω), we have demonstrated that π → π * and σ → σ * transitions show two types of behaviour. For small wavevector Q ≈ 0 they indeed behave like pure interband (SP) transitions, because the screening is completely absent, but as Q increases and the screening becomes prominent, they rather suddenly acquire a collective character. In this Q region one can indeed say that π and π + σ plasmons exist in graphene, though they are always Landau damped and appear as broad structures in the spectra. We also demonstrated that because of strong Landau damping, the dielectric function never crosses zero, even in the region of collective excitations. These conclusions, based on exact numerical results, are also confirmed by a brief analytic discussion.
We have partially confirmed the seemingly conflicting results of Ref. [11] and Refs. [12,14]. The claim in [11] that π and σ excitations are unscreened SP transitions is indeed correct in the low Q region. This result is expected because in 2D materials, like graphene, macroscopic screening E 1 is significantly reduced in comparison with the bulk system like graphite [5,32]. Our results agree very well with the experimental data of Ref. [12], however we partially disagree that π and σ electrons behave as 2D electrons with √ Q dispersion relation. We showed that π and σ excitations change their character: in Q ≈ 0 region they have SP character with Q 2 dispersion, as Q increases they acquire collective character with √ Q like dispersion, and for even higher Q's they again acquire SP character with linear dispersion. In this way we have presented a complete description of this quite intriguing and previously controversial problem, explaining the changing character of electronic excitations in graphene.
By presenting our conclusions we would like to emphasize once again the importance of the in-depth interpretation of experimental measurements and how they relate to the theoretical results. So, to understand the nature of plasmon dispersion it is crucial to perform measurements for all values of the wavevector Q (e.g. with optical absorption spectra and EELS) and to give a careful theoretical interpretation in each of the limits (i.e. Q ≈ 0 and Q = 0).