Electronic excitations in quasi-2D crystals: What theoretical quantities are relevant to experiment?

The ab initio theory of electronic excitations in atomically thin [quasi-two-dimensional (Q2D)] crystals presents extra challenges in comparison to both the bulk and purely 2D cases. We argue that the conventionally used energy-loss function $-$Im $1/\epsilon({\bf q},\omega)$ (where $\epsilon$, ${\bf q}$, and $\omega$ are the dielectric function, the momentum, and the energy transfer, respectively) is not, generally speaking, the suitable quantity for the interpretation of the electron-energy loss spectroscopy (EELS) in the Q2D case, and we construct different functions pertinent to the EELS experiments on Q2D crystals. Secondly, we emphasize the importance and develop a convenient procedure of the elimination of the spurious inter-layer interaction inherent to the use of the 3D super-cell method for the calculation of excitations in Q2D crystals. Thirdly, we resolve the existing controversy in the interpretation of the so-called $\pi$ and $\pi+\sigma$ excitations in monolayer graphene by demonstrating that both dispersive collective excitations (plasmons) and non-dispersive single-particle (inter-band) transitions fall in the same energy ranges, where they strongly influence each other.


Introduction
Quasi two-dimensional (Q2D) crystals (the systems periodic and infinite in two dimensions while of atomic-size thickness) are presently attracting much of attention because of their potential in future nano-technologies. In particular, the electronic excitations in Q2D crystals, such as electron-hole pairs generation, plasmons, excitons, etc., are of great importance, since they determine the response to external actions, and, ultimately, shape the materials' functionality (see, e.g., [1] for a recent review). At the same time, the theoretical approaches as well as the computational methods to deal with the excitations in Q2D systems become much more involved compared to both the three-dimensional (3D) and purely 2D cases, and are still far from finally established.
Indeed, in both purely 3D and purely 2D cases (the latter referring to model systems of zero thickness), one usually introduces the dielectric function ε(q, ω), where q is the wave-vector of the corresponding dimensionality, and ω is the frequency. The definition of ε reads ‡ φ ext (q, ω) = ε(q, ω)φ tot (q, ω), where φ ext and φ tot are the externally applied and the total (external plus induced) scalar potentials of the electric field, respectively. The definition (1) is not, however, straightforwardly transferable to the Q2D case, which can be easily appreciated by considering that φ tot depends on the z (perpendicular to the layer) coordinate, even when φ ext is uniform in z direction. Therefore, we need to redefine the dielectric function (energy-loss function, conductivity, etc.) of a Q2D crystal in accordance to the system's structure. Importantly, the new definitions must be consistent with the quantities measured experimentally, the violation of which requirement having caused much of confusion in recent literature.
Secondly, the widest used at this time method of dealing with Q2D systems numerically is the super-cell geometry approach. This is to artificially replicate the system periodically in z direction, choosing the distance d between the layers large enough to prevent the inter-layer interactions. The resulting fictitious system is 3Dperiodic, and well developed 3D solid-state methods as well as the existing 3D software can be used to efficiently handle it. In the ground-state calculations, the interpretation of the super-cell results is straightforward: Provided that the inter-layer distance is such that the corresponding wave-functions do not overlap, properties of an isolated single layer coincide with those of a single super-cell.
However, if we are concerned with the dynamic response to electric field, the situation changes. In the case of a Q2D crystal, the electric field decays into vacuum as ∼ e −q|z| . Clearly, whatever large is the separation d between the layers, at sufficiently small q the inter-layer interaction persists, and results of the super-cell calculation cannot be literally transferred to a single layer. A simple method of taking d very large does not work, because small q-s are often of special interest, and the computational load grows rapidly with the increase of d. Therefore, we need to establish relations between the response functions of the fictitious 3D array system and those of the single layer of interest. The paramount requirement to these relations (test for sanity) is that they should produce the same response functions of a single layer from the differing 3D response functions calculated with different separations d.
In this paper, in the framework of the time-dependent density-functional theory (TDDFT), we construct a scheme satisfying the above requirements. As its practical application, we calculate and perform a detailed analysis of the excitation spectrum of a single-layer graphene in the 1-30 eV energy range. Resolving the recent controversy in the interpretation of the π and π + σ peaks as plasmons [2,3,4,5] or single-particle inter-band transitions [6], we show that, while the problem is far more complicated than thought before, there do exist prominent π and π + σ plasmons in graphene, although they overlap with the inter-band transitions. The isolation of the both kinds of excitations becomes possible by following the wave-vector dependence of spectra (spatial dispersion). We also show that much of the previous confusion was due to using unsuitable theoretical quantities for comparison with experiment, as well as to the uncritical use of the results of the super-cell calculations in application to Q2D systems.
This paper is organized as follows. In Sec. 2 we re-examine and revise the definitions of the response-functions (dielectric function, energy-loss function, conductivity) pertinent to measurements on Q2D crystals. Sec. 3 is devoted to the methods of obtaining the response-functions of Q2D crystals from the 3D super-cell calculations. In Sec. 4 we present results of calculations performed for graphene as an important example of a Q2D crystal, discuss the differences arising from the use of our theory, and compare with experiment. Section 5 contains conclusions, and in the Appendix we give the derivation of some formulas presented in Sec. 3. We use the atomic units (e 2 = = m e = 1) unless explicitly specified otherwise.

Linear response of quasi-2D crystals: Definition of quantities
Quasi-2D crystals are neither periodic three-dimensionally nor strictly two-dimensional within the mathematical plane: They rather have a finite atomic or nano-scale thickness. Accordingly, one has to exercise caution in the choice of quantities (dielectric function, conductivity, density-response function, etc.) which characterize the dielectric response, giving them rigorous, suitable for Q2D systems, meaning. The safe way to do this is to start from the general 3D case without assuming the periodicity (translational invariance) in any direction. The density-response function χ(r, r ′ , ω) is defined in the real-space representation as where φ ext is the scalar potential of the externally applied electric field and ρ is the charge-density induced in the system in response to φ ext . The inverse generalized dielectric function ε −1 (r, r ′ , ω) is defined as where φ is the total (external plus induced) potential, and it is expressed through χ as By Fourier-transform, Eqs. (2)-(4) can be equivalently written in the reciprocal-space representation, of which we will need only For Q2D crystals, we will use the mixed reciprocal-and real-space representation, in the xy plane and z direction, respectively. Then the definition (2) becomes where G and G ′ are the 2D reciprocal lattice vectors, and q is the 2D wave-vector in the first Brillouin zone. We will also need the expression for the induced potential We distinguish between different kinds of experimental setups, where it will be shown that the measurable quantities correspond to different definitions of the energyloss function.

Electron energy-loss spectroscopy
The principal experimental method of probing q -and ω-dependent response of Q2D systems is the Electron Energy-Loss Spectroscopy (EELS). We will consider separately the reflection and transmission EELS.

Reflection EELS.
In reflection EELS, the quantity, which characterizes the inelastic electron scattering, is the so-called g-function [7] §. Its definition reads: Let the external potential be Then asymptotically at z → +∞ With the use of Eqs. (6) and (7), the g-function can be written as [11,12,7] When q and ω are substituted with the incident electron's loss of the parallel momentum and energy, respectively, the differential cross-section of the inelastic electron scattering in the reflection geometry is determined by the energy-loss function −Im g(q , ω), i.e., where dΩ is the element of the solid-angle of the detection of the scattered electron, , p and p ′ are the momenta of the incident and scattered electron, respectively. Equation (10) expresses g through the density-response function χ.

Transmission EELS.
For the transmission EELS of an arbitrary (possibly, non-periodic) target, the differential cross-section of the inelastic electron scattering is proportional to the energy-loss function (see, e.g., [13]) We emphasize, that in the right-hand side of Eq. (12) ε −1 of Eq. (5) enters with q and q ′ both substituted with the same momentum transfer q.
We consider the case of normal incidence. Then and, since p is large, .
(15) § In this paper, we restrict ourselves to the dipole scattering regime [7]. Impact-scattering regime is more involved, requiring the solution of the inelastic scattering problem within the distorted-waves method [8,9,10], which is outside the scope of this paper.
Everywhere below, only the z dependence of quantities is written down explicitly. The implied r and t dependence is always e i(q ·r −ωt) .
Then, by Eqs. (12) and (5), and restricting ourselves to the momentum transfer q within the first Brillouin zone, we have Equation (16) is the valid energy-loss function to be used in the theory of the transmission EELS. We point out that it replaces − 4π q 2 Im 1 ǫ(q,ω) [14], the latter not being defined for Q2D crystal.
2.1.3. In-plane electronic transport. Let us calculate the conductivity and the energy dissipation in a Q2D crystal under the action of the uniform in z direction external potential Then, by Eq. (6), and for the 2D particle-density averaged in the xy plane Using the continuity equation where j 2D is the 2D current-density, we can write and since, by Eq. (17), where σ ext 2D is the 2D conductivity with respect to the external field, which is given by For the energy dissipation Q per unit time per unit cell area under the action of the external potential (17) we can write ¶ ¶ Because Q is the quadratic rather than linear property, we must temporarily return to the (r, t) representation.
where j(r, t) and E(r, t) are well defined 3D current-density and electric field, respectively, and A is the 2D unit cell of the Q2D crystal. Since carriers do not commit work on themselves, Eq. (24) can be rewritten as and simplified to Using Eqs. (6) and (17), we can rewrite Eq. (26) as The important conclusion of Sections 2.1.2 -2.1.3 is that the energy losses in the reflection EELS, transmission EELS, and the energy dissipation in the in-plane transport are governed by the quantities e q (z+z ′ ) χ 00 (z, z ′ , q , ω)dzdz ′ , e iqz(z ′ −z) χ 00 (z, z ′ , q , ω)dzdz ′ , and χ 00 (z, z ′ , q , ω)dzdz ′ , respectively, which, for a Q2D crystal, are clearly different quantities. We note that in the limiting case of the purely 2D (zero-thickness) crystal, χ(z, z ′ , q , ω) contains δ(z)δ(z ′ ), and, consequently, all these three quantities coincide.
2.1.4. Dielectric function of Q2D crystal. It has already been noted in Sec. 1 that the dielectric function cannot be rigorously introduced for a Q2D crystal in a simple multiplicative way via Eq. (1). This is, in the first place, due to the uncertainty at which z φ tot (z, q , ω) must be used. It is, however, convenient to keep the notion of the dielectric function on the intuitive level, resorting to an analogy with the purely 2D system. In the latter case where χ 2D G G ′ (q , ω) is the 2D density-response function and the relation holds where ǫ 2D (q, ω) is the dielectric function of the 2D system. By Eqs. (23) and (28) and, using Eq. (29) For a Q2D crystal, we define the dielectric function using Eq. (31), where in the right-hand side we substitute σ ext 2D of Eq. (23), which is a well defined quantity for the Q2D crystal. Therefore By comparing Eqs. (16), (10), and (32), we see that the inelastic scattering of electrons in EELS of Q2D system is not defined by the loss-function −Im 1 ǫ Q2D (q ,ω) . It will be shown that this difference is of a quantitative importance in the interpretation of the measurements on Q2D crystals.
Below we will see that it is convenient to have our results in the full reciprocal-space representation. Using the Fourier expansion in the − d 2 , d 2 interval, where at this point we consider d an arbitrary sufficiently large distance, we can write Then Eq. (16) can be written as Eq. (10) as and Eq. (32) as 1 ǫ Q2D (q , ω) = 1 + 2πd q χ 00,00 (q , ω).
Having established expressions for the observables via the density-response function of a single layer χ G G ′ (z, z ′ , q , ω), we turn to the methods of the calculation of the latter.

Excitation of Q2D system from the super-cell geometry calculation
By far the widest used at this time way of the calculation of the electronic response of atomically thin Q2D crystals is the super-cell geometry approach. This is, instead of considering a single layer [schematized in Fig. 1, a)], to artificially periodically replicate it in z-direction [ Fig. 1, b)], and study the resulting 3D periodic system in place of the original single-layer one. It is, however, known [15,16,17] that the 3D dielectric function obtained in the super-cell geometry has nothing to do with that of a singlelayer system of interest, unless a special procedure of extracting the 2D response from the 3D super-cell geometry calculation (i.e., the elimination of the spurious inter-layer interaction) is applied. In Appendix A we give a detailed derivation of the expression for the density-response function χ of the single-layer system fromχ -the density-response function of the array system. It reads in the matrix form + where the matrix C GG ′ is given by The density-response function of the array systemχ G G ′ (q , ω) can be routinely obtained with the 3D solid-state periodic codes such, e.g., as Elk, abinit, Wien2k, etc.. Whileχ is different for different layers' separations d, χ obtained through Eq. (38) does not depend on d, which is crucial for its being a true characteristic of a stand-alone + In Eq. (38), the full 3D reciprocal space representation for bothχ and χ is used, which is natural for the former, but may look artificial for the latter. We, however, note that an arbitrary function of z and z ′ can be expanded in the Fourier series in z, z ′ ∈ [− d 2 , d 2 ], and then, if necessary, the real-space representation can be retrieved by the inverse Fourier transform.
2D system * . We note that Eqs. (38)-(40) are valid not only in RPA [16], but also in TDDFT with (semi-)local f xc (such as, e.g., ALDA), when different layers interact by the Coulomb potential only.
We conclude this section by noting that the ignoring the fundamental difference between the single-layer response of a Q2D crystal and that of the array of those layers has led to the misinterpretation of the π and π + σ peaks as single-particle inter-band transitions rather than plasmons [6].  We have conducted the super-cell geometry calculation for the monolayer pristine graphene followed by the application of Eq. (38) for the extraction of the single-layer response (elimination of the spurious inter-layer interaction) from the 3D calculation. In the DFT super-cell calculation, we used the full-potential linear augmented plane-wave * Strictly speaking, this is

Results of calculations and discussion
does formally depend on d, which is of no physical consequence.  (FP-LAPW) code Elk [19]. The z-axis period of the super-cell was taken 20 a.u. The k-point grid of 512 × 512 × 1, 30 empty bands, and the damping parameter of 0.002 a.u. were used in both the ground-state and the linear-response calculations. The former was conducted within the local-density approximation for the exchange-correlation potential, while the latter was the random-phase approximation (RPA) one (i.e., the exchangecorrelation kernel f xc was set to zero). The dielectric matrix of the 3D system was of the size 55 × 55, which was inverted to obtain ε 3D (q , ω; d) with the local-field effects included both in the perpendicular and in-plane directions.
A clear criterion for an excitation to be classified as plasmon is the requirement of the real-part of the dielectric function to cross zero at the corresponding energy. Otherwise, if the peak at the energy-loss function is due to a peak at the imaginary part of the dielectric function, the excitation is a single-particle transition. In our case, experiment -Im g -Im 1/ε Q2D Figure 4. The reflection EELS energy-loss function −Im g(q , ω) of Eq. (36) (solid lines) and the energy-loss function −Im 1 εQ2D (q ,ω) of Eq. (37) (dashed lines) for monolayer (left) and bilayer (right) graphene in the energy-range of π (upper) and π + σ (lower) plasmons. Symbols are experimental reflection EELS of Ref. [18]. The theoretical spectra have been roughly normalized to the experimental ones.
starting from greater wave-vectors (from q ≈ 0.05 and 0.1 a.u., for π and π + σ peaks, respectively) the real part of the dielectric function does cross zero, and, thereby, from those q-s up the peaks are plasmons unambiguously. Are they also plasmons at smaller q-s, or they turn into single-particle transitions, although keeping full similarity to their higher-q plasmon counterparts ? To answer this question, we need to establish whether, although Re ǫ 2D does not cross zero at smaller wave-vectors, the features at the energyloss function are due to Re ǫ 2D approaching zero, or they are due to peaks at Im ǫ 2D (cf., [17]).
As can be seen in Fig. 2, middle panel, Im ε 2D has prominent non-dispersive (standing in place with the variation of q ) feature at ω ≈ 4 eV. A prominent dispersive peak in the energy-loss function (bottom panel) follows the position of not only zero (at greater values of q ) but also the minimum (at smaller q -s) of Re ε 2D in the range of 4 -6 eV, and must, therefore, in both cases be classified as π plasmon.
The situation is similar, though less transparent, in the case of the π + σ peak. Here, the broad feature at L(q , ω), which spreads from ω ≈ 11 eV upward, is due to the overlap of the single-particle and collective excitations. Up to ω ≈ 14 eV, the excitation is purely due to the single-particle transitions and is non-dispersive. Then, from 14 eV upward, a dispersive feature due to zero or small values of Re ε 2D begins.
Because of Re ε 2D having a low slope and remaining small in absolute value in a wide ω range, this feature becomes very broad. It, however, is clearly plasmon as well. Interestingly, in two points, at ω ≈ 8 and 14 eV, Re ε 2D becomes nearly unity for all values of q , so that all the curves in Fig. 2, upper panel, intersect at these two points. The physical reason for this is still to be understood.
In Fig. 3, the reflection EELS-related energy-loss function −Im g(q , ω), calculated by Eq. (36), is compared with the loss-function −Im 1 ε Q2D (q ,ω) of Eq. (37), the latter relevant to the energy dissipation in the in-plane transport . While the π features (∼ 5 eV) are similar for both loss-functions, the π + σ features are considerably different, the EELS loss-function having a richer structure in this energy range. This can be understood by noting that, according to Eq. (37), 1/ε is determined by the head element of the density-response matrix χ 0000 (q , ω), while, by Eq. (36), the g-function also includes contributions from χ 0Gz0G ′ z (q , ω) with nonzero G z and/or G ′ z . The role of the latter increases with the increase of ω, transitions into the higher bands becoming allowed, which leads to the differences in the spectra.
In Fig. 4 we compare our theory with reflection EELS experiment of Ref. [18] on monolayer and bilayer graphene. In the energy range of π plasmon (two upper graphs in Fig. 4), g-function and the energy-loss (L) function are practically identical and both compare to experiment rather well. On the contrary, in the range of π + σ plasmons, the g-and L-functions acquire differences for both monolayer and, particularly, for bilayer graphene (lower graphs). For bilayer graphene π + σ peak at g-function is narrower than at the L-function, the former being closer to experiment (lower right). Remarkably, for monolayer graphene both g-and L-functions find themselves in shear disagreement with experiment (lower left). Obviously, the main source of the disagreement between the theory and experiment can be the presence of the SiC (0001) substrate in the latter and the absence thereof in the former. The substrate may cause additional screening leading to the shift of plasmon peaks and change the dispersion law [20,21,22,23]. Moreover, the interaction of graphene with a substrate can cause the deformation of the graphene sheet, leading to the confinement of the plasmon oscillations between the ripples [24].
Let us briefly discuss the modifications to the theory necessary to include Q2D crystals supported on substrates. For the ab initio theoretical treatment of surfaces, either pure or adsorbates covered, the super-cell method is widely applied as well.
In variance with the stand-alone Q2D systems, for surfaces, rather thick slabs of the material alternated with the similar thickness vacuum slabs are used. Equations (38)-(40) are readily applicable in this case too, producing as an output the density-response function of a single slab. This slab having two surfaces is not, however, the semi-infinite system of interest. Neither the array of such slabs, from which we have started, is. Which system, the array of slabs or a single slab, is better to model a Q2D crystal on top of a semi-infinite substrate must be decided in each particular case. An alternative method may be the inclusion of the substrate by means of the effective background dielectric function after having calculated the response of a stand-alone Q2D system. Although it may be useful in some situations, this would necessarily be a phenomenological approach, probably inconsistent with the rigorousness of the calculation for the Q2D crystal alone. Finally, the ultimate theory will be that asymptotically matching Q2D crystal on the surface with the 3D periodic bulk of the substrate, which is a very challenging task.
The neglect of the many-body exchange-correlation effects [25,26], both static and dynamic, may further contribute to the discrepancy. Both the accurate inclusion of the substrate in the theoretical setup and the inclusion of the many-body effects by means of using nonlocal exchange-correlation kernel of TDDFT remain the major challenges of the present-day theory of Q2D materials. Figure 5 shows the dispersion of π and π + σ plasmons with the q vector varied in the ΓM direction calculated in the EELS setup. The π + σ plasmon has a complex structure and is split into distinct sub-excitations.
As noted in Ref. [6], in a 2D system there can be no plasmon at q = 0. In the ordinary 2D electron gas [27] and in doped graphene in the case of 2D plasmon due to electrons in the π * band [28], this is ensured by the plasmon disappearance via its energy tending to zero as √ q when q → 0. However, as can be seen from Fig. 2, lower panel, for π and π + σ plasmons in graphene the same realizes in a different way: The energy positions of these plasmons converge to finite values (4 and 14 eV for π and π + σ plasmons, respectively) as q → 0, while the amplitudes of the corresponding peaks go to zero. In this regard, as in many others, graphene is also special compared to the ordinary 2D electron gas. The zero limit of the plasmons' amplitude can be understood from simple arguments. Indeed, quite generally, in a 2D crystal, the relation between the dielectric function ε 2D and the conductivity σ 2D is Since in the single-layer graphene, σ 2D (q = 0, ω) is finite [and known [29] to tend to e 2 /4 as ω → 0], ε 2D (q = 0, ω) is unity identically, leading to the zero energy-loss function. Our calculated conductivity of the monolayer pristine graphene at q = 0 versus the frequency is shown in Fig. 6, where the static limit of 1/4 is shown by the dotted line. Re σ 2D Im σ 2D Figure 6. Frequency-dependent conductivity of pristine monolayer graphene at q = 0. When ω → 0, the conductivity must converge to its analytic two-bands model value of 1/4, which it does, except for too small ω (at the very tip of the Dirac cone), where the accuracy of the calculation is limited by the k-point grid.

Conclusions
We have identified and dealt with the extra challenges the ab initio theory of electronic excitations in quasi-2D crystals presents in comparison to both the bulk and purely 2D cases. In particular, we have re-examined the problem of the correspondence between the theoretically calculated quantities and the observables in the measurements on quasi-2D crystals and found that the energy-loss function −Im 1 ǫ(q ,ω) , conventionally used for the interpretation of the EELS data, is not, generally speaking, the right quantity to be compared with this kind of the experiment. Instead, the EELS-related energy-loss functions, specific for both transmission and reflection experimental setups, must be used, which has been shown to be both qualitatively and quantitatively different in the case of quasi-2D systems. In the limit of the zero crystal thickness, our theory reduces to the conventionally used one.
We have addressed the problem of the classification of the π and π + σ peaks in monolayer graphene. By the use of the accurate procedure of extracting the dielectric function of a 2D crystal from the 3D super-cell geometry calculation, we have shown conclusively that there exist prominent π and π + σ collective exitations (plasmons) in graphene, although they are accompanied by interband transitions in the close energy ranges. Plasmons and single-particle interband transitions can be distinguished from each other provided the wave-vector-resolved analysis is performed.
We have also demonstrated the importance of the correct interpretation of the super-cell geometry calculation results, the latter taken simplistically, can lead to erroneous qualitative conclusions for materials of reduced dimensionality. We expect that the present theory of the electronic excitations in Q2D crystals, by correctly accounting for the atomic or nano-thickness of the systems studied, will further boost the theoretical and experimental work, improving the methods of comparison between the two.
Using Eqs. (A.8) and (A.9), the definition ofχ as the density-response function of the array system n G (q , ω) = G ′χ GG ′ (q , ω)φ ext G ′ (q , ω), (A. 10) and in view of the arbitrariness of φ ext , we can write in the matrix form χ(q , ω) = χ(q , ω) + χ(q , ω)Cχ(q , ω). (A.11) Finally, inverting Eq. (A.11), we arrive at Eq. (38). It may be useful to consider approximations to the above exact scheme. First, since G d ≫ 1 for G = 0, this is practically exact to keep in Eq. (A.2) the G = 0 term only. Further, if qa ≪ 1, where a is the effective width of the layer, Eq. (38) can be shown to yield a convenient relation between the corresponding 3D and 2D dielectric functions [17] 1 ε Q2D (q , ω) = 1 + 1 2 where ε 3D (q , ω; d) is the dielectric function of the fictitious periodic 3D system comprised of the layers separated by the distance d.