{'fraction_non_alphanumeric': 0.0389967998809258, 'fraction_numerical': 0.014363325146982213, 'mean_word_length': 4.1948453608247425, 'pattern_counts': {'":': 0, '<': 31, '<?xml version=': 0, '>': 15, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Electromagnetic properties of a double layer graphene system with electron-hole pairing

We study electromagnetic properties of a double layer graphene system in which electrons from one layer are coupled with holes from the other layer. The gauge invariant linear response functions are obtained. The frequency dependences of the transmission, reflection and absorption coefficients are computed. We predict a peak in the reflection and absorption at the frequency equals to the gap in the quasiparticle spectrum. It is shown that the electron-hole pairing results in an essential modification of the spectrum of surface TM plasmons. We find that the optical TM mode splits into a low frequency undamped branch and a high frequency damped branch. At zero temperature the lower branch disappears. It is established that the pairing does not influence the acoustic TM mode. It is also shown that the pairing opens the frequency window in the subgap range for the surface TE wave.

We study electromagnetic properties of a double layer graphene system in which electrons from one layer are coupled with holes from the other layer. The gauge invariant linear response functions are obtained. The frequency dependences of the transmission, reflection and absorption coefficients are computed. We predict a peak in the reflection and absorption at the frequency equals to the gap in the quasiparticle spectrum. It is shown that the electron-hole pairing results in an essential modification of the spectrum of surface TM plasmons. We find that the optical TM mode splits into a low frequency undamped branch and a high frequency damped branch. At zero temperature the lower branch disappears. It is established that the pairing does not influence the acoustic TM mode. It is also shown that the pairing opens the frequency window in the subgap range for the surface TE wave.

I. INTRODUCTION
Electron-hole pairing in a system of two conducting layers with oppositely charged carriers was predicted in 1,2 . Such a system below the Kosterlitz-Thouless transition temperature may support a superflow of electron-hole pairs. It should result in zero resistance under a flow of oppositely directed and equal in modulus electrical currents in the adjacent layers.
Zero counterflow resistance has not been observed yet, but a number of experiments confirm that the pairing occurs. The increase in the interlayer drag resistance at low temperatures was detected in a double quantum well in AlGaAs heterostructures 3,4 and in hybrid double layer systems comprising a single-layer (bilayer) graphene in close proximity to a quantum well created in GaAs 5 . The upturn of the drag resistivity witnesses for the approaching the paired state 6,7 . The electron-hole pairing was also predicted 8-10 for quantum Hall bilayers where the overall filling factor of the lowest Landau level of two layers is close to unity. The pairing was confirmed by transport experiments [11][12][13] where an exponential increase of the longitudinal counterflow conductivity and the vanishing of Hall resistance were discovered. Recent observation of a perfect longitudinal interlayer drag in the Corbino disk geometry 14 also witnesses for the electron-hole pairing.
After the experimental discovery of graphene the double layer graphene systems are considered as perspective ones for achieving the electron-hole pairing and counterflow superconductivity [15][16][17] . Potentially, such systems have a number of advantages. In graphene the electron and hole spectra coincide with each other at low energies and the condition of nesting of the electron and hole Fermi surfaces is fulfilled automatically. There is no gap between the electron and hole subbands and the density of carriers in the electron and hole layers can be easily controlled by external gates. At last, conducting electrons in graphene do not undergo a localization at low density of carriers.
The main obstacle in realizing the electron-hole pairing in graphene is the screening of the Coulomb interaction between electrons and holes 18,19 . Strictly spearing, the same obstacle emerges for the carriers with the quadratic spectrum, but in the latter case relatively high temperature of pairing can be achieved in the low density limit. In this limit electrons and holes bind in small size pairs (smaller than the average distance between the pairs) and the screening is suppressed. The Dirac carriers are not coupled in small size pairs, but they may form large size electron-hole pairs analogous to the Cooper pairs. Then the screening is also suppressed. Depending on the material parameters two situations are possible. If the bare interaction is weak, the screening will be strong and the critical temperature will be very low. Such a case was analyzed in 19 . On the contrary, if the bare interaction is strong, the pairing will be accompanied by an essential weakening of screening and it will occur at high temperature 20,21 . The latter possibility can be realized if the Coulomb interaction strength α ef f exceeds certain critical value α c and the distance between the graphene layers d related to the inverse Fermi wave number is small: dk F < ∼ 1. The interaction strength (also called the effective fine structure constant) is defined as α ef f = e 2 /hεv F , where v F is the Fermi velocity and ε is the dielectric constant of the surrounding media. The parameter α c grows up under increase in the number of Dirac components N (number of valleys times the number of spin components). For graphene (N = 4) the dynamical screening theory 20 yields α c ≃ 1.5. Thus a double layer graphene system in a vacuum (ε = 1 and α ef f = 2.19) can be in a paired state at relatively high temperatures (> 1 K). Other known two-dimensional Dirac crystals (silicene, germanene 22,23 , α-graphyne 24 ) have the same number of Dirac components, but smaller v F . Having the parameter α ef f approximately in two times larger than one for graphene, double layer silicene, germanene and α-graphyne look more promising for achieving high critical temperature. In particular, a double layer silicene(germanene) system embedded into a nanoporous or nanostructured matrix with ε < ∼ 2 25,26 can also demonstrate the electron-hole pairing at high temperature. Another perspective system is a topological insulator (TI) 27 . Two opposite surfaces of a TI plate serve as two adjacent conducting layers. The spectrum of the surface states of TI has an odd number of Dirac cones. For N = 1 the critical parameter α c ≃ 0.5 20 and the pairing at high temperature is possible in a TI plate embedded in a dielectric media with ε < ∼ 5 28 . One can also mention a possibility 29 of reducing the Fermi velocity in graphene and shifting the system parameters into the strongly coupled pairing regime by a periodic magnetic field.
Electron-hole pairing was also predicted for double layer graphene systems subjected to a strong uniform magnetic field directed perpendicular to the layers [30][31][32][33] . The important advantage of the quantum Hall state in Dirac spectrum systems is a large energy gap between the zeroth and the first Landau level. The screening in this case is considerable reduced even without pairing. Double layer system made of graphene sheets with the band gap (induced, for instance, by hydrogenation) is also a promising candidate for obtaining the electron-hole superfluidity 34 . The presence of the gap makes the formation of small size local pairs possible. One can expect that in that case the screening will be of a minor importance.
A direct manifestation of the electron-hole pairing would be an observation of zero counterflow resistance. But since in two-dimensional superfluid systems a flow causes unbinding of vortex pairs, a small voltage appears in any case. The presence of areas where the pairing is suppressed (in the case of system composed of puddles of the superconductive phase 35 ), can be another source of nonzero resistance. Therefore, it is desirable to have an independent indicator of pairing. For usual superconductors the Meissner effect can serve as such an indicator. Diamagnetic response of a double layer system with electron-hole pairing is very small 1,36 and it cannot be used for detection of the pairing. The pairing may reveals itself in a considerable enhancement of tunneling conductivity in the vicinity of the critical temperature 37 . This phenomenon interpreted as a fluctuational internal Josephson effect is the general one for electron-hole bilayers and it can be used as an indicator of pairing. The pairing reduces screening that can be seen by measuring the electric field of a test charge located near graphene layers. A strong increase of this field under lowering in temperature can be considered as a universal hallmark of the electron-hole pairing 36 , but it requires sensitive sensors with high spatial resolution.
In this paper we consider microwave methods of indirect observation of the pairing with reference to a double layer graphene system made of two gapless monolayer graphene sheets in zero magnetic field. We study the influence of the electron-hole pairing on the transmission and reflection characteristics of double layer graphene systems. We also analyze how the pairing changes the spectrum of surface plasmon-polaritons in double layer graphene systems. In Sec. II we develop the approach based on the generalized Nambu formalism. In difference with the original Nambu approach 38 we use the matrix Green's function of dimension 4×4. In Sec. III gauge invariant linear response functions are obtained. In Sec. IV an influence of the pairing on the transmission, reflection and absorption of electromagnetic waves in a double layer graphene system is analyzed. In Sec. V surface TM and TE waves in a double layer graphene system with the electron-hole pairing are studied.

II. EXTENDED NAMBU FORMALISM FOR A DOUBLE LAYER GRAPHENE
In the low-energy approximation the conducting electrons in a graphene layer are described by the matrix Hamiltonian where σ i are the Pauli matrices, σ 0 is the identity matrix, k is the wave vector counted from the Dirac point, µ is the chemical potential, α is the valley index, σ is the spin index, ψ + kασ and ψ kασ are the electron creation and annihilation operators that have the spinor structure The components of the pseudospinors (2) are the operators of creation and annihilation of electrons in the graphene sublattices A and B. A double layer system with the electron-hole pairing is described by the mean-field Hamiltonian presented in terms of 4-component spinors 17 where and 1 and 2 are the layer indexes. In the Hamiltonian (3)γ αβ are the 4 × 4 matrices defined through the direct productγ We imply the electron density in the layer 1 equals the hole density in the layer 2. It corresponds to µ 1 = −µ 2 = µ. Eq. (3) contains the order parameter of the electron-hole pairing ∆ k and the Hartree-Fock potential χ k . These quantities have to be found from the self-consistence equations.
In the general case 39 the order parameter is a 2 × 2 matrix∆ XX ′ , components of which describe the pairing of an electron in the sublattice X and a hole in the sublattice X ′ . This matrix can be expressed through the Pauli matriceŝ ∆ = µ=0,x,y,z ∆ µ σ µ . In this paper we consider the state, where only ∆ z is nonzero. It corresponds to the maximum energy gap in the quasiparticle spectrum and the minimum of the ground state energy 17 .
We would note the difference between the formalism used here and the approach 15,16,20,21,36 in which the standard Nambu notation 38 can be used 36 . In the latter case the order parameter is defined as an average ∆ k,λ = c + 1,k,λ c 2,k,−λ , where c + i,k,λ is the operator of electron creation in the state with the wave vector k in the Dirac subband λ in the layer i. The single layer Hamiltonian (1) written in terms of operators c i,k,λ has the scalar form. Then the Hamiltonian of a double layer system can be written in terms of two-component spinors.
In the extended Nambu formalism the Green's function is a 4 × 4 matrix Here and below we seth = 1. The Green's function (6)  It is convenient to present the Green's function (6) in the form where and ξ λ,k = λv F k − µ − χ k . The matrices A λ,k , B λ,k are expressed through theγ µν matrices: where θ k is the angle between the wave vector k and the x axes. The interaction part of the Hamiltonian written in terms of four-component spinors (4) reads where S is the area of the system, the notation : Ψ +γ ΨΨ +γ Ψ : means the normal ordering, V ± (q) = V S (q) ± V D (q), and V S (q), V D (q) are the Fourier-components of the intralayer and interlayer Coulomb interaction, respectively. We specify the case of a uniform dielectric environment that corresponds to the same intralayer interaction potential V S in both layers. The second sum in Eq. (11) compensates the mean-field terms in the Hamiltonian (3). The self-consistence condition requires the nullifying of the lowest order self-energy correction to the mean-field Green's function (6). The self-energy Σ(k, ω) is the 4 × 4 matrix and the self-consistence equation has the matrix form equivalent to 16 scalar equations. Most of them are satisfied by symmetry. Nontrivial ones are the self-consistence equation for the Hartree-Fock potential χ k and the equation for the order parameter ∆ k . The main effect of the Hartree-Fock potential is an unessential shift of the chemical potential that can be included in the definition of µ. The order parameter satisfies the equation where T is the temperature. We note that Eq. (12) does not contain the factor with cosine in difference with the self-consistence equation used in 15,16,20,21,36 . The origin of this difference is that the order parameter ∆ k introduced in the Hamiltonian (3) does not depend on λ. In the approach 15,16,20,21,36 the order parameter depends on λ. In that case the self-consistence equation has the form The consistence of Eq. (13) with the condition ∆ +1,k = ∆ −1,k requires zero contribution of the term with cosine into the integral in Eq. (13). Then Eq. (13) reduces to Eq. (12). The screening can be taken into account by replacing the bare interaction V D with the screened one. In the random phase approximation the Fourier-component of the screened interaction reads where Π ±,00 (k, ω) are the density-density response functions defined below.
In the static screening approximation the quantity V D (k) in Eq. (12) is replaced with V scr D (k, 0). In the dynamical screening approximation the self-consistence equation is modified to 20 where ω n = πT (2n + 1) are the odd Matsubara frequencies.

III. WARD IDENTITY AND GAUGE INVARIANCE OF THE RESPONSE FUNCTIONS
The interaction of the double layer graphene system with the electromagnetic field is described by the Hamiltonian whereĵ ±,0 are the charge density operators,ĵ ±,i -are the current density operators (i = x, y), A ±,i -are the vector potentials, A ±,0 = −cϕ ± , ϕ ± are the scalar potentials, and c is the light velocity. The index " + " (" − ") notates the sum (difference) of the corresponding quantities in the layers 1 and 2. The charge and current density operators are given by the equationĵ where the vertexes have the matrix form In the linear response approximation the currents j ±,µ are the linear functions of the potentials A ±,ν : The response functions Π ±,µν (q, ω) are obtained as the analytical continuation of the imaginary frequency currentcurrent correlators where N = 4 is the number of Dirac components,ĵ ±,µ (q, τ ) = e H mf τĵ ±,µ (q)e −H mf τ , β = 1/T , and T τ is the τ ordering operator.
Neglecting the interaction (11) one obtains the mean-field response functions It is well known 38,40 that the mean-field response functions in the Bardeen-Cooper-Schrieffer (BCS) theory are not gauge invariant. It is the result of that in the mean-field approximation the interaction is not taken into account properly. The interaction is included in the self-energy part while the renormalization of vertexes is neglected.
Any gauge transformation of the potentials A ±,µ should leave the currents (19) unchanged. This condition corresponds to the following equation for the response functions If Π ±,µν are not gauge invariant, the current (19) will not satisfy the continuity equation. Then considering an electromagnetic problem we would obtain one answer with the use of the boundary condition for the electric displacement field and another answer with the use of the boundary condition for the magnetic field. It would make the theory inconsistent. Therefore, for the problems we study in the next two sections the gauge invariance of the response functions is mandatory. Gauge invariance can be restored by dressing the vertexes. The response functions will be gauge invariant if the dressed vertexes Γ ± i satisfy the generalized Ward identity i=x,y This statement can be proven by the direct substitution of Eqs. (23) and (24) into Eq. (22). One can also prove that, as in the BCS theory 40 , the ladder approximation yields the vertexes that satisfy the Ward identity. For finding the vertex function in the ladder approximation one should solve a matrix integral equation. It is extremely cumbersome problem. Fortunately, the vertexes that satisfy the identity (24) can be found in a simpler way.
In the absence of interaction the Hartree-Fock potential χ k and the order parameter ∆ k are equal to zero and the Ward identity is fulfilled for the bare vertexes Γ ± µ = γ ± µ . If the Hartree-Fock potential is taken into account only as a shift of the chemical potential, the bare vertexes will satisfy the Ward identity in the normal state (∆ k = 0).
For the paired state we consider the gauge invariance problem in the constant gap (∆ k = ∆) and constant Hartree-Fock potential approximation. The constant Hartree-Fock potential is included below into the definition of µ. Then Eq. (24) reduces to One can see that the bare vertexes Γ + µ = γ + µ satisfy the Ward identity (25). Thus in the constant gap approximation the mean-field response functions Π +,µν (21) The functionsΓ − µ (q, iΩ) should have a pole at Ω → 0 and q → 0. This pole corresponds to the Anderson-Bogoliubov(AB) mode 41,42 . The AB mode is connected with fluctuations of the phase of the order parameter. In a three-dimensional superconductor the spectrum of this mode is Ω q = v F q/ √ 3. A genuine AB mode emerges only in neutral superfluids. In charged system the AB mode is coupled to the scalar potential and it renormalizes the electromagnetic response function. In the bilayer system the phase of the order parameter is coupled to the potential ϕ − and the AB mode renormalizes the response functions Π −,µν . One can show that in a two dimensional system the spectrum of the AB mode is modified as Ω q = v F q/ √ 2. Basing on the arguments given above, we seek for a solution of Eq. (27) The functions f µ (q, iΩ) are assumed to be regular at Ω → 0 and q → 0. Then, in the linear in q and Ω order we find f i = 2s 2 q i and f 0 = 2iΩ. Finally, we obtain the following renormalized vertexes The vertexes (29) satisfy the Ward identity and guarantee the gauge invariance of the response functions (23). We emphasize that Eqs. (29) are derived in the low frequency long wavelength limit.

IV. REFLECTION, TRANSMISSION AND ABSORPTION IN THE TERAHERTZ RANGE
The coefficient of transmission T for an electromagnetic wave going through an undoped graphene is almost frequency independent in a wide frequency range. The absorption is caused by the transition between filled and empty states in the Dirac cones. The coefficient of absorption is equal to A = πα, where α = e 2 /hc is the fine structure constant [43][44][45] . The reflectivity of graphene R is extremely small, proportional to α 2 . In the electron doped graphene at zero temperature the low energy states (counted from the Dirac point) are all filled. Therefore, there is no absorption 46 in the frequency range 1/τ ≪ ω < 2µ, where τ is the relaxation time. A double layer electron-hole graphene system with spatially separated layers should demonstrate the same behavior in the normal state. If the tunneling between the layers is negligible, an electron from one layer cannot transit to the empty state in the other layer and the layers contribute independently. In the paired state the quasiparticles do not belong to a given layer. At the same time, the quasiparticle spectrum has a gap 2∆. Therefore, one can expect the shifting of the absorption edge from ω = 2µ to ω = 2∆.
To obtain the frequency dependence of the transmission, reflection, and absorption coefficients we consider a ppolarized incident wave (the magnetic field of the wave is parallel to the graphene layers). At the frequency of the wave ω ∼ µ and for the interlayer distance d < k −1 F the ratio of the interlayer distance to the wavelength is d/λ < ∼ 10 −3 . Since d is much smaller than the wavelength one can consider a double layer graphene system as a single boundary. The electric current at that boundary is the sum of currents in two graphene layers j + . The parallel current conductivity tensorσ + relates the current j + with the electric field E + : j +,x = σ +,xx E +,x . For d/λ ≪ 1 the fields where the E x is the tangential component of the electric field at the boundary, and E +,x = 2E x .
Let the incident wave has the wave vector k = (k x , 0, k z ). Then the boundary condition for the magnetic field is Here the boundary at z = 0 is implied, and the indexes 1 and 2 stand for the z > 0 and z < 0 half-spaces. The boundary condition (30) together with Maxwell equations determine the following relations between the amplitudes of the electric field of the incident (i), reflected (r) and transmitted (t) waves where θ is the incident angle. The amplitudes of the magnetic components of the waves satisfy the same relations. In Eqs. (31) k pl = i x k x is the tangential component of the wave vector, and i x is the unit vector along the x axis. From the relations (31) one finds the transmission, reflection, and absorption coefficients for the normal incidence (θ = 0) where σ + (ω) = σ +,xx (0, ω) is the uniform high-frequency conductivity. Using the boundary condition for the normal component of the electric field (E 1,z − E 2,z = 4πρ + ) one obtains the same expressions for T , R and A. To get them one should take into account the continuity equation and the relation that comes from Eq. (22). The conductivity σ + (ω) is given by the response function Π +,xx The relation (33) allows to get σ + (ω) from the response function Π +,00 , as well: Using the expansion (7), computing the traces, summing over the imaginary frequencies, and doing the analytical continuation we get from (21) the following expressions for response functions Here and below the abbreviated notations E ≡ E λ,k , E ′ ≡ E λ ′ ,k+q , f ≡ (exp(E λ,k /T ) + 1) −1 , and f ′ ≡ (exp(E λ ′ ,k+q /T ) + 1) −1 are used. The factors F , P and L in Eqs. (36) and (37) have the form The factors (39) coincide with the coherence factors in the BCS theory 38,40 . We note that the four-component spinor formalism yields the same result for the response functions Π +,µν as the two-component spinor formalism 36 . But the important advantage of the formalism used in this paper is that it is compatible with the algorithms of checking and restoring of the gauge invariance developed in the BCS theory.
The integral in the expression for Π +,xx (q, ω) diverges. This divergence is unphysical one. It is connected with the linear approximation for the spectrum at large k. The same problem emerges for the monolayer graphene 47,48 . To escape this problem one should regularize the expression for Π +,xx (q, ω). Since the constant magnetic field cannot induce electrical currents in the normal system, the regularized response function Π r +,xx should be zero at ω = 0 and ∆ = 0 . This condition is fulfilled if one adds to the function (37) the compensating term that depends only on the cutoff wave vector k m : Here Π +,xx is given by Eq. The rule of integration of singularities in Eqs. (36), (37) is fixed by the standard substitution ω → ω + iη, where η = +0. For the computation we imply a finite η = 10 −3 µ considering it as a phenomenological scattering parameter.
We have computed σ + (ω) from the regularized response function Π r +,xx (qi x , ω) and from the response function Π +,00 (qi x , ω) and have got the same result. It confirms that the response functions obtained numerically are gauge invariant.
The real and imaginary parts of the conductivity σ + (ω) at the temperature T = 0.1µ and two different ∆ are presented in Fig. 1. For the comparison the conductivity in the normal state is also shown. The frequency dependence of the transmission, reflection and absorption coefficients for the system with pairing and for the normal system are given in Fig. 2. One can see that the pairing results in the appearance of a sharp peak in the reflection and absorption at ω = 2∆. This peak is accompanied with a deep minimum in the transmission. Below the gap (ω < 2∆) the system with pairing is completely transparent, while at ω > 2∆ it demonstrates significant absorption. In this range of frequencies the absorption coefficient decreases under increase in frequency up to ω = 2 ∆ 2 + µ 2 . At that frequency a step-like increase of the absorption coefficient occurs. It is connected with opening of an additional channel of absorption caused by the transition between λ = +1 and λ = −1 states (8). Normal systems demonstrate similar step-like features at ω = 2µ, but at much lower temperatures. At T = 0.1µ a step-like peculiarity is completely smeared out in the normal system (Fig. 2, left panel). In Fig. 2 the spectral characteristics are shown for ∆ = 0.5µ and ∆ = 0.2µ. One can see that lowering of ∆ results in a shift of the peak position. It the same time, the height and the width of the peaks remains practically unchanged. The computations show that the peaks caused by pairing are detectable already at ∆ = 0.01µ.
The gap in the excitation spectrum may emerge not only due to electron-hole pairing, but due to the interlayer tunneling, as well. Absorption and reflection coefficients are not sensitive to the gap origin and the tunneling can mimic the effect of pairing. To identify the origin of the gap one can take into account that the gap caused by pairing depends on temperature. For the systems with pairing we expect a red shift of the peak and their disappearance under increase in temperature.
Strong concentration mismatch of graphene layers may also modify the spectral characteristics of double layer systems. Such a mismatch results in an appearance of two step-like singularity at ω = 2µ 1 and ω = 2µ 2 . The amplitude and the shape of these singularities differ significantly from ones for the peaks caused by pairing. Besides, in the normal system the step-like singularities are smeared out at rather small temperatures. It allows to distinguish the effect of the mismatch and of the pairing.
Concluding the section we note that a similar problem, an impact of an excitonic gap in the high-frequency conductivity, was considered in 49 with reference to a single graphene layer where the interaction may cause an opening of a gap in the quasiparticle spectrum in zero or in a finite magnetic field 50 . It was shown in 49 that the low ω part of the interband contribution to the conductivity cuts off at 2µ or 2∆ whichever is the largest. For the bilayer system with the pairing we predict different behavior. The cutoff frequency is the minimum (not maximum) of 2µ and 2∆. In addition we find another step-like feature in the conductivity at the frequency ω = 2 µ 2 + ∆ 2 . Thus electron-hole pairing can be detected through an observation of strong reflection and absorption at the frequency ω = 2∆. Taking dk F ≈ 0.1, v F = 10 8 cm/s, and d = 10 nm we evaluate µ ≈ 6 meV and ν ∼ µ/h ≈ 1.6 THz. This estimate shows that for ∆ ∼ µ the features we have described can be observed in the terahertz spectral range.

V. SURFACE PLASMON-POLARITONS IN THE SYSTEM WITH PAIRING
Surface plasmon-polariton modes in graphene are now the subject of intensive study [51][52][53][54][55] . In particular, the interest to graphene is connected with the possibility of modifying the energy spectrum by external gates. The latter effect can be utilized in the plasmonic transformation optics 56 . There are two kinds of surface plasmon-polariton waves in graphene, the TM 57,58 and TE 59 ones. In the double layer system with the interlayer distance d smaller than the inverse Fermi wave number two TM modes, the symmetric (optical) and antisymmetric (acoustic) ones 60 , and one (symmetric) TE mode can propagate. The question we consider is how the pairing influences the spectrum and damping of these modes.

A. TM modes
We consider a double layer graphene system in a vacuum. The graphene layers 1 and 2 are located in the z = +d/2 and z = −d/2 planes and separated by a spacer with the dielectric constant ε d . The wave vector of the plasmon mode q is directed along the x axis. The electric field of the surface TM mode has the longitudinal (x) as well as the transverse (z) component, and the magnetic field has only the transverse (y) component. The electric and magnetic fields are exponentially decaying in both directions away from the graphene layers.
The solution of Maxwell equations with the boundary conditions (41) yields the dispersion equations for the sym- ) TM modes: where In a wide range of frequencies the wave vector of the surface TM mode is much larger than the wave vector of an electromagnetic wave with the same frequency in a free space (q ≫ ω/c). This condition is violated only for the symmetric mode at very small q. This range of q is not considered here. Implying also that q ≫ ε d ω/c we replace κ 1 and κ 2 in the dispersion equations (43), (44) with q and reduce Eqs. (43), (44) to the form where are the two-dimensional(2D) dielectric functions. The functions ε ± determine the screening of the scalar potentials ϕ ± = ϕ 1 ± ϕ 2 of the test charges q 1 and q 2 = ±q 1 located in the graphene layers 1 and 2 one above the other. The functions (46) can be presented in the form where One can see that Eq. (47) corresponds to the random phase approximation for the dielectric functions. Eq. (45) determines the dispersion of 2D plasmons. The difference between 2D plasmons and three-dimensional (3D) plasmons is the following. 3D plasmons are the longitudinal excitations of the electric field. Plasmons in 2D conductors in a 3D space are TM waves. The electric field of that wave has the longitudinal as well as the transverse component. The magnetic field is also nonzero but small in the limit q ≫ ω/c.
We specify the case dk F ≪ 1 and consider the range of wave vectors q < k F . Expanding V ± (q) in the small parameter qd and neglecting the higher order terms we obtain Let us first analyze the dispersion equation for the symmetric TM wave. According to Eqs. (45) and (49) this mode can propagate in the frequency range where Im [σ +,xx (qi x , ω)] > 0. The qualitative analysis can be done by replacing σ +,xx (qi x , ω) with σ + (ω). Then from Fig. 1 we see that in the system with the pairing at T = 0 the dispersion equation ε + (q, ω) = 0 may have two solutions, one is in the range ω < 2∆, and the other, in the range ω > 2∆. The real part of σ +,xx determines Landau damping. Since Re [σ + (ω)] is extremely small at ω < 2∆ one can expect that the low frequency solution corresponds to weakly damped plasmons. For ω > 2∆ the real part of conductivity is rather large and the high frequency solution will correspond to strongly damped plasmons.
The influence of pairing on Landau damping can be understood from the spectrum of electron-hole single particle excitations. In the system with pairing all states with negative energies (E = −E k,λ ) are filled at T = 0, and all states with positive energies (E = +E k,λ ) are empty. Therefore, the continuum of electron-hole excitations is determined by the inequality where the minimum is taken over all k, λ and λ ′ . It follows from this inequality that E e−h (q) ≥ 2∆. In the normal system the electron-hole excitation continuum starts from E e−h (q) = 0. Fig. 3 demonstrates the modification of the electron-hole excitation continuum under the pairing. One can see that in the paired state Landau damping for the modes with energies ω(q) < 2∆ is suppressed, while the modes with ω(q) > 2∆ may suffer from a strong damping. We compute the spectrum of the symmetric TM wave from the equation where ε + (q, ω) is given by Eq. (49) in which the dependence of σ +,xx on q is taken into account. The damping rate is evaluated as where ω(q) is the solution of Eq. (52). Since Eqs. (52), (53) are valid at small damping we do not consider the solutions of Eq. (52) with γ(q) > ∼ ω(q). The results of computation for ∆ = 0.5µ at T = 0.1µ and T = 0.2µ are presented in Fig. 4. One can see that in the state with the pairing the symmetric TM mode splits into two branches (the spectrum in the normal state is also shown in Fig. 4). In the long wavelength range the lower branch is a weakly damped one. There is a critical wave vector q c above which the solution of Eq. (52) that corresponds to the lower branch disappears. At q approaching q c the damping rate for the lower branch increases sharply. Under increase in temperature the frequency of this mode and the critical wave vector q c grows up . The lower branch is a thermally activated mode. At T = 0 this mode does not exist. It can be seen from the explicit expression for the response function Π +,00 (q, ω) (Eq. (36)). At T = 0 this response function is real and negative at ω < 2∆. Therefore Im [σ +,xx (q, ω)] < 0 and Eq. (52) has no solution in the frequency range ω < 2∆ at zero temperature. The behavior of the lower mode allows considering it as one connected with plasmon oscillations of the normal component decoupled from the superfluid one.
The frequency of the upper branch is restricted from below by the inequality ω > 2∆. The damping rate for this mode is much higher than for the lower branch. This branch exists in the wave vector range q > q c1 . At q close to q c1 the frequency of this mode approaches 2∆ and its damping rate increases sharply (the mode becomes overdamped).
In Fig. 5 the spectrum and the damping rate for the lower and upper branches at ∆ = 0.2µ and T = 0.1µ are shown. Comparing Fig. 5 and Fig. 4 we conclude that lowering of ∆ (at the same T ) results in an expansion of the wave vector range for two branches and in a decrease of the frequency and the damping rate of the upper branch.
The plasmon mode spectrum can also be extracted from the energy loss function This function determines the losses of energy of a pair of test charges oscillating with the frequency ω located in the adjacent graphene layers opposite to each other. The losses are connected with the excitation of plasmons at this frequency. A sharp peak in the loss function corresponds to a weakly damped plasmon mode. In Fig. 6 the dependence of L + (q, ω) on ω at fixed q is presented for T = 0 and T = 0.1µ. At T = 0.1µ this dependence contains a sharp peak that corresponds to the lower mode and a wide peak that corresponds to the upper mode. The positions of the peaks depend on ∆. At T = 0 the peak that corresponds to the lower mode disappears, while the upper mode peak remains unchanged. There is only one peak at ∆ = 0. It corresponds to the symmetric TM mode in the normal state. The behavior of the loss function demonstrates splitting of the TM wave into the lower and upper branches. The lower one is a weakly damped and thermally activated mode. The upper one is a strongly damped mode, and practically unsensitive to the temperature in the temperature range considered. Concluding this analysis we note that the splitting of the plasmon mode into a number of branches may also for the system is subjected by a strong magnetic field directed perpendicular to the graphene layers (without electron-hole pairing) 61 .
The spectrum of the antisymmetric TM mode computed from the equation Re[ε − (q, ω)] = 0 with ε − (q, ω) given by Eq. (50) is shown in Fig. 8. The parameters ε d = 4 and k F d = 0.2 are used for the computations. Fig. 8 shows that the antisymmetric TM mode remains the acoustic one in the state with the pairing and its velocity slightly reduces comparing to one in the normal state. In Fig. 9 the loss function L − (q, ω) is presented. Fig. 9 demonstrates that the antisymmetric mode is a weakly damped one in the normal as well as in the paired state. We remind that the renormalized vertex function (29) was obtained in the limit ω ≪ ∆ and its applicability at frequencies ω > ∼ ∆ is questionable.
The fact that the pairing does not influence the antisymmetric TM mode can be understood as follows. This mode corresponds to out-of-phase oscillations of electron densities in graphene layers. Out-of-phase oscillations in electron densities are equivalent to in-phase oscillations in electron and hole densities. It is quite natural that the electron-hole pairing does not suppress such oscillations.

+
where κ = q 2 − ω 2 /c 2 and k = ε d ω 2 /c 2 − q 2 . Eq. (59) corresponds to the symmetric TE wave, and Eq. (60), to the antisymmetric TE wave. At σ ± = 0 (no graphene layers) and the general d Eqs. (59), (60) have a number of solutions which correspond to the symmetric and antisymmetric waveguide modes. At qd < ∼ 1 only the lowest symmetric mode survives. The presence of graphene layers influences only the modes with the phase velocities close to the light velocity c. We restrict our analysis with frequencies comparable or smaller than the chemical potential µ. The mode with the phase velocity v ≈ c and the frequency ω < ∼ µ has the wave vector q < ∼ k F v F /c ≪ k F . Since the pairing occurs at k F d < 1, the strong inequality qd ≪ 1 is fulfilled. In this limit the real part of Eq. (59) is reduced to where α ≈ 1/137 is the fine structure constant,σ + (ω) is the conductivity normalized to e 2 /4h,ω = ω/µ and C = dk F (ε d − 1)/2πα ef f is the material parameter. In Eq. (61) we take into account that q is in 3 orders smaller than k F and replace the finite wave vector conductivity with the uniform conductivity σ + (ω). Eq. (61) determines the phase velocity v as the function of ω. Eq. (61) has a solution in the frequency range, where Im(σ + (ω)) < Cω. For typical parameters (dk F < 0.2, α ef f = 2.2, and ε d = 4) the constant C is very small (C < 0.05) and it can be neglected. In fact, the TE wave frequency range is determined by the same condition as one for the monolayer graphene (Im(σ + (ω)) < 0). In the normal state the TE mode exists only at ω > 1.67µ. One can see from Fig. 1b that the pairing opens the low frequency window ω min < ω < 2∆ for the TE mode. At temperature T → 0 the condition Im(σ + (ω)) < 0 is fulfilled for 0 < ω < 2∆ (see Fig. 10) and the lower edge ω min goes to zero.

VI. CONCLUSION
In conclusion, we have shown that the electron-hole pairing significantly changes spectral properties of double layer graphene systems in the terahertz range. The pairing causes the appearance of sharp high peaks in the absorption and reflection at the frequency ω = 2∆ and a rather large (much larger than in an undoped graphene) absorption at ω > 2∆.
The pairing influences essentially the surface symmetric TM mode. This mode splits into the lower and upper branches. The lower branch has the frequency ω < 2∆. It is practically undamped and appears only at nonzero temperatures in the long wavelength range. The spectrum of the lower branch is strongly temperature dependent. The upper branch frequency is in the range ω > 2∆, and this mode is strongly damped.
The influence of pairing on the antisymmetric TM mode is inessential. It is established that in the paired state a low frequency TE mode can propagate. In the normal state the frequency range for such a mode is restricted from below by the inequality ω > 1.67µ. In the paired system the additional frequency window ω min < ω < 2∆ opens, where ω min goes to zero at T → 0.