Extended frequency range of transverse-electric surface plasmon polaritons in graphene

The dispersion relation of surface plasmon polaritons in graphene that includes optical losses is often obtained for complex wave vectors while the frequencies are assumed to be real. This approach, however, is not suitable for describing the temporal dynamics of optical excitations and the spectral properties of graphene. Here, we propose an alternative approach that calculates the dispersion relation in the complex frequency and real wave vector space. This approach provides a clearer insight into the optical properties of a graphene layer and allows us to find the surface plasmon modes of a graphene sheet in the full frequency range, thus removing the earlier reported limitation (1.667<$\hbar\omega/\mu$<2) for the transverse-electric mode. We further develop a simple analytic approximation which accurately describes the dispersion of the surface plasmon polariton modes in graphene. Using this approximation, we show that transverse-electric surface plasmon polaritons propagate along the graphene sheet without losses even at finite temperature.


I. INTRODUCTION
Surface plasmon polaritons (SPPs) are collective excitations of charge density coupled to electromagnetic waves that can travel along a conductor-dielectric interface [1]. Interestingly, even an atomically thin conducting layer, such as a graphene sheet, can support SPPs [2]. In a graphite intercalated compound that contains multiple non-interacting two-dimensional graphene layers, SPPs were modelled theoretically and observed experimentally [3]. The first discovery of graphene [4] triggered more theoretical studies on SPPs in this material [5][6][7]. Since the first experimental observation of SPPs in graphene [8], more efforts have been made to use the SPPs for controlling the optical properties of graphene, such as electronically tuned extraordinary transmission [9], mode confinement by gap plasmons [10], and resonant absorption by an antidot array [11].
Although most studies of SPPs have focused on transverse-magnetic (TM) polarization, a single graphene layer, unlike normal-metal sheets, can also support transverse-electric (TE) SPP modes [12]. Technically, the existence of each type of SPP modes can be confirmed by solving a dispersion relation between the frequency and the wave vector following from Maxwell's equations and the conductivity model. For instance, the SPP dispersion for a dielectic/metal interface or a thin film of a Drude metal allows only TM modes which have zero magnetic field components normal to the interface and along the propagation direction. By analyzing the dispersion relation for complex wave numbers and real frequencies, Mikhailov et al. have shown [12] that TE SPP modes in graphene exist only in the range 1.667 < ω/µ < 2 at zero temperature, where ω and µ are, respectively, the light frequency and the chemical potential. Here, the lower limit corresponds to the zero * OhS2@cardiff.ac.uk of the imaginary part of the optical conductivity, while the upper limit is given by the minimum of the interband transition energy at zero temperature.
Importantly, the lower limit (≈1.667) for the normalized frequency ω/µ of the TE SPP mode can vary because the imaginary part of the optical conductivity of graphene may change depending on temperature and gate voltage. In addition, the range for the TE SPP mode frequency ω itself can be tuned by changing the chemical potential µ which in turn may be controlled by applying gate voltage or external magnetic field [13]. In the literature, the TE SPP mode solution has been studied for a graphene layer sandwiched between two dielectric media, and it was found that the range of the TE SPP mode frequencies can be modified by changing the permittivity contrast [14]. Furthermore, the range of TE SPP modes can be reversed when the surrounding material has a negative refractive index [15]. TE SPP modes in a graphene sheet placed on top of a nonlinear material substrate were also found to be limited to a similar frequency range [16].
In this work, we focus on the dispersion of SPP modes in a graphene layer with finite temperature and non-zero chemical potential and show that a complex-frequency analysis developed in this paper removes both the upper and the lower limits for the TE SPP mode in graphene.
While it is well-known how to calculate the dispersion relations of SPPs in graphene for given optical conductivities [12], this has been done assuming that any SPP mode has a real frequency but complex wave number q = q ′ + iq ′′ [12,16,17]. This corresponds to a continuous-wave excitation of a SPP which has a finite propagation length within the graphene layer of the order of 1/q ′′ . This picture is more suited for describing electromagnetic waves propagating in inhomogeneous waveguides not conserving the in-plane component of the wave number q. In contrast, uniform waveguides conserve q which can naturally be taken real. In our approach, assuming the excitation of the system is limited in time, the temporal evolution of SPP modes is described by a complex frequency ω = ω ′ + iω ′′ , with typically ω ′′ < 0 corresponding to a temporal decay due to radiative losses or absorption. The main advantage of this approach is that it provides a direct access to the optical spectra of the system where the imaginary parts of complex frequencies of isolated modes usually correspond to the half width at half maximum of the resonance peaks. A. Page et al. [18,19] have recently shown that the complex frequency approach can be used to describe optical gain (for ω ′′ > 0) in a non-equilibrium inverted graphene system where TM SPP modes were considered.
The complex-frequency approach has recently become widespread and broadly used in optics, owing to the useful concept of resonant states [20,21], also known in the literature as quasi-normal modes [22]. These are the eigen solutions of Maxwell's equations satisfying outgoing wave boundary conditions. They present a rigourous and powerful tool for analyzing optical spectra, such as scattering and transmission [23,24], with the real part of the complex frequency of the resonant state typically corresponding to the frequency position of a spectral line and the imaginary part to the half of its linewidth. Physically, resonant states form as a result of constructive interference of multiply reflected electromagnetic waves from the boundaries or inhomogeneities within optical systems. They have been studied in the literature both in finite optical systems, such as dielectric [25][26][27] and plasmonic nanoparticles [28][29][30], and in infinitely extended systems, such as planar waveguides [31] and photonic crystals [32][33][34]. Although in atomically thin films, such as a single graphene layer, these states do not normally exist, the formalism of complex-frequency modes can still be very useful, as we show in this paper.
The paper is organized as follows. Section II introduces the optical conductivity of a single graphene layer. Section III A describes the secular equations determining the dispersion of SPPs in TM and TE polarizations. Sections III B and III C present the main results of the paper, including both exact numerical and approximate analytical solutions of the secular equations in the complex frequency plane, their dependence on the propagation constant, temperature, and the chemical potential, and elimination of both the lower and the upper boundaries for the TE SPP mode frequencies. The temperature dependence of the threshold frequencies is discussed in Sec. III D. Section IV summarizes the results of the paper. Appendices A-E provide details on derivations of the optical conductivity of graphene and secular equations for the TM and TE SPP modes, and supply an additional material on our study of the TE mode near the lower threshold frequency and on the SPP dispersion at a finite damping.

II. OPTICAL CONDUCTIVITY OF GRAPHENE
The conductivity of graphene has an "interband" term in addition to the usual metallic Drude term, also referred to as "intraband". It is the interband term of the conductivity which gives rise to TE SPP modes found in a graphene layer in contrast to a Drude-metal layer where such modes do not exist. Following Refs. [12,17,35], we derive in Appendices A and B the two-dimensional (2D) optical conductivity of a homogeneous graphene sheet. In the long-wavelength limit, the expression for the 2D optical conductivity is given as a function of the light frequency ω by where is the Fermi-Dirac distribution, Ω = ω/µ is the normalized frequency, β = (k B T ) −1 is the inverse temperature, α = e 2 /( c) is the fine-structure constant, and Γ is a phenomenological damping. The first term in the curly bracket in Eq. (1) arises from intraband transitions, whereas the second one comes from interband transitions near the Dirac point in graphene dispersion. Importantly, Eq. (1) has an analytic dependence on Ω which can be continued into the complex Ω-plane without changing the formula. Note that for real Ω and Γ = 0, the integrand encounters a pole on the integration path. This requires that the diverging integral is split into a principal-value part and a half-pole contribution that can technically be achieved by keeping Γ positive infinitesimal in Eq. (1). For Im Ω = 0, this is no longer needed. However, the interband term is represented by a multi-valued function having a logarithmic nature. Therefore, for the analytic continuation, one has to choose the right Riemann sheet which provides the proper values of the integral at real Ω.

Figures 1 (a) and (b)
show the graphene conductivity as a function of normalized real frequency Ω = ω/µ and damping Γ at zero temperature. Increasing temperature smoothes out the Fermi distribution which is reflected by the smearing the step-like feature in the real part and the logarithmic divergence in the imaginary part of the interband conductivity at Ω = 2. The intraband Drude-like component of Eq. (1) has a zero-frequency pole which is moving away from the real axis as the damping increases.      Fig. 1 can also be understood as complex-frequency plots of the conductivity, treating Ω as the real and Γ (or its portion) as the imaginary part of frequency.
In the following we will be using, where it is convenient, the frequency Ω, wave number Q, and temperature (µβ) −1 normalized with respect to the chemical potential µ, treating the latter as a natural scaling parameter for, respectively, the frequency ω, wave number q, and temperature T . The chemical potential itself can be controlled e.g. by charge carrier concentration of graphene [36] and is also linked to the damping Γ. In the absence of voltage or in an undoped graphene layer, the valence band is fully filled and the conduction band is empty at zero temperature, implying that µ = 0 and Γ = 0. The concentration of free carriers -electrons in the conduction band and holes in the valence band -is zero, n 0 = 0. In a doped graphene layer or in the presence of voltage, the concentration of free carriers n 0 determines the chemical potential via the following equation [36,37]: where V is the electron Fermi velocity in graphene. At zero temperature, the above integral is equal to 1/2, which gives n 0 = µ 2 /(π 2 V 2 ). At a non-zero temperature, this integral depends on µ and T , in accordance with Eq. (2), thus introducing a temperature-dependent correction to the above expression for n 0 . The damping Γ contributing to Eq. (1) increases with density of impurities [17].

A. Secular equations for SPP modes
In this subsection, we present the secular equations determining in each polarization the dispersion relation between the SPP mode frequency ω and the propagation wave number q. We briefly describe their derivation for an infinitely thin graphene sheet. To gain physical insight of SPP modes in graphene, we discuss, at the end of this subsection, a comparison between SPP modes in graphene and in a Drude metal.
For a very thin planar conducting layer, the dispersion relation of SPP modes can be obtained in two ways. One may start deriving from Maxwell's equations and boundary conditions a secular equation for the SPP modes for a finite-thickness conducting material with a bulk conductivity. This secular equation can then be simplified in the limit of an infinitesimal film thickness. Alternatively, one may obtain the secular equation and SPP dispersion relation by assuming an infinitesimal layer with a surface conductivity σ(ω). The two ways lead to identical results. We follow the second approach and derive in Appendix C the secular equations for both polarizations, assuming the conducting layer is placed at z = 0 and is surrounded by vacuum, which is expressed by the permittivity in the entire space where δ(z) is the Dirac delta function. Using for brevity the units in which the speed of light in vacuum c = 1, the secular equations for SPP modes in both polarizations are given by where k(ω) = ω 2 − q 2 and q are, respectively, the normal and in-plane components of the light wave vector in vacuum, and ω is the complex light frequency. The dimensionless surface conductivity σ(ω) of a graphene layer to be used in the above equations is given by Eq. (1). Equations (5) and (6) determine the dispersion relations between the real in-plane wave number q and complex frequency ω. Note that the same equations were used in the literature [12,38] for finding SPP modes at a fixed real frequency of light, thus determining instead from Eqs. (5) and (6) complex propagation constants of SPPs. Now, before moving on to graphene, we briefly describe the solution of the secular equations (5) and (6) for a Drude metallic sheet. In this case, the conductivity would only consist of the first term in Eq. (1), so that σ(ω) ∝ i/ω, and Eq. (5) results in a TM mode having a square-root dispersion, ω ∝ √ q. We find in the next subsection a similar SPP mode for graphene in the frequency range dominated by the Drude conductivity. It can be seen that Eq. (6) for TE polarization has no solution with a non-zero real part of frequency for the Drude conductivity. We found, however, a TE mode of a Drude metallic sheet which has a purely imaginary frequency: Re ω = 0. This TE mode is also purely decaying in time: Im ω < 0. However, with the full graphene conductivity, there is also a propagating SPP mode in TE polarization [12], with Re ω = 0, which we discuss in more depth in the rest of this section.

B. SPP dispersion: Exact results and analytic approximations
In this subsection we present graphene SPP dispersion for both TM and TE modes which are found by solving the secular equations (5) and (6) numerically and comparing with a developed analytic approximation. We show in particular that the mode spectral linewidth (given by Im ω) is controlled by the temperature µβ contributing via the Fermi function Eq. (2). We analyze the earlier reported [12] limited frequency range 1.667 < ω/µ < 2 for the TE mode and remove both the upper and the lower boundaries for this mode. Figure 2 shows the dispersion of the TM mode, which includes both the complex ω plots and dependencies on q of the real and imaginary parts of ω and k, all shown at high (µβ = 10) and low (µβ = 10 5 ) temperatures. The TM dispersion lies well below the light line, ω = q, as can be seen from Figs. 2 (c) and (d) by comparing magnitudes of Re ω and q. In this frequency regime, the conductivity is dominated by intraband transitions within graphene band structure, described by the Drude-like term in Eq. (1). At low temperature, both Im ω and Re k are exponentially small, i.e. proportional to exp(−µβ). Note that Re ω asymptotically approaches the temperature dependent limit Ω TM up for large q. The magnitude of Im ω increases swiftly with q, with the rate determined by temperature, as more charge carrier vacancies become available below the Fermi level, compare Figs. 2 (a) and (b). The normal component of the wave number k shown in Figs. 2 (e) and (f) demonstrates that the SPP mode is localized in the z direction, due to Im k ≈ q > 0. At the same time, since Re k < 0 there is a propagation of light towards the conducting sheet, although it is much smaller than q.
pendix C 1: where σ ′ is the real part of the conductivity (taken at the mode frequency, ω ≈ √ ω 0 q). The approximation is based on the fact that in this frequency range, the conductivity is dominated by the Drude term, so that its imaginary part is given by σ ′′ ≈ ω 0 /ω whereas the real part can be treated as a small correction, i.e. |σ ′ | ≪ |σ ′′ |. Figure 3 presents a dispersion of the TE mode in graphene, again showing both the complex ω plots and dependencies ω(q) and k(q). The SPP mode in this polarization of light is unique to the graphene conductivity [12] due to the interband part of Eq. (1) and does not exist in a normal (e.g. Drude) metallic sheet. The real part of the dispersion curve, Re ω, lies close to the light   line ω = q, while the imaginary part, Im ω, is a few orders of magnitude smaller than the real part, as it is clear from Figs. 3 (c) and (d). Interestingly, both the real and imaginary parts of k, albeit being small compared to q, are now comparable to each other, in contrast to the TM SPP mode (compare Figs. 2 (f) and 3 (f)).
An approximate analytic solution was also obtained for the TE mode, by using the fact that away from the ω = 0 pole, the graphene conductivity is small, |σ(ω)| ≪ 1, as it is proportional to the fine-structure constant which is a small number, α ≪ 1. In this limit, the solution to Eq. (6) takes the form: where σ ′ and σ ′′ are both taken at ω = q, see Appendix C 2 for derivation. These complex frequency ω(q) and wave number k(q) are plotted as dashed curves in Figs For completeness, we also show in Appendix E the SPP mode dispersion both in TM and TE polarizations for the graphene conductivity at a non-zero damping of Γ = 0.05.

C. Removing the boundaries for the TE mode
The TE mode dispersion in graphene has been studied in Ref. [12] at zero temperature, with the mode frequency reaching but never exceeding the upper boundary at Ω TE up = 2. Our complex-frequency analysis allows us to eliminate this boundary, even for low temperatures. In fact, we see from Fig. 3 that the TE mode exists both below (Re Ω < Ω TE up ) and above (Re Ω > Ω TE up ) the threshold. Furthermore, we observe, by comparing Figs. 3 (c) and (d), that the temporal loss, which is given by −Im Ω, increases (decreases) with temperature below (above) the threshold. This can be understood simply as a smearing effect of the electronic distribution over the graphene band structure as temperature rises. Note that the threshold frequency is exactly twice the Fermi level E F of graphene, see the insets in Figs. 3 (c) and (d) which provide sketches of the electronic dispersion. At zero temperature, interband absorption only takes place above Ω TE up since no charge carriers occupy the electronic bands above the Fermi level, and so the losses are high for Re Ω > Ω TE up and zero for Re Ω < Ω TE up . When temperature is finite, some charge carriers occupy energy states just above the Fermi level, thus the interband absorption decreases for Re Ω > Ω TE up , thus reducing the losses. At the same time, vacancies of charge carriers are formed below the Fermi level at finite temperature and interband absorption can take place also for Re Ω < Ω TE up , so for this region the losses increase with temperature. This increase of losses is also reflected by the real part of the conductivity. In fact, it is clear from Fig. 1 (c) that the real part is smeared around Re Ω = Ω TE up as the temperature increases. No similar effects are observed for the TM mode. Figure 3 also demonstrates the earlier reported in the literature [12] lower threshold Ω TE low for the TE mode, which was observed at Ω TE low ≈ 1.667 at zero T . It is shown in Figure 3 by blue dots in the complex frequency plane and by red dots on the imaginary parts of the mode frequency ω and the normal component of the wave number k. Both imaginary parts change their sign at this threshold due to the change of sign of σ ′′ [see Eqs. (9) and (10) and also Fig. 1 (c)]. Physically, this threshold frequency corresponds to a condition that the intraband and interband electronic transitions are in balance. From a technical viewpoint, however, a positive imaginary part of the frequency, Im ω > 0, observed below the threshold, implies an exponential growth in time of the electric and magnetic fields (at any given point in space), whereas a negative imaginary part of the wave number, Re k < 0, also observed below the threshold, means an exponential growth of the field in space away from the graphene layer. While the latter is typical for radiative modes [21] and thus seems acceptable, the former usually corresponds to a gain [30], which is obviously not present in this system. One could therefore conclude that the TE mode does not exist below the threshold. However, the complex-frequency analysis allows us again to understand the properties of the SPP mode near the threshold and to eliminate this lower boundary.
To see that the TE mode exists both above and below the lower threshold at Ω TE low , we consider the spatial and temporal behaviour of the electric field together, having the following explicit form E(x, z; t) = E 0 e i(qx+kz−ωt) , where E 0 is a constant, and we have taken z > 0 for definiteness (see Appendix C for the analytic form of the fields). Now, separating the real and imaginary parts of the frequency, ω = ω ′ +iω ′′ , and of the normal component of the wave number, k = k ′ + ik ′′ , and using the fact that the in-plane component of the wave number q is real, we can separate the oscillating part of the field, E 0 e i(qx+k ′ z−ω ′ t) , from its amplitude, which is either exponentially decaying or exponentially growing in time and space. Using the fact that σ ′ > 0, which is equivalent to the positive imaginary part of the permittivity, we find from Eq. (10) that k ′ < 0, which physically corresponds to a plane wave propagating in vacuum towards the graphene sheet. The in-plane wave number q is however much larger, q ≫ |k ′ |, so that the direction of the electromagnetic wave is almost parallel to the sheet. At the same time, the z-coordinate of any point sitting on the wavefront and moving together with it can be approximately described as a function of time by z(t) = z(0) + k ′ t/q. Substituting this into Eq. (11), we find the amplitude of the field at the selected point on the wavefront to be where we have introduced a temporal decay rate Now, substituting here the dispersion for the TE mode, given by Eqs. (9) and (10), we find which implies that the TE mode has no losses in reality. In contrast, for the Drude-like TM mode, the decay rate found from Eq. (13) and the approximation Eqs. (7) and (8) is given by which implies an absorption, since the amplitude A(t) of the wave front decays with time in this case. However, for a Drude conductivity without damping, we obtain γ TM = 0, since the real part of the conductivity σ ′ vanishes, again implying that there are no losses in the system. The same is true for the graphene conductivity at zero temperature, since in this case σ ′ = 0 for Ω < 2. We would like to emphasize, however, that Eq. (14) is obtained for the TE mode at a non-zero temperature when the optical losses are present in the conductivity, since σ ′ > 0 (equivalent to Im ε > 0).
We see that the result for the TE mode, Eq. (14), is the same below and above the lower threshold at Ω TE low , and the TE mode demonstrates a fully physical behaviour from the energy conservation viewpoint on both sides of the threshold, as discussed above. Since below the threshold, the TE mode has, without any gain, an exponential growth with time, albeit at a very small rate compared to the mode frequency, this mode possesses a unique property, which has never been reported in the literature, to the best of our knowledge.

D. Temperature dependence of the threshold frequencies
Finally, we study the temperature dependence of the TE threshold frequency, Ω TE low . As discussed earlier in relation to Eqs. (9) and (10), this frequency threshold corresponds to a simultaneously change of sign of the imaginary part of the conductivity and the mode frequency. The red dashed line in Fig. 4 showing the threshold frequency Ω TE low as function of temperate is thus a solution of the equation Im σ(ω) = 0. Figure 4 clearly shows that the threshold frequency has a minimum, taking the value of Ω TE low = 1.6225 at 1/µβ = 0.0824, which is the result of a trade-off of the two terms in the graphene conductivity, corresponding to the intraband and interband transitions. More details on the properties of the TE mode near the lower threshold are provided in Appendix D, where an analytic equation determining its value at zero temperature (Ω TE low = 1.667) is derived. Interestingly, the TM modes has an upper threshold Ω TM up which coincides with Ω TE low at zero temperature but deviates for non-zero temperatures. It is also shown in Fig. 4, by a blue solid curve. This threshold, however, has a different physical meaning as it plays the role of an asymptote for the TM mode frequency in the limit of q → ∞. It can be seen that in this limit, Eq. (5) simplifies to σ(ω) = 0. Clearly, at zero temperature and frequencies below Ω TE up = 2, this coincides with the equation Im σ(ω) = 0 for Ω TE low , as can be seen also from Fig. 1 (a) demonstrating that Re σ(ω) = 0 in this frequency range.

IV. CONCLUSIONS
In summary, we have numerically calculated the complex-frequency dispersion of the surface plasmon polariton (SPP) modes in a homogeneous graphene layer both in transverse-magnetic (TM) and transverse-electric (TE) polarizations of light. We have further developed a simple analytic approximation which agrees well with the numerically exact solution of the secular equations for the modes in both polarizations.
We have shown that the TM SPP mode is determined by the Drude-like intraband part of the optical conductivity, demonstrating a square-root dispersion of the mode frequency with respect to the propagation wave number. In this polarization, the temporal decay of the electromagnetic field, which is given by the imaginary part of the mode frequency, monotonously increases with temperature. The TE SPP mode is in turn determined by both the intraband and the interband part, the latter being crucial for its existence. Unlike the TM mode, its dispersion is close to the light line, and the temporal decay demonstrates a nontrivial dependence on temperature and the propagation wave number. We have observed, in particular, that at finite temperature and chemical potential the TE mode exists above the upper threshold for the normalized frequency at Ω = 2, posed by an asymptotic behaviour of the dispersion at zero temperature. The temporal decay rate of the TE mode increases below the threshold, and decreases otherwise as the temperature increases. This is explained by considering occupation of electronic energy bands in graphene near the K-point at different temperatures.
We have also proven that the TE mode exists both above and below the lower threshold Ω = Ω TE low (taking the value of Ω TE low = 1.667 at zero temperature) and have studied its behavior near the threshold. This threshold is caused by a change of the sign of the imaginary part of the graphene conductivity which in turn causes a simultaneous change of the sign of the imaginary part of the SPP complex eigenfrequency and the normal component of the wave number. By investigating its spatial and temporal evolution, we have shown that the TE SPP in graphene presents a unique optical mode, as it can have below the threshold a positive imaginary part of the eigenfrequency without introducing gain into the material. Furthermore, we have demonstrated that in spite of the positive real part of the conductivity implying the positive imaginary part of the permittivity and hence an absorption, the TE SPP mode propagates along the graphene sheet without losses even at non-zero temperatures. This is correct at least up to second order in the conductivity.
Appendix A: Intraband conductivity of graphene Following [12,17,35], one can derive the twodimensional (2D) optical conductivity of a homogeneous graphene sheet by using the Kubo formula [6,[39][40][41], namely, by expanding the thermodynamic average of the current to first order in the amplitude of an external electric field. In the long wavelength limit, i.e. for a small light wave number compared to that of the electron, one may neglect effects of the spatial dispersion. The intraband conductivity is then given by [12] σ intra αβ (ω) = − where g s = 2 and g v = 2 are, respectively, the spin and valley degeneracies, S is the sample area, is the electron dispersion in graphene, with l = 1 and 2, is the Fermi-Dirac distribution function, and V is the electron Fermi velocity. Since we neglect the spatial dispersion, the conductivity tensor has a diagonal symmetric form, σ xx = σ yy = σ with σ xy = σ yx = 0. The expressions for the partial derivatives are given by which we substitute into (A1) to obtain after introducing g(E) = f (E) − f (−E) and performing integration over the angle in polar coordinates. The last integral can be evaluated analytically using integration by parts, which follows from the fact that for the Fermi function, Applying the limits of integration, we obtain and finally At zero temperature, µβ → ∞, the intraband conductivity simplifies to

Appendix B: Interband conductivity of graphene
The interband conductivity is derived in a similar way, i.e. again using the Kubo formula, which leads in the long wavelength limit to the following expression [12]: where 0 + is a positive infinitesimal, and corresponding to its eigenvalues Eq. (A2) and having the following explicit form: The matrix elements in Eq. (B1) then take the form Again using the symmetry of the conductivity tensor in the absence of the spatial dispersion, we obtain with the help of Eq. (A2) and after integration over the angle in polar coordinates Appendix C: SPP modes of a thin conducting sheet in vacuum To derive the secular equations (5) and (6) for the SPP modes in a graphene layer, let us consider a model of an infinitely thin sheet with 2D optical conductivity σ(ω), placed at z = 0. Choosing x as the propagation direction of light, so that the parallel components of the wave number are k x = p and k y = 0, Maxwell's equations are split into two blocks of first-order partial differential equations, separating TM and TE polarizations [34]: TE: where the speed of light in vacuum is taken c = 1 for brevity, ∂ z ≡ ∂/∂z, ε(ω; z) and µ(ω; z) are, respectively, the frequency and spatially dependent permittivity and permeability, E = A 0 (E x , E y , E z ) and H = A 0 (H x , H y , H z ) are, respectively, the electric and magnetic fields, ω is the light frequency, and A 0 (x, t) = e i(qx−ωt) is a common factor representing the temporal behaviour and the spatial dependence of the fields in the propagation direction. Clearly, Eqs. (C1) and (C2) can be obtained from each other by simultaneous swapping ε ↔ µ and E ↔ iH [34]. Let us further assume that where is a 2D susceptibility and σ(ω) is the 2D electrical conductivity of the graphene sheet, consisting of the intraband and interband components calculated in Appendices A and B. Note that the factor of 2 affecting the definition of σ is introduced in Eq. (C5) for convenience.

TM polarization
The set of equations (C1) for TM polarization then simplifies to For z = 0, we find the wave equation for H y and relations between the field components: To include the z = 0 point, we use Eq. (C4) and integrate Eqs. (C7), (C8), and (C9) over an infinitesimal interval including z = 0. Then we obtain where 0 + (0 − ) is a positive (negative) infinitesimal. From the above equations (C10) -(C15) we find after applying outgoing or incoming wave boundary conditions to the electro-magnetic field. Here A is a normalization constant, and k is the normal component of the wave number in vacuum satisfying the light dispersion (C20) Finally, using Eq. (C14), we obtain a secular equation for the SPP mode: which can be written more explicitly, using Eq. (C5), as identical to Eq. (5).
To obtain an approximate analytic solution to Eq. (C22), let us first consider the limiting case of very small frequencies when the conductivity is dominated by intraband transitions. This results in the standard SPP mode of an undamped Drude metal sheet. In fact, in this case where see Eq.

TE polarization
The secular equation for TE polarization is obtained in a similar way. Using Eqs. (C3) and (C4), the TE block given by Eq. (C2) can be written as For z = 0, we obtain and for z = 0, we integrate Eq. (C34) around this point, obtaining at the same time having both E y and H z continuous across z = 0. A solution satisfying outgoing or incoming wave boundary conditions then takes the form: and Eq. (C40) provides a secular equation for the SPP mode: which can be written more explicitly, using Eq. (C5), as identical to Eq. (6).
To obtain an approximate analytic solution of Eq. (C44), let us use the fact that |σ(ω)| ≪ 1 if the frequency is not too small. This is due to the fact that σ is proportional to the fine-structure constant α which is a small number. Combining Eqs. (C20) and (C44), we obtain Extracting the real and imaginary parts we then find Interestingly, k ′ < 0, since σ ′ (ω) > 0, at least for a real frequency -the same as in TM polarization. This implies that the light in the SPP mode propagates both along and towards the graphene sheet. At the same time, the amplitude of the wave exponentially decreases (increases) with distance from the sheet when σ ′′ (ω) < 0 (σ ′′ (ω) > 0). We also see that the sign of k ′′ and ω ′′ changes simultaneously at the lower threshold frequency as we discuss in detail in Sec. III B.

Appendix D: Lower threshold frequency for TE polarization
In this appendix, we discuss in more detail the condition for the lower threshold frequency Ω TE low of TE mode and derive an equation determining the threshold value of Ω TE low ≈ 1.667 at zero temperature. Let us first note that, in deriving Eq. (6), the electromagnetic field of the wave coupled to charge oscillations is proportional to where k = ω 2 − q 2 , q is the in-plane wave number in the direction of travel, and z is distance to the graphene sheet. For a bounded solution, we therefore assert that Im k > 0 .
In the opposite case, Im k < 0, the electromagnetic field would grow exponentially away from the graphene layer, and the threshold frequency is define as a value at Im k = 0. As an example for finite temperature µβ = 10, we see in Fig. 5 that the signs of both Im k and Im ω change simultaneously at the threshold frequency of Ω TE low ≈ 1.625. This is in agreement with the analytic approximation given by Eqs. (9) and (10).
The value of the threshold frequency depends on the temperature (and the chemical potential) as demonstrated by Fig. 4. It satisfies a general equation which can be easily obtained from Eq. (6) by using the fact that Im k = 0 at the threshold, and hence the mode (D5) The numerical solution of this equation gives Ω TE low ≈ 1.667. Figure 4 shows that Ω TE low approaches this value in the limit T → 0.  and TE polarizations. Figures 6 (a-d) show the dispersion of SPP TM mode for Γ = 0.05. Here we observe that, compared to zero damping, Im ω is mainly lifted up by Γ/2 throughout the range of q. Whereas Im k, shown as red curves in Figs. 6 (e) and (f), is unchanged compared to the case of Γ = 0 (Fig. 2), Re k, shown as blue curves in the same figures, changes its behaviour and its values significantly, increasing fast around q = 0 and then gradually decreasing with q.
In contrast, the SPP TE mode shown in Figs. 7 (a-d), demonstrates more changes in Im ω. In comparison with the dispersion for zero Γ (Fig. 3), the SPP TE mode has now smoother dependencies of Im ω and k. Im k, shown in Figs. 7 (e) and (f) starts from a negative value at q = 0. The TE SPP mode frequency still has a positive imaginary part in the range below the threshold frequency Ω TE low , similar to what we have seen for Γ = 0.