Optical absorption and conductivity in quasi-two-dimensional crystals from first principles: Application to graphene

This paper gives a theoretical formulation of the electromagnetic response of the quasi-two-dimensional (Q2D) crystals suitable for investigation of optical activity and polariton modes. The response to external electromagnetic field is described by current-current response tensor $\Pi_{\mu\nu}$ calculated by solving the Dyson equation in the random phase approximation (RPA), where current-current interaction is mediated by the photon propagator $D_{\mu\nu}$. The irreducible current-current response tensor $\Pi^0_{\mu\nu}$ is calculated from the {\em ab initio} Kohn-Sham (KS) orbitals. The accuracy of $\Pi^0_{\mu\nu}$ is tested in the long wavelength limit where it gives correct Drude dielectric function and conductivity. The theory is applied to the calculation of optical absorption and conductivity in pristine and doped single layer graphene and successfully compared with previous calculations and measurements.


I. INTRODUCTION
Understanding the interaction between light and electrons in a crystal has always been an attractive topic, and its extensive study led to the realization of many devices, such as lasers, semiconducting solar cells or LED diodes. More recently it led to new discoveries, such as subwavelength light transmission [1], waveguiding [2], hybrid solid state/organic solar cells [3], etc. There are many theoretical models which successfully deal with these phenomena, mostly based on solving Maxwell's equations at the boundaries of the crystals of different shapes [4] and different dielectric properties, calculated at different levels of approximations, e.g. within the Drude dielectric model [5,6], or from first principles and beyond the random phase approximation (RPA) [7].
However, what is still missing is a theoretical approach where the interaction between light and crystal electrons would be calculated fully microscopically, so that the electronic structure is calculated using ab initio methods (usually in the simplified tight-binding or subband models [8][9][10][11][12]), the polarizability of the system is described by the current-current response tensor (usually by the density-density response function [13]) and where the electron-electron interaction is mediated by photons (usually described by instantaneous Coulomb interaction). Inclusion of these effects could be crucial if one wanted to explore new optically active (radiative) modes or self-sustainable electromagnetic modes (polaritons) in crystals.
The aim of this paper is to give a theoretical formulation of the interaction between electromagnetic field and electronic excitations in quasi-two-dimensional (Q2D) crystals (consisting of one or few atomic layers), suitable for investigation of optically active electronic modes and polaritons. This formulation was partially developed in Ref. [14] where it was applied to investigate electromagnetic modes in the jellium metallic slab. Here this formulation is extended to the case where the ground state electronic structure is calculated from first principles. To test the theory we calculate the optical absorption and conductivity in the self standing graphene monolayer and compare our results with the recent experimental and theoretical works [15][16][17][18][19][20].
Optical properties of graphene have already been extensively investigated, both from the experimental and theoretical viewpoints. Apart from the above mentioned works, the information about optical activity of π or π+σ plasmons or single particle excitations was also extracted from electron energy loss (EELS) experiments and corresponding calculations [21][22][23][24][25][26][27]. In some cases the dispersion relations of 2D polaritons and conductivity in doped graphene are calculated using RPA density-density response function [26,[28][29][30] or using additionally a phenomenological relaxation-time approximation to account for the damping effects [6,31]. In Ref. [31,32] optical properties and conductivity in graphene are investigated at a high level of accuracy, beyond RPA, however the orbital and band structure are described within the tight binding approximation (TBA). On the other hand, in Ref. [15] optical properties of graphene are calculated from first principles including quasiparticle corrections and solving Bethe-Salpeter equation (so-called GW-BSE scheme) for the polarizability tensor. This GW-BSE scheme includes excitonic effects properly, resulting in a nice agreement of ultraviolet (UV) active π → π * peak with the experimental measurements, while our RPA theory underestimates this experimental value. Nevertheless, our theory provides very good agreement with the experiment in the infrared (IR) and visible regions and is capable to calculate the optically inactive (evanescent) modes such as surface plasmon polaritons (SPP).
It is worth mentioning, due to recent interest [33][34][35][36], that methodology presented in this paper can be adopted to obtain the charge-charge response function χ(Q, ω) and the dielectric function ε(Q, ω) of Q2D material by connecting the former with current-current correlation function Π µν (Q, ω) using the gauge invariance and the conservation of local charge density [31,37]. Numerically this can be very convenient for obtaining ε(Q, ω) within RPA response, because the Q 2 divergence of the bare Coulomb interaction is automatically cancelled due to Q prefactor in current vertex function and the special care for Q = 0 case [38] is not needed.
In Sec. II A we first present the general formulation of the problem, description of the system and the derivation of the Dyson equation for the screened current-current response tensor Π µν . In Sec. II B we formulate the Dyson equation for a specific geometry of the system in terms of Kohn-Sham (KS) electronic wave functions. In Sec. II C the expressions for optical absorption and conductivity are given in terms of the tensor Π µν . During ab initio calculation of nonlocal paramagnetic and local diamagnetic terms in Π µν certain numerical problems arise which are discussed in Sec. II D and resolved using an alternative expression for the current-current response tensor. In Sec. II E we prove that alternative expression of Π 0 µν in the long wavelength limit leads to the Drude dielectric function and conductivity. In Sec. III we apply the developed methodology to calculate the optical absorption spectrum and conductivity in doped and pristine graphene and compare it with recent experimental and theoretical results. Sec. III A gives details of the computational procedure. We use density functional theory (DFT) ground state calculation to get the crystal KS orbitals φ nK and band structure E nK . For the ground state calculation we use PWSCF method which means that our crystal should be 3D periodic. Here the superlattice consists of periodically repeated supercells containing several atomic layers (crystal slab). We also explain how to avoid the effect of interaction with the neighboring supercells. In the second stage of the calculation we solve the Dyson equation for current-current response tensor Π µν where irreducible current-current response tensor Π 0 µν enters and electromagnetic interaction is mediated by the photonic propagator D µν [14]. Here we restrict our calculations to RPA where irreducible current-current response tensor Π 0 µν can be obtained directly from the crystal KS states. A problem arises from the fact that we want to investigate optical properties of single Q2D crystal slab, while electronic structure is calculated for the entire 3D superlattice. This problem can not be solved simply by increasing the vertical separation between slabs because now each of them radiates electromagnetic field which spreads across the entire space and interaction between slabs is unavoidable. We solve this by cutting off the vertical range of the photon propagator D µν and allowing propagation only within one Q2D crystal slab. This procedure allows smaller vertical distances between adjacent slabs (enough to avoid the overlap between their charge densities) while still cancelling the interaction between them. Similar method is successfully utilized for calculation of EELS spectrum in graphene [26,27]. In Sec. III B we present and discuss in detail results for optical absorption spectra in pristine and doped graphene. In Sec. IV we present the conclusion.

II. FORMULATION OF THE PROBLEM
A. Derivation of the current-current response tensor In this section we will first derive the Dyson equation for the screened current-current response tensor in the Q2D crystal consisting of one or few atomic layers. We consider independent electrons which move in a local (DFT) crystal potential and interact with the electromagnetic field described by the vector potential operator A(r, t). Then the Hamiltonian of the system can be written as where describes non-interacting electrons in some local potential. Here c † nK is the creation operator of an electron in the Bloch state {n, K}, with the wave function φ nK (r) and energy E nK , where n is the band index and K is the electron wavevector in the plane parallel to crystal layers. H EM 0 is the Hamiltonian of free electromagnetic field. In the Φ = 0 gauge the interacting Hamiltonian consists of the paramagnetic part and the diamagnetic part Here the electron current operator is electron density operator is and electron field operators are The key quantity that will give the physical properties (e.g. optical properties) of the system of electrons interacting with the 'internal' electromagnetic field (electronelectron interaction mediated by photons) will be the current-current response tensor Π µν (r, r , t, t ), which in the 0-th order of the perturbation expansion over interaction V in contains two terms: Here the paramagnetic term is: and Ψ 0 e is the ground state of the Hamiltonian H e 0 . The diamagnetic term is Π dia µν (r, r , t, t ) = − e 2 mc n(r)δ(t − t )δ(r − r )δ µν (11) where n(r) = Ψ 0 e |ρ(r)| Ψ 0 e represents the ground state electron density. In the lowest order these two terms in the expansion of the irreducible current-current response tensor are represented diagrammatically in Fig. 1. Next step is to provide perturbation expansion for the tensor Π which now includes the higher order terms in the interaction V in . If we restrict our consideration within RPA the perturbation expansion of the current-current response tensor (9) becomes: where the symbol ⊗ denotes the convolution or integration over spatial and time variables (r, t) in addition to the matrix multiplication over indices µ = x, y, z. Also for clarity we omit to write spatial and time variables. From the expansion (12) it is obvious that the 'screened' current-current response tensor can be calculated by solving the Dyson equation where the only inputs are the non-interacting currentcurrent response tensor Π 0 given by (9) and the retarded free photon propagator given by where Ψ 0 EM is the photon vacuum (ground state of H EM 0 ), and the operator A 0 µ is defined as: The perturbation expansion of the current-current response tensor Π is diagrammatically presented in Fig. 2, where the green wavy line represents the external field A ext , that induces current fluctuations in the crystal. Blue wavy lines represent the propagator of electromagnetic field D 0 , that mediates electromagnetic interaction within the crystal. It should be mentioned here that this expansion actually goes beyond RPA. Namely, current-current response function (9) can be calculated by means of single-particle Green's functions. On the other hand electron-electron interaction is partially (depending on the DFT approximation used) included in the Bloch states used to construct the single-particle Green's function. Therefore the electron-electron interaction is included already in the a lowest order of the expansion (12) in the form of self-energy corrections to the irreducible polarizability Π 0 . The response tensor Π µν contains information about spectroscopic properties of the system (single-particle or collective electromagnetic modes in the system) or information about the response to the external electromagnetic field. Suppose that the crystal is exposed to the external (classical) electromagnetic field described by the vector potential A ext (r, t). In this case the total Hamiltonian (1) gets an additional term V ext which has the form analogous to (3)(4)(5) except that now A should be replaced by A ext . If V ext is a small perturbation it is  (16). Π 0 represents non-interacting current-current response tensor (9), Π interacting currentcurrent response tensor within RPA (13), and D 0 free photon propagator (14).
sufficient to keep only the term linear in A ext . In this case the current induced by the external field becomes as shown schematically in Fig. 2. The induced charge density fluctuations is similarly given by where we introduce the density-current response function which is (in the lowest order in expansion over V in ) given by Induced current and density fluctuations are connected by the continuity equation B. Calculation of the screened current-current response tensor In this calculation we shall first exploit the symmetry of the system which leads to the conservation of the Q vector parallel to the surface. Symmetry of the system also enables division to s and p polarizations. Suppose that the layered crystal ground state electronic density is in the perpendicular z direction restricted to the region −L/2 < z < L/2, as sketched in Figs. 3(b) and 3(c). Crystal periodicity is broken in the z direction but remains in xy plane, so it is appropriate to perform the Fourier transform of the Dyson equation (13) in the xy plane: where G = (G x , G y ) are 2D reciprocal lattice vectors. In (20) we also simultaneously performed the Fourier transform in ω-space. The Fourier transform of the free-photon propagator (14) is given by [14] D 0 where [39] D 0 (Q, ω, z, z ) = − 4πc Here the unit vectors are adapted to the geometry of the system such that e s = Q 0 × z and e p = c ω [−β sgn (z − z ) Q 0 + Qz] ( where Q 0 is the unit vector in the Q direction) represent directions of s(TE) and p(TM) polarized fields, respectively.
We see that the z integration in (20) is restricted exactly from −L/2 to L/2 which implies that the current fluctuation created in the region −L/2 < z 1 < L/2 can interact via photon propagator D 0 (Q, ω, z 1 , z 2 ) (even though the induced electromagnetic filed spreads over the all space) only with the current fluctuation in the region L/2 < z 2 < L/2. This restriction guarantees that Π contains information only about the electromagnetic modes characteristic for the electronic system limited to the region −L/2 < z < L/2 (e.g., Q2D systems).
The Dyson equation (20) can be additionally Fourier transformed in the z direction, so that it becomes a full matrix equation represents the Fourier transform of the current-current response tensor (9), and where the full Fourier transform of photon propagator can be obtained using (21) Using (6), (8) and (11) the Fourier transform of the paramagnetic contribution to the current-current response tensor becomes: where the current vertices are and where Here the 3D position vector is r = (ρ, z). Using (7) and (8) the Fourier transform of the diamagnetic contribution to the current-current response tensor becomes Here Ω = S ×L is the normalization volume, S is the normalization surface and f nK = (e (E nK E F )/kT + 1) −1 is the Fermi-Dirac distribution at temperature T . Integrations in (27) and (30) are performed over the normalization volume Ω. Plane wave expansion of the wave function has the form where the coefficients C nK are obtained by solving the KS equations self-consistently.

C. Optical absorption spectrum and conductivity tensor
If the layered system is in interaction with external electromagnetic field described by the vector potential A ext (r, t), the power at which the external electromagnetic energy is absorbed in the system can be obtained from the classical expression where in the Φ = 0 gauge the external electrical field can be calculated from After inserting the induced current (16) into (31) the absorption power becomes [40] Suppose now that the crystal layers plane lie parallel to the xy plane, as shown in Fig. 3(a), and the incident electromagnetic field is a plane wave of unit amplitude and polarization e: where the incident wave vector is q = (Q, β) and Q = (Q 0x , Q 0y ) is the wave vector parallel to the xy plane. The dispersion relation for electromagnetic waves in vacuum ω = |q|c leads to the relation β = ω 2 c 2 − Q 2 where Q = |Q|. This implies that for ω > Qc the perturbing field has radiative character with respect to the z-axis, as sketched in Fig. 3(b), and for ω < Qc it has evanescent character as sketched in Fig. 3(c). The radiative field can excite optically active modes, such as bright excitons, while the evanescent field is suitable for excitation of collective modes such as polaritons. The unit vector e represents the polarization of incident electromagnetic field. After combining Eqs. (32), (33) and (34), we do the Fourier transform of the current-current response tensor and the expression for the absorption power becomes Fourier transforming equations (16) and (32) in ω space and combining them we obtain (36) Using the fact that the induced current can also be calculated in terms of electrical conductivity tensor as we obtain the useful connection between the conductivity tensor and current-current response tensor Now we want to exploit our results for the tensor Π µν,GG (Q, ω). Fourier transforming it to the real space: and inserting in (35) the absorption power per unit area becomes where the spectral function is and the form factors are It is more convenient to deal with absorption coefficient A(Q, ω) = P (Q, ω)/|P| where the incident flux of electromagnetic energy (Poynting vector) is given by P = c 4π E×B. For unit amplitude incident electrical field (34) the incident flux is |P| = c 8π and the absorption coefficient is The Fourier transform of the conductivity tensor (38) can be obtained directly from the Fourier transform of current-current response tensor as The current which is induced by a homogeneous electric field E = e cos ωt (which corresponds to electromagnetic field (34)) which is incident normally (Q = 0) and for β 1/L can be, using (37) and (43), written as Here, as we are not interested in the current variation within the unit cell in the parallel xy direction, we retained only G = G = 0 components. If the field is directed along the x or y axis the current flow per unit thickness of the sample can be obtained by z integration in (44), when it becomes: which corresponds to the experimentally measurable conductivity.

D. Alternative expression for the current-current response tensor
For numerical reasons a straightforward calculation of the current-current response tensor (24) from Eqs. (26)(27)(28)(29)(30) can lead to non-physical results, so here we shall derive an alternative expression which avoids this problem. Namely, the expression for Π para (26) includes summation over all unoccupied bands, which is not the case for Π dia (29), so if the calculation is not performed with the same high precision the result for Π (24) might be erroneous.
Suppose for the moment that the system interacts only with the external electromagnetic field A ext (r, t) and that the interaction V in is neglected. Then the induced current and charge distributions are given by (16) and (17), respectively. After inserting (16) and (17) into the continuity equation (19) and performing the Fourier transform in (q, ω) space it becomes [37] where (q x , q y , q z ) = (Q + G , G z ) and ν = x, y, z. The Fourier transform of the current-current response tensor Π µν is given by (24), (26) and (29) and the Fourier transform of the charge-current response tensor Π 0ν , defined by (18), is given by where the current and charge vertices are defined by (27) and (30), respectively. In the static limit (ω = 0), when the current flow becomes stationary (also called the direct current or DC limit) the continuity equation (46) becomes It is important to note that Π 0 µν should be calculated so that condition (48) is satisfied very accurately, otherwise it can lead to some incorrect physical conclusions. For example, using the definition of the conductivity tensor (38), the DC conductivity becomes We see that if the condition (48) is not satisfied (49) could lead to the wrong conclusion about, e.g., the existence of superconducting state (σ → ∞). Also, in the normal metal state it could affect the Drude plasmon frequency. As already mentioned, the problem arises from the numerical calculation of the paramagnetic term (26) which includes summation over all unoccupied bands and will never be capable to cancel exactly the diamagnetic contribution (29) which includes only the summation over occupied states and can be calculated very accurately. This problem can be solved so that we calculate the paramagnetic term at some appropriate level of accuracy and then require the diamagnetic term to satisfy the continuity equation (48) After using (24), (26) and (50) the redefined currentcurrent response tensor becomes Now we shall demonstrate that the current-current response tensor given by (51) satisfies the continuity equation in the whole frequency range, i.e. it satisfies the equation (46). We start from the operator form of the continuity equation, which for the system of independent electrons and without interaction with the external field can be written as  (6) and (7) and the Hamiltonian H e 0 is given by (2). After Fourier transformation of equation (52), i.e.
dre −iqr [ρ(r, t), H e 0 ] = i ∇j(r, t) , and some commutation relations manipulation it becomes: By equating left and right sides in (53) we obtain a useful connection between the charge and current vertices After inserting charge vertices (54) into charge-current response function (47) and then into the left hand side of the continuity equation (46) it becomes exactly equal to its right hand side in which the current-current response tensor (51) is inserted. Therefore, current-current response tensor (51) satisfies the continuity equation for any ω. Important aspect of the new current-current response tensor (51) is appearance of the ω/ (E nK − E mK+Q ) prefactor which ensures that Π 0 µν (Q, ω) → 0 when ω → 0 and also naturally compensates the ω −1 divergence in the Kubo conductivity formula. E. The long wavelength limit, Q ≈ 0 We shall now analyze the long wavelength limit, Q ≈ 0, of the redefined current-current response tensor (51). First it is convenient to decompose the current-current response tensor into its intraband (n = m) and interband (n = m) contributions, For Q ≈ 0 and n = m we have that E nK − E mK+Q ≈ 0, and we can write the intraband contribution as [31,32], where we changed f nK → f (E nK ), and for simplicity we write Π 0 µν,GG (Q ≈ 0, ω) ≡ Π 0 µν,GG (ω). This intraband term leads to the Drude conductivity formula, as will be shown below. For n = m we have, Using some simple manipulation we can bring the expression (57) into the following form [41], One can easily see that the above expression for the interband term has a different behaviour in the static limit (ω = 0) than expression (57). The expressions (56), (57) and (58) will be used to calculate the adsorption coefficient (42) and conductivity (45) for optical wavevectors Q ∼ Q light . In addition, we shall verify that (51) leads to the correct dielectric function and conductivity in the three dimensional electron gas in the long wavelength (q ≈ 0) limit. For a polarisable system of arbitrary symmetry (in the linear response approximation) the electric field E and electric displacement D can be related as where ε is the non-local dielectric tensor and ⊗ is matrix multiplication and convolution in real space. After combining the definition (59) with the Maxwell and Dyson equations for electrical field E in the presence of a polarisable system we can obtain a general relationship between the dielectric tensor and the current-current response tensor [14,42] ε µν (r, r , ω) = δ(r − r )δ µν + 4πc ω 2 Π 0 µν (r, r , ω). (60) If we consider a 3D homogeneous electron gas the Fourier transform of (60) can be written as where Π 0 µν (q, ω) can be obtained from (51), G = G = 0, K → k becomes a 3D wave vector and Q → q becomes 3D momentum transfer. Also, after using the fact that in a homogeneous electron gas there is only one (n = m = 1) parabolic band (E k = 2 k 2 2m ) and the wave functions are plane waves (Φ k (r) = 1 √ Ω e ikr ) the current-current response tensor (51) becomes where we also used the definition of the current vertices (27) and (28). One can easily show that such a currentcurrent response tensor in the long wavelength limit becomes where n = 1 Ω k f (k) is the electron density. After inserting (63) into (61) we get the well known Drude dielectric function: where ω p = 4πne 2 /m is the bulk plasmon frequency. Similarly, combining (63) and the definition (38) leads to the Drude conductivity tensor [43] σ µν (q ≈ 0, ω) = i ne 2 m From the definition of the simple Drude DC conductivity σ = ne 2 m τ it is obvious that η plays the role of the inverse relaxation time, i.e. η = 1/τ .

III. RESULTS AND DISCUSSION
In order to illustrate the advantages of our theoretical approach we shall apply it to calculate optical properties of a free-standing single layer graphene.

A. Computational details
The first part of the calculation consists of determining the KS ground state of the single layer graphene and the corresponding wave functions and energies. For the unit cell constant we use the experimental value of a = 4.651 a.u. [44], 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 [45]. The core-electron interaction was approximated by the norm-conserving pseudopotentials [46], and the exchange correlation (XC) potential by the Perdew-Zunger local density approximation (LDA) [47]. To calculate the ground state electronic density we use 30 × 30 × 1 Monkhorst-Pack K-point mesh [48] of the first Brillouin zone (BZ) and for the plane-wave cut-off energy we choose 50 Ry. In order to achieve better resolution in the low energy and the static limit (ω → 0) the current response tensor (51) is evaluated from the wave functions φ nK (r) and energies E n (K) calculated for the 402 × 402 × 1 Monkhorst-Pack K-point mesh, and band summation is performed over 30 bands. In the calculation we use two kinds of damping parameters: η intra for transitions within the same band (n → n), and η inter for transitions between different bands (n → m). These two damping energies will be variable parameters. We take G = G = 0 in all the calculations because the crystal local field effects in the crystal layer plane (xy) are negligible in the optical limit (Q ≈ 0). However, broken symmetry in z direction results in big inhomogeneity of induced currents and fields in that direction. This requires inclusion of the crystal local field effects in z direction which we describe with 21 G z Fourier components. After solving the Dyson equation (23) we obtain the screened current-current response tensor Π which enters in the absorption coefficient (42). The conductivity tensor (43) can be, due to negligible screening for Q ≈ 0, calculated from the unscreened current-current response tensor (51).

B. Optical absorption
In this section we study the absorption of incident electromagnetic field (34) in the free-standing single layer graphene. We fix the parallel component of the incident wavevector Q = (Q 0x , Q 0y ) and change the incident frequency ω. Due to the relation ω = |q|c, for ω < Qc the perpendicular component of incident wavevector β is imaginary, and the incident field has evanescent character (in z direction), as sketched in Fig. 3(b). For ω ≥ Qc β is real and the incident field has radiative character in all three directions [ Fig. 3(c)]. In the latter case the incident wavevector q is inclined relative to the graphene surface by an angle θ given by as sketched in Fig. 3 (a). Let us discuss some specific cases. For example, for ω = Qc (on the light cone) the incident electromagnetic field is a plane wave which propagates parallel to the xy plane (θ = 0), s polarization is in xy plane, and p polarization is wholly in the perpendicular z direction. Also, if Q = 0, the electromagnetic field has radiative character in the whole frequency range (ω > 0), incidence is normal to the graphene surface (θ = π/2) and s and p polarizations become equal.
Retardation effects are most pronounced for very small wavevectors, so we shall divide our discussion of optical absorption into two parts: for small Q ≈ Q light , and large Q Q light , wavevectors.
1. Q ≈ Q light Figure 4(a) shows the optical absorption coefficient (42) of s(x) or p(y) polarised electromagnetic fields in pristine graphene for normal incidence Q = 0. Optical absorption onset appears already at ω = 0 which is due to the gapless dipole active π → π * interband transitions near K point of the BZ. In the infrared (IR) and visible regions (ω < 3 eV) absorption monotonically increases. The first absorption maximum, which appears in the ultraviolet (UV) region at ω = 4 eV, is a consequence of transitions between π and π * bands along the MM and MΓ directions of Q, as discussed in detail in [27]. This resonance absorbs about 12% of incident electromagnetic energy. In the far UV region ω > 6 eV the spectrum shows more structures, which are due to optically active transitions between σ and σ * bands, with the main peak at ω = 13.9 eV. This very strong excitation mode absorbs 30% of the incident electromagnetic energy.
Black solid line in the inset of Fig. 4(a) shows the details of optical absorption in IR, visible and UV regions. In the whole IR region the absorption is close to the universal value of πα = 2.293% (denoted by the horizontal dashed line), as predicted experimentally in Refs. [17,18]. In the far IR region (ω < 200 meV) the absorption begins to decrease faster until it reaches A(ω = 0) value which is about half of the universal value πα. However, the A(ω = 0) value strongly depends on the damping constants η intra and η inter used in the calculation. Blue dashed line is the theoretical result taken from Ref. [15], and red dotted line is the experimental result taken from Ref. [16]. We see that our absorption maximum is at a substantially lower energy (4.05 eV) then the 4.62 eV peak which appears in both, theoretical and experimental spectra. This is because the authors in Ref. [15] in their GW-BSE calculation included quasiparticle corrections of the DFT band structure and also the exitonic effects, i.e. electron-hole interaction. Here the spectrum is calculated within the RPA which includes screening effects, although this screening is very inefficient in the optical limit (Q ≈ 0), as discussed in Ref. [27]. We see quite a nice agreement of our results with the theoretical result from Ref. [15] in the whole IR and visible region, which will be the region of our main future interest. However, for ω < 0.5 eV both theoretical results start disagreeing with the experimental spectrum which suddenly decreases below the universal value. This is probably due to the weak doping which causes a shift of optical absorption onset from ω = 0 to ω = 2E F . Fig.  5 shows our theoretical results obtained using Eq. (57) (yellow line) and Eq. (58) (blue line) for the interband channel, which are then compared with the measured absorption spectrum for a graphene sample (red solid lines) taken from Ref. [18]. We see very good agreement in the frequency region 0.3 eV < ω < 1.2 eV (especially when using Eq. (58)), however we have slightly doped our graphene (E F = 0.1 eV) in order to achieve this re-   sult.
We want to emphasize here that Π zz (Q ≈ 0, ω) component of the full current-current response tensor is negligible in graphene [21], as seen in Fig. 4(b), so it is sufficient to use only x and y components for investigating electromagnetic response in graphene. Now we show the results for the two electromagnetic modes in doped graphene, appearing within the window constrained by the intraband and interband continua: the usual TM mode (2D plasmon-polariton or 2DPP), shown for Q ∼ 10 −7 and ω ≈ 0 in Figs. 6(a) and (b), and the TE mode [9,49,50], shown for Q ∼ 10 −4 and ω ∼ E F in Figs. 7(a) and (b). It can be seen that the 2DPP mode can be observed only if the incident EM wave is p polarised, while for the same (Q, ω) values and s polarised EM wave the only feature appearing is the Drude peak intersected with the light line (Qc). In the insets of Figs. 6(a) and (b) we compare the dispersion of the 2DPP with the √ Q 2D plasmon dispersion (when c → ∞). Here we see how coupling to the light waves influences the longitudinal plasmon mode and how its dispersion changes from √ Q to Qc for very small Q and ω when the electronic excitations of graphene are coupled to EM waves [51]. The effects of this coupling can also be seen in the second, TE, mode for ω ∼ E F . From Figs. 7(a) and (b) we see that the energy peak of this mode is always just below Qc, making it almost undistinguishable from the light line. Insets of these figures emphasise this even more.
In addition we present absorption intensities for ω ≈ 4 eV and p polarised EM wave, where the peak due to π → π * transitions appears [ Fig. 7(c)]. It is clear from the figure that there is no coupling between these transitions and EM waves, which is an additional proof in the pristine graphene plasmon debate, that there is no so-called π plasmon for Q ≈ 0 while for the large Q wavevectors the plasmon is formed [27,33]. where ∆Q = 0.0039 a.u.. The brown arrows indicate direction of increasing wavevectors. It should be noted that for n = 0 the whole frequency range in Fig. 8 corresponds to the radiative region, but for all nonzero wavevectors, n > 0 (e.g. for n = 1 Qc = 14.5 eV), the whole frequency range corresponds to the evanescent region.
Figs. 8(a), (b) and (c) show absorption spectra for s(x) polarized incident light for three different dopings E F = 0, 0.5 and 1 eV, respectively. For pristine graphene (E F = 0) far-IR absorption (ω < 200 meV) slowly decreases as Q increases and shows positively dispersive peaks, with the most intense one for n = 1 when it reaches almost 5% absorption. As Q increases the intensities of these peaks decrease. The black dots in the Fig. 8(a) inset show the dispersion relation ω s (Q) obtained by following the energies of these low energy peaks as functions of the wavevector Q. We see that ω s (Q) nicely follows the lower edge of the π → π * interband continuum, (v F Q), for pristine graphene, shown by black dashed line. The interband π → π * peak at ω ≈ 4 eV remains dispersionless.
For doped graphene the π → π intraband excitation channel is open and absorption spectra for s(x) polarized light get strong maxima in the far-IR region, as shown in Figs. 8(b) and (c). These absorption peaks are most intense for small Q's (n = 1, 2) and rapidly decrease as Q increases. The dispersion relations of this low energy peaks ω s (Q) are shown in the inset of Figs. 8(b) and (c). The upper edge of the intraband π → π electronhole continuum is also shown (black dashed lines) for comparison. Around the Dirac cone it can be approximated by v F Q, where for the Fermi velocity we took v F = √ 3ta/2 and where the hopping parameter of tightbinding model for graphene is t ≈ 2.7 eV (as in [26]). We can notice that both peaks are linearly dispersive, where for E F = 1 eV ω s (Q) follows exactly the v F Q line. Because of the Pauli exclusion principle, as the doping increases the lower threshold for interband π → π * transitions also increases. In other words, as these transitions are optically active the increased doping will open the optical gap and move the absorption onset towards higher energies. In Figs. 8(b) and (c) we can notice this wide optical gap in the whole IR and visible region (depending on the doping) with the optical absorption onset at 2E F , denoted by the vertical dashed lines.
The absorption of p(y) polarised incident light is shown in Figs. 8(d), (e) and (f). The parallel wave vector Q is chosen also to be in the y (ΓM) direction, so the response to this external electromagnetic field (for Q = 0) can be considered as longitudinal. It is evident from the results that this polarization gives generally more dispersive absorption spectrum in pristine graphene than the s(x) polarization. More specifically, we see that as Q increases the value of A(ω ≈ 0) rapidly decreases, and the π → π * peak (π plasmon) becomes blue shifted and less intense. In doped graphene the Dirac cones are partially filled and Q2D plasma is formed. As already mentioned, this plasma supports longitudinal self-sustaining oscillations called 2D plasmons [26], or in the electrodynamic limit, 2DPP. These 2DPP have evanescent character, i.e. electrical fields which they produce in z direction decay exponentially, as shown in Fig. 3(b). This implies that these modes exist in the ω < Qc region and cannot be excited (directly) by incident electromagnetic field which has fully radiative character. Figs. 8(e) and (f) show absorption spectra of p(y) polarised incident field for doped graphene with E F = 0.5 and 1 eV, respectively. For Q ≥ ∆Q the incident field has evanescent character in the shown frequency range and becomes capable to excite 2DPP, which gives rise to the strong peaks in the optical gap v F Q < ω < 2E F . Appearance of 2DPP causes strong screening which rapidly increases with Q. One consequence of such screening is that for p(y) polarization there is no intraband π → π absorption maximum in the far-IR region (as was the case for s(x) polarization). The intensity of intraband π → π transitions is strongly reduced by the 2DPP screening field. The 2DPP dispersion relations ω 2DPP (Q) shown by black dots in the inset of Figs. 8(e) and (f), are compared with the simple long wavelength approximation √ 2E F Q [26] (red solid line). The apparent disagreement between these two dispersion relations for ω > 0.8E F is because in our calculations we considered both intraband and interband transitions, while to get the simplified dispersion one only accounts for the intraband transitions. In the region below 0.8E F the agreement is better because there the interband contributions to 2DPP dispersion are negligible [26,31,32].

C. Optical conductivity
In this section we analyse optical conductivity obtained from (45) in the system exposed to the homogeneous electrical field directed in the µ = x or y direction and the wavevector is Q = 0. For Q = 0 the response is completely transverse, s(x) and p(y) polarizations are equivalent, the screening is inactive, so the screened Π and unscreened Π 0 current-current response tensors (51) become equal. As the crystal local field effects in the z di- rection are negligible, only the G z = G z = 0 component in Π 0 is nonzero. Also, because in the entire analysed frequency region ωL c → 0, the form factor in (41) becomes F (G z = 0) → 1. Under these conditions the absorption (42) becomes proportional to optical conductivity (45), i.e. A(ω) ∼ σ µµ (ω), and therefore the discussion of the optical conductivity corresponds to the discussion of absorption spectra which were thoroughly analysed in Figs. 4 and 5 of the previous section. However, here we shall emphasize the Drude limit (ω ≈ 0) and we focus on the case of pristine graphene only. Fig. 9 shows the interband, Eqs. (57) and (58), and intraband, Eq. (56), contributions to optical conductivity in pristine graphene. Our numerical method (being limited by finite K point sampling) is not especially appropriate to study different limits (e.g. T → 0, ω → 0 and η → 0) when approaching the ballistic minimum conductivity, as explained in Refs. [52][53][54], but it can be very suitable to explore the Drude regime. We use the finite temperature (T = 300 K) and energy damping (η > 0) and in this case the intraband channel in pristine graphene is open [ Fig. 9(c)]. Interband contributions to optical conductivity Eqs. (57) and (58) give two different behaviours for ω 0 as seen in Fig. 9(a) and (b). The first gives finite values up to ω = 0, while the second goes to 0 for ω = 0. In our case, where intraband channel presented in Fig. 9(c) is open, it is easily seen that Eq. (58) is more appropriate to use for the interband contribution to optical conductivity. Combining these two contributions one obtains good description of the experimental result from Ref. [18].

IV. CONCLUSION
In this paper we presented a microscopic theory of electromagnetic response in Q2D layered crystals, in terms of dynamically screened current-current response tensor. In this approach the explicit knowledge of the electromagnetic field propagator can give us information about different polariton modes, both radiative and nonradiative, their spectra and intensities, their coupling to external fields and to other excitations in the crystal (e.g. phonons). Specifically it is straightforward to evaluate optical properties of such crystals (absorption, reflection, transmission and conductivity) and also to calculate higher-order many body processes. The key physical quantity is the current-current response tensor, calculated from first principles which implies inclusion of the realistic crystal structure, wave functions and electronic band structure. In order to test the developed formulation we calculated optical absorption and conductivity in a single layer graphene and compared with recent experimental and theoretical results. The obtained results agree well with the measurements and experiments in IR and visible frequency regions, though the use of RPA is not capable to give correct excitonic effect. The theory is therefore suitable for the study of the layered nanostructures in the IR frequency range which is nowadays of great importance in plasmonics, which is a growing branch of theoretical but also of applied physics.