The Casimir Force between Two Graphene Sheets: 2D Fresnel Reflection Coefficients, Contributions of Different Polarizations, and the Role of Evanescent Waves

Abstract

We consider the Casimir pressure between two graphene sheets and contributions to it determined by evanescent and propagating waves with different polarizations. For this purpose, the derivation of the 2-dimensional (2D) Fresnel reflection coefficients on a graphene sheet is presented in terms of the transverse and longitudinal dielectric permittivities of graphene with due account of the spatial dispersion. The explicit expressions for both dielectric permittivities as the functions of the 2D wave vector, frequency, and temperature are written along the real frequency axis in the regions of propagating and evanescent waves and at the pure imaginary Matsubara frequencies using the polarization tensor of graphene. It is shown that in the application region of the Dirac model nearly the total value of the Casimir pressure between two graphene sheets is determined by the electromagnetic field with transverse magnetic (TM) polarization. By using the Lifshitz formula written along the real frequency axis, the contributions of the TM-polarized propagating and evanescent waves into the total pressure are determined. By confronting these results with the analogous results found for plates made of real metals, the way for bringing the Lifshitz theory using the realistic response functions in agreement with measurements of the Casimir force between metallic test bodies is pointed out.

1. Introduction

By now, graphene has assumed great importance in the field of fundamental physics and its numerous applications, where it plays a broad spectrum of roles [1,2]. The distinctive characteristic features of graphene, as opposed to ordinary bodies, are the 2-dimensional (2D) crystal structure of carbon atoms and massless quasiparticles described not by the Schrödinger equation, but by the Dirac equation, where the speed of light is replaced with the much smaller Fermi velocity. As a result, at energies below a few eV, the electrical and optical properties of graphene are well described by the relatively simple Dirac model [1,2,3,4,5,6]. This enables one to investigate the main features of graphene not by using some phenomenological approach, which is the usual practice in condensed matter physics, but on the solid basis of thermal quantum field theory and, more specifically, quantum electrodynamics at nonzero temperature.

The subject of this paper is the Casimir force [7], which acts between any two uncharged closely spaced material bodies owing to the zero-point and thermal fluctuations of the electromagnetic field. In his original publication [7], Casimir calculated the force acting between two ideal metal planes kept at zero temperature. At a later time, E. M. Lifshitz [8,9,10] developed the general theory expressing the Casimir force between two plates at any temperature via the reflection coefficients written in terms of the frequency-dependent dielectric permittivities of plate materials. In recent years, the Casimir force continues to grow in popularity due to the role it plays in quantum field theory, elementary particle physics, condensed matter, atomic physics, and even cosmology (see the monographs [11,12,13]).

Experiments measuring the Casimir force between metallic test bodies faced problems when comparing the measurement data with theoretical predictions of the Lifshitz theory. It turned out that if the low-frequency response of metals is described by the universally used dissipative Drude model, the obtained theoretical predictions are excluded by the measurement data. If, however, the low-frequency response is described by the dissipationless plasma model, which should not be applicable at low frequencies, the theory gives results in agreement with the experiment (see [12,14,15,16,17] for a review). Quite recently, it was shown [18] that the roots of the problem are not in the Drude model as a whole, but only in its possible incorrectness in the restricted area of transverse electric evanescent waves where it has no sufficient experimental confirmation.

The response functions of metals, including the Drude model, are of more or less phenomenological character. In this regard, of special interest is the Casimir effect in graphene systems, which has drawn the attention of many authors. At the early stages of investigation, the response of graphene to the electromagnetic field was also described by means of phenomenological methods based on the 2D Drude model, density functional theory, Boltzmann transport equation, random phase approximation, Kubo theory, hydrodynamic model, etc., and the obtained results were used to calculate the Casimir force in graphene systems [19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39]. In doing so, it was found that in the framework of the Dirac model the spatially nonlocal response of graphene at the pure imaginary Matsubara frequencies can be described by the polarization tensor in (2+1)-dimensional space-time and calculated precisely from the first principles of thermal quantum field theory [40,41]. These results were generalized to the entire plane of complex frequencies including the real frequency axis [42,43]. In such a manner, the reflection coefficients of electromagnetic fluctuations on a graphene sheet were expressed directly via the components of the polarization tensor.

The results of first-principles calculations of the Casimir force between two graphene sheets using the polarization tensor were compared [44] with those obtained using various phenomenological methods, and serious limitations of the latter were demonstrated. What is more, the measurement data of two experiments measuring the Casimir force in graphene systems were compared with the predictions of the Lifshitz theory using the reflection coefficients on graphene expressed via the polarization tensor and found to be in excellent agreement [45,46,47,48]. Specifically, the most precise measurements [47,48] confirmed the theoretical prediction of [24] that for graphene systems a big thermal effect in the Casimir force arises at much shorter separations than for metallic or dielectric bodies.

Thus, in the case of graphene, the Lifshitz theory does not suffer from a problem arising for metallic plates whose electromagnetic response was determined on partially phenomenological grounds (we recall that the experimental data for the complex index of refraction of metals are available only in the frequency region above some minimum frequency and are usually extrapolated by the Drude model to below this frequency [49]). One can conclude that graphene supplies us with some kind of road map on how to correctly describe the Casimir force between metallic plates. Because of this, it is important to compare both theoretical descriptions in parallel, including the form of reflection coefficients, the contributions of different polarizations of the electromagnetic field, and the propagating and evanescent waves.

In the current study, we underline that the reflection coefficients on a graphene sheet expressed via the polarization tensor are nothing more than the 2D Fresnel reflection coefficients expressed via the spatially nonlocal longitudinal and transverse dielectric permittivities. It is stressed that for a 2D graphene sheet, as opposed to the 3D Casimir configurations, the spatial dispersion can be taken into account exactly on a rigorous theoretical basis. Then, it is shown that in the application region of the Dirac model the Casimir force between two pristine graphene sheets is completely determined by the transverse magnetic polarization of the electromagnetic field. In doing so, at short separations up to hundreds of nanometers, both the propagating and evanescent waves make essential contributions to the Casimir force, whereas at larger separations the total force value is mostly determined by the contribution of evanescent waves. This is compared with the case of metallic plates where the evanescent waves play an important role in the problem of disagreement between the predictions of the Lifshitz theory using the Drude model and the measurement data.

The paper is organized as follows. In Section 2, we present the detailed derivation of the 2D Fresnel reflection coefficients on a graphene sheet in terms of the dielectric permittivities of graphene with an accurate account of the spatial dispersion. In Section 3, the explicit expressions for the transverse and longitudinal dielectric permittivities of graphene at any temperature are presented. Section 4 contains the results of the numerical computations of contributions to the Casimir pressure between two graphene sheets due to different polarizations of the propagating and evanescent waves. Section 5 provides a discussion of both similarities and distinctions in the Casimir pressures between metallic plates and graphene sheets. Finally, Section 6 contains our conclusions.

2. Fresnel Reflection Coefficients on a Two-Dimensional Sheet

It is known that the Casimir force between two parallel plates is expressed by the Lifshitz formula through the amplitude reflection coefficients of the electromagnetic waves. For the ordinary three-dimensional plates, these are the familiar Fresnel reflection coefficients written in terms of the frequency-dependent dielectric permittivity of the voluminous plate material. Graphene is a two-dimensional sheet of carbon atoms. Its dielectric permittivity is spatially nonlocal and essentially depends not only on the frequency, but on the wave vector and also on temperature. The expressions for the two-dimensional analogues of the Fresnel coefficients presented in terms of the dielectric permittivity of a 2D material are not quite known (see, e.g., [50,51,52], where the transverse magnetic coefficient [50,51] and both reflection coefficients [52] were expressed in terms of the 2D conductivity with no account of spatial dispersion, or [13,31] where they are presented with no detailed derivation).

Below, we demonstrate in detail that the reflection coefficients on a 2D sheet are obtainable in close analogy to the standard 3D Fresnel reflection coefficients, but with due account of the spatial dispersion.

Let the graphene sheet be in the plane $z=0$, where the z-axis is directed downward in the plane of Figure 1 and the y-axis is directed upward perpendicular to it.

There are empty half spaces I and II on both sides of the graphene sheet. The fluctuating electromagnetic field induces some surface charge density, ${\varrho }^{2\mathrm{D}}(\mathit{\rho },t)$, and current density, ${\mathit{j}}^{2\mathrm{D}}(\mathit{\rho },t)$, on the sheet where $\mathit{\rho }=(x,y)$ and t denotes the time. Then, the Maxwell equations in the 3D space take the form

$$\begin{array}{ccc}& & \nabla \mathit{D}(\mathit{r},t)=4\pi {\varrho }^{3\mathrm{D}}(\mathit{r},t),\hfill \\ & & \nabla \mathit{B}(\mathit{r},t)=0,\hfill \\ & & \nabla \times \mathit{E}(\mathit{r},t)+\frac{1}{c}\frac{\partial \mathit{B}(\mathit{r},t)}{\partial t}=0,\hfill \\ & & \nabla \times \mathit{H}(\mathit{r},t)-\frac{1}{c}\frac{\partial \mathit{D}(\mathit{r},t)}{\partial t}=\frac{4\pi }{c}{\mathit{j}}^{3\mathrm{D}}(\mathit{r},t),\hfill \end{array}$$

(1)

where $\mathit{r}=(x,y,z)=(\mathit{\rho },z)$, c denotes the speed of light, D is the electric displacement, B is the magnetic induction, and E and H are the electric and magnetic fields, respectively. The 3D charge and current densities in Equation (1) are given by [13,50]

$${\varrho }^{3\mathrm{D}}(\mathit{r},t)={\varrho }^{2\mathrm{D}}(\mathit{\rho },t)\delta \left(z\right),{\mathit{j}}^{3\mathrm{D}}(\mathit{r},t)={\mathit{j}}^{2\mathrm{D}}(\mathit{\rho },t)\delta \left(z\right).$$

(2)

Note that we use the Gaussian units in Equation (1) and below. In these units, ${\mathit{j}}^{3\mathrm{D}}$ has the dimension of ${\mathrm{g}}^{1/2}{\mathrm{cm}}^{-1/2}{\mathrm{s}}^{-2}$, whereas the dimension of ${\mathit{j}}^{2\mathrm{D}}$ is ${\mathrm{g}}^{1/2}{\mathrm{cm}}^{1/2}{\mathrm{s}}^{-2}$.

The standard electrodynamic boundary conditions on the plane $z=0$ are given by

$$\begin{array}{ccc}& & [{\mathit{D}}_{\mathrm{II}}(\mathit{\rho },0,t)-{\mathit{D}}_{\mathrm{I}}(\mathit{\rho },0,t)]·\mathit{n}=4\pi {\rho }^{2\mathrm{D}}(\mathit{\rho },t),\hfill \\ & & [{\mathit{B}}_{\mathrm{II}}(\mathit{\rho },0,t)-{\mathit{B}}_{\mathrm{I}}(\mathit{\rho },0,t)]·\mathit{n}=0,\hfill \\ & & [{\mathit{E}}_{\mathrm{II}}(\mathit{\rho },0,t)-{\mathit{E}}_{\mathrm{I}}(\mathit{\rho },0,t)]\times \mathit{n}=0,\hfill \\ & & [{\mathit{H}}_{\mathrm{II}}(\mathit{\rho },0,t)-{\mathit{H}}_{\mathrm{I}}(\mathit{\rho },0,t)]\times \mathit{n}=-\frac{4\pi }{c}{\mathit{j}}^{2\mathrm{D}}(\mathit{\rho },t),\hfill \end{array}$$

(3)

where $\mathit{n}=(0,0,1)$ is the unit vector directed along the z-axis (see Figure 1).

Below we assume that all fields have the form of monochromatic plane waves, e.g.,

$$\mathit{E}(\mathit{r},t)={\mathit{E}}^0{e}^{i(\mathit{k}\mathit{r}-\omega t)},\mathit{H}(\mathit{r},t)={\mathit{H}}^0{e}^{i(\mathit{k}\mathit{r}-\omega t)},\mathit{B}(\mathit{r},t)={\mathit{B}}^0{e}^{i(\mathit{k}\mathit{r}-\omega t)}.$$

(4)

Here, ${\mathit{E}}^0$, ${\mathit{H}}^0$, and ${\mathit{B}}^0$ are the amplitudes; $\mathit{k}=(k_x,k_y,k_z)\equiv (\mathit{q},k_z)$ is the 3D wave vector, and $\omega$ is the wave frequency.

For a derivation of the Fresnel reflection coefficients on a 2D sheet, it is suffice to restrict our consideration to the third line of the Maxwell equations (1) and the third and fours lines in the boundary conditions (3).

Substituting Equation (4) into the third line of Equations (1) and (3), it is readily seen that in both regions I and II

$$\mathit{k}\times {\mathit{E}}^0-\frac{\omega }{c}{\mathit{B}}^0=0$$

(5)

and

$$({\mathit{E}}_{\mathrm{II}}^0-{\mathit{E}}_{\mathrm{I}}^0)\times \mathit{n}=0,$$

(6)

where ${\mathit{E}}_{\mathrm{I}}^0$ and ${\mathit{E}}_{\mathrm{II}}^0$ are the field amplitudes in the regions I and II, respectively.

Now we look at the fourth line in the boundary conditions (3). Taking into account that the graphene sheet is a spatially nonlocal material, the 2D current density in the fourth line of Equation (3) takes the form

$${\mathit{j}}^{2\mathrm{D}}(\mathit{\rho },t)=\underset{-\infty }{\overset{t}{\int }}d{t}^{\prime }\int d^2{\mathit{\rho }}^{\prime }{\sigma }^{2\mathrm{D}}(\mathit{\rho }-{\mathit{\rho }}^{\prime },t-t^{\prime }){\mathit{E}}_{\mathrm{lat}}({\mathit{\rho }}^{\prime },t^{\prime }).$$

(7)

Here, ${\sigma }^{2\mathrm{D}}(\mathit{\rho },t)$ is the 2D conductivity of a graphene sheet (it has the dimension cm/s) and ${\mathit{E}}_{\mathrm{lat}}$ is the projection of the electric field on the plane of graphene calculated at $z=0$:

$${\mathit{E}}_{\mathrm{lat}}(\mathit{\rho },t)=\mathit{E}(\mathit{\rho },0,t)-\mathit{n}\left(\mathit{E}(\mathit{\rho },0,n)·\mathit{n}\right)=\mathit{n}\times [\mathit{E}(\mathit{\rho },0,t)\times \mathit{n}].$$

(8)

Substituting Equations (4) and (7) into the fourth line of Equation (3), one obtains

$$({\mathit{H}}_{\mathrm{II}}^0-{\mathit{H}}_{\mathrm{I}}^0)\times \mathit{n}=-\frac{4\pi }{c}{\sigma }^{2\mathrm{D}}(\mathit{q},\omega ){\mathit{E}}_{\mathrm{lat}}^0,$$

(9)

where ${\sigma }^{2\mathrm{D}}(\mathit{q},\omega )$ is the Fourier image of ${\sigma }^{2\mathrm{D}}(\mathit{\rho },t)$ in the 2D space and time, q is the 2D wave vector, and ${\mathit{E}}_{\mathrm{lat}}^0$ is the amplitude of the quantity (8):

$$\begin{array}{ccc}& & {\mathit{E}}_{\mathrm{lat}}(\mathit{\rho },t)={\mathit{E}}_{\mathrm{lat}}^0{e}^{i(\mathit{q}\mathit{\rho }-\omega t)},\hfill \\ & & {\mathit{E}}_{\mathrm{lat}}^0=\mathit{n}\times [{\mathit{E}}^0\times \mathit{n}].\hfill \end{array}$$

(10)

Note that by introducing ${\sigma }^{2\mathrm{D}}(\mathit{q},\omega )$ we have used the translational invariance in the plane of a graphene sheet. In the standard Casimir problems, where the plates are made of 3D materials separated by a gap, there is no translational invariance in the 3D space and it is impossible to rigorously introduce the conductivity ${\sigma }^{3\mathrm{D}}(\mathit{k},\omega )$ (and the dielectric permittivity) depending on the 3D vector k. Because of this, for taking into account the effects of spatial dispersion, it is necessary to use some approximations, such as the suggestion of specular reflection [53,54].

We recall also that the spatially dispersive materials are characterized by the two independent conductivities, in our case, ${\sigma }^{2\mathrm{D},\mathrm{L}}(\mathit{q},\omega )$ and ${\sigma }^{2\mathrm{D},\mathrm{Tr}}(\mathit{q},\omega )$, depending on whether ${\mathit{E}}_{\mathrm{lat}}^0$ in Equation (9) is parallel or perpendicular to the wave vector q, respectively [55,56]. These conductivities are called the longitudinal and transverse ones.

We are coming now to the derivation of the amplitude reflection coefficients on a graphene sheet for two independent polarizations of the electromagnetic field using Equations (5), (6) and (9).

Let us start with the case of transverse electric polarization when the amplitudes of the electric field of the incident, ${\mathit{E}}_0^0$, transmitted, ${\mathit{E}}_2^0$, and reflected, ${\mathit{E}}_1^0$, waves are perpendicular to the plane of incidence $(x,z)$ and directed along the positive direction of the y-axis (see Figure 2). The corresponding wave vectors are ${\mathit{k}}_0$, ${\mathit{k}}_2$, and ${\mathit{k}}_1$, and the amplitudes of the magnetic field, which lie in the plane of incidence, are ${\mathit{H}}_0^0$, ${\mathit{H}}_2^0$, and ${\mathit{H}}_1^0$.

Taking into account that the 2D sheet spaced in the plane $(x,y)$ or, equivalently, $z=0$ is spatially homogeneous, one finds $k_{0x}=k_{1x}=k_{2x}$. Considering also that $k_0^2=k_1^2=k_2^2={\omega }^2/c^2$ because the space outside of a graphene sheet is empty, one obtains

$$sin{\theta }_0=\frac{k_{0x}}{k_0}=sin{\theta }_1=\frac{k_{1x}}{k_1}=sin{\theta }_2=\frac{k_{2x}}{k_2},$$

(11)

i.e., in our case, all the three angles are equal.

According to Figure 2,

$${\mathit{E}}_{\mathrm{I}}^0={\mathit{E}}_0^0+{\mathit{E}}_1^0,{\mathit{E}}_{\mathrm{II}}^0={\mathit{E}}_2^0,$$

(12)

where ${\mathit{E}}_0^0=(0,E_{0y}^0,0)$, ${\mathit{E}}_1^0=(0,E_{1y}^0,0)$, and ${\mathit{E}}_2^0=(0,E_{2y}^0,0)$.

Taking this into account, the boundary condition (6) reduces to

$$E_{0y}^0+E_{1y}^0=E_{2y}^0.$$

(13)

The boundary condition (9), where ${\mathit{E}}_{\mathrm{lat}}^0$ is defined in Equation (10), is more complicated. In view of Equation (6), both ${\mathit{E}}_{\mathrm{I}}^0$ and ${\mathit{E}}_{\mathrm{II}}^0$ can be substituted into Equation (10) in place of ${\mathit{E}}^0$. We choose ${\mathit{E}}_{\mathrm{II}}^0$ for the sake of brevity. Then the condition (9) takes the form

$$({\mathit{H}}_{\mathrm{II}}^0-{\mathit{H}}_{\mathrm{I}}^0)\times \mathit{n}=-\frac{4\pi }{c}{\sigma }^{2\mathrm{D},\mathrm{Tr}}(\mathit{q},\omega )[\mathit{n}\times [{\mathit{E}}_{\mathrm{II}}^0\times \mathit{n}]].$$

(14)

Here, we took into account that the electric field is perpendicular to q.

From Figure 2, it follows that

$${\mathit{H}}_{\mathrm{I}}^0={\mathit{H}}_0^0+{\mathit{H}}_1^0,{\mathit{H}}_{\mathrm{II}}^0={\mathit{H}}_2^0,$$

(15)

where ${\mathit{H}}_0^0=(H_{0x}^0,0,H_{0z}^0)$, ${\mathit{H}}_1^0=(H_{1x}^0,0,H_{1z}^0)$, and ${\mathit{H}}_2^0=(H_{2x}^0,0,H_{2z}^0)$.

Substituting Equations (12) and (15) into the boundary condition (14), one obtains after the elementary algebra:

$$H_{2x}^0-H_{0x}^0-H_{1x}^0=\frac{4\pi }{c}{\sigma }^{2\mathrm{D},\mathrm{Tr}}(\mathit{q},\omega )E_{2y}^0.$$

(16)

From the Maxwell equation (5), written for the incident wave in free space where ${\mathit{B}}^0={\mathit{H}}^0$, one obtains

$${\mathit{H}}_0^0=\frac{c}{\omega }[{\mathit{k}}_0\times {\mathit{E}}_0^0].$$

(17)

With account of ${\mathit{k}}_0=(k_{0x},0,k_{0z})$ and ${\mathit{E}}_0^0=(0,E_{0y}^0,0)$, this reduces to

$$H_{0x}^0=-\frac{c}{\omega }{k}_{0z}{E}_{0y}^0,$$

(18)

where, considering Equation (11) and using that in this case $k_{0x}$ plays the role of q,

$$k_{0z}=\frac{\omega }{c}cos{\theta }_0=\sqrt{\frac{{\omega }^2}{c^2}-q^2}.$$

(19)

In a similar way, from Equation (5) written for the reflected and transmitted waves, one finds

$$H_{1x}^0=-\frac{c}{\omega }{k}_{0z}{E}_{1y}^0,H_{2x}^0=-\frac{c}{\omega }{k}_{0z}{E}_{2y}^0.$$

(20)

Substituting Equations (18) and (20) into Equation (16), we finally obtain

$$E_{0y}^0-E_{1y}^0-E_{2y}^0=\frac{4\pi \omega }{c^2{k}_{0z}}{\sigma }^{2\mathrm{D},\mathrm{Tr}}(\mathit{q},\omega )E_{2y}^0.$$

(21)

By solving this equation together with Equation (13), the transverse electric (TE) reflection coefficient is found in the form

$$r_{\mathrm{TE}}(\mathit{q},\omega )=\frac{E_{1y}^0}{E_{0y}^0}=-\frac{2\pi \omega {\sigma }^{2\mathrm{D},\mathrm{Tr}}(\mathit{q},\omega )}{c^2{k}_{0z}+2\pi \omega {\sigma }^{2\mathrm{D},\mathrm{Tr}}(\mathit{q},\omega )}.$$

(22)

Note that in Ref. [52] this reflection coefficient was obtained in the special case of normal incidence with ignored spatial dispersion.

By taking into account that for a 2D sheet, the spatially nonlocal dielectric permittivity is expressed via the conductivity as [13,31]

$${\epsilon }^{2\mathrm{D},\mathrm{Tr}}(\mathit{q},\omega )=1+\frac{2\pi i{\sigma }^{2\mathrm{D},\mathrm{Tr}}(\mathit{q},\omega )q}{\omega },$$

(23)

and using Equation (19), we rewrite the reflection coefficient (22) in the final form

$$r_{\mathrm{TE}}(\mathit{q},\omega )=-\frac{{\omega }^2[{\epsilon }^{2\mathrm{D},\mathrm{Tr}}(\mathit{q},\omega )-1]}{i{c}^{2}q\sqrt{\frac{{\omega }^2}{c^2}-q^2}+{\omega }^2[{\epsilon }^{2\mathrm{D},\mathrm{Tr}}(\mathit{q},\omega )-1]}.$$

(24)

This is the transverse electric Fresnel reflection coefficient on a 2D graphene sheet expressed via the spatially nonlocal transverse dielectric permittivity of graphene.

We now proceed to a derivation of the transverse magnetic reflection coefficient on a graphene sheet. In this case, the amplitudes of the magnetic field of the incident, ${\mathit{H}}_0^0$, transmitted, ${\mathit{H}}_2^0$, and reflected, ${\mathit{H}}_1^0$, waves are perpendicular to the plane of incidence and directed along the positive direction of the y-axis (see Figure 3). The amplitudes of the electric field, ${\mathit{E}}_0^0$, ${\mathit{E}}_2^0$, and ${\mathit{E}}_1^0$, lie in the plane of incidence.

According to Figure 3,

$${\mathit{H}}_{\mathrm{I}}^0={\mathit{H}}_0^0+{\mathit{H}}_1^0,{\mathit{H}}_{\mathrm{II}}^0={\mathit{H}}_2^0,$$

(25)

where ${\mathit{H}}_0^0=(0,H_{0y}^0,0)$, ${\mathit{H}}_1^0=(0,H_{1y}^0,0)$, and ${\mathit{H}}_2^0=(0,H_{2y}^0,0)$.

Taking into account that, in this case,

$$E_{\mathrm{II}}^0={\mathit{E}}_2^0=(E_{2x}^0,0,E_{2z}^0),$$

(26)

one obtains

$$\begin{array}{ccc}& & [{\mathit{E}}_2^0\times \mathit{n}]=(0,-E_{2x}^0,0),\hfill \\ & & \mathit{n}\times [{\mathit{E}}_2^0\times \mathit{n}]=(E_{2x}^0,0,0)=(E_2^{0}cos{\theta }_0,0,0),\hfill \end{array}$$

(27)

where $cos{\theta }_0$ is defined in Equation (19).

The boundary condition (9), (10) takes the form

$$({\mathit{H}}_{\mathrm{II}}^0-{\mathit{H}}_{\mathrm{I}}^0)\times \mathit{n}=-\frac{4\pi }{c}{\sigma }^{2\mathrm{D},\mathrm{L}}(\mathit{q},\omega )[\mathit{n}\times [E_{\mathrm{II}}^0\times \mathit{n}]]$$

(28)

The longitudinal conductivity ${\sigma }^{2\mathrm{D},\mathrm{L}}$ appears in this equation because the 2D wave vector q is now parallel to ${\mathit{E}}_{\mathrm{lat}}^0$.

Substituting Equations (25) and (27) into the boundary condition (28), one finds

$$H_{2y}^0-H_{0y}^0-H_{1y}^0=-\frac{4\pi {\sigma }^{2\mathrm{D},\mathrm{L}}(\mathit{q},\omega )}{c}{E}_2^{0}cos{\theta }_0.$$

(29)

For the transverse magnetic polarization, the boundary condition (6) reduces to

$$E_{0x}^0+E_{1x}^0-E_{2x}^0=0.$$

(30)

With account of Equation (11), which is valid for both polarizations of the electromagnetic field, Equation (30) is equivalent to

$$E_0^{0}cos{\theta }_0-E_1^{0}cos{\theta }_0-E_2^{0}cos{\theta }_0=0$$

(31)

and finally to

$$E_0^0-E_1^0-E_2^0=0.$$

(32)

Let us now use the Maxwell equation (5) for the incident wave ${\mathit{H}}_0^0={\mathit{B}}_0^0=(0,H_{0y}^0,0)$. Then it takes the form of Equation (17). By using ${\mathit{E}}_0^0=(E_{0x}^0,0,E_{0z}^0)$, one obtains from Equation (17) with the help of Equations (11) and (19)

$$H_{0y}^0=\frac{c}{\omega }(k_{0z}{E}_{0x}^0-k_{0x}{E}_{0z}^0)=\frac{c}{\omega }(k_{0z}cos{\theta }_0+k_{0x}sin{\theta }_0)E_0^0=E_0^0.$$

(33)

In a similar way, from the Maxwell Equation (17) applied to the reflected and transmitted waves, one obtains

$$H_{1y}^0=E_1^0,H_{2y}^0=E_2^0.$$

(34)

Substituting Equations (33) and (34) into Equations (29) and (32), one finds

$$\begin{array}{ccc}& & H_{2y}^0-H_{0y}^0-H_{1y}^0=-\frac{4\pi {\sigma }^{2\mathrm{D},\mathrm{L}}(\mathit{q},\omega )}{c}{H}_{2y}^{0}cos{\theta }_0,\hfill \\ & & H_{0y}^0-H_{1y}^0-H_{2y}^0=0.\hfill \end{array}$$

(35)

By solving these equations together, we derive the transverse magnetic (TM) reflection coefficient on a 2D graphene sheet:

$$r_{\mathrm{TM}}(\mathit{q},\omega )=\frac{H_{1y}^0}{H_{0y}^0}=\frac{2\pi {\sigma }^{2\mathrm{D},\mathrm{L}}(\mathit{q},\omega )cos{\theta }_0}{c+2\pi {\sigma }^{2\mathrm{D},\mathrm{L}}(\mathit{q},\omega )cos{\theta }_0}.$$

(36)

The result (36) was obtained in Refs. [50,51,52] in the spatially local case. In the presence of spatial dispersion, both the results (22) and (36) are presented in Refs. [13,31] with no detailed derivation.

By using an expression for the longitudinal dielectric permittivity of a 2D sheet through its conductivity [13,31]

$${\epsilon }^{2\mathrm{D},\mathrm{L}}(\mathit{q},\omega )=1+\frac{2\pi i{\sigma }^{2\mathrm{D},\mathrm{L}}(\mathit{q},\omega )q}{\omega },$$

(37)

and Equation (19), one obtains

$$r_{\mathrm{TM}}(\mathit{q},\omega )=\frac{[{\epsilon }^{2\mathrm{D},\mathrm{L}}(\mathit{q},\omega )-1]\sqrt{\frac{{\omega }^2}{c^2}-q^2}}{iq+[{\epsilon }^{2\mathrm{D},\mathrm{L}}(\mathit{q},\omega )-1]\sqrt{\frac{{\omega }^2}{c^2}-q^2}}.$$

(38)

This is the transverse magnetic Fresnel reflection coefficient on a 2D graphene sheet expressed via the longitudinal dielectric permittivity of graphene.

3. Spatially Nonlocal Dielectric Permittivities of Graphene and the Polarization Tensor

It is common knowledge that at low energies (smaller than approximately 3 eV [57]) graphene is well described by the Dirac model as a set of massless quasiparticles satisfying the Dirac equation, where the speed of light c is replaced with the Fermi velocity $v_F\approx c/300$ [1,2,3]. In the framework of the Dirac model, it is possible to derive explicit expressions for the polarization tensor of graphene, which describes the response of a graphene sheet to the electromagnetic field [40,41,42,43], and thus find both the transverse and longitudinal permittivities of graphene starting from the first principles of quantum electrodynamics. The dielectric permittivities obtained in this way depend on the wave vector, on the frequency, and also on temperature.

The polarization tensor of graphene in (2+1)-dimensional space-time is notated as ${\mathsf{\Pi }}_{\mu \nu }(\mathit{q},\omega )$, where $\mu ,\nu =0,1,2$ and the dependence on temperature is implied (here, we consider the pristine graphene sheet with no energy gap in the spectrum of quasiparticles and perfect hexagonal crystal lattice). The transverse dielectric permittivity of graphene is expressed as [58]

$${\epsilon }^{2\mathrm{D},\mathrm{Tr}}(\mathit{q},\omega )-1=-\frac{c^2}{2\hslash q{\omega }^2}\mathsf{\Pi }(\mathit{q},\omega ),$$

(39)

where the quantity $\mathsf{\Pi }$ is the following combination of the components of the polarization tensor:

$$\mathsf{\Pi }(\mathit{q},\omega )\equiv q^2{\mathsf{\Pi }}_{\mu }^{\mu }(\mathit{q},\omega )+\left(\frac{{\omega }^2}{c^2}-q^2\right){\mathsf{\Pi }}_{00}(\mathit{q},\omega ),$$

(40)

ℏ is the reduced Planck’s constant and the summation is made over the repeated indices.

The longitudinal dielectric permittivity of graphene is immediately expressed via the 00 component of the polarization tensor [58]

$${\epsilon }^{2\mathrm{D},\mathrm{L}}(\mathit{q},\omega )-1=\frac{c^2}{2\hslash q}{\mathsf{\Pi }}_{00}(\mathit{q},\omega ),$$

(41)

The polarization tensor of graphene along the real frequency axis was obtained in Ref. [42]. It was considered for the propagating waves, which satisfy the condition

$$q⩽\frac{\omega }{c},$$

(42)

and for the evanescent waves, which satisfy either the condition

$$\frac{\omega }{c}<q⩽\frac{\omega }{v_F}\approx 300\frac{\omega }{c}$$

(43)

(the so-called plasmonic region [59]) or the condition

$$q>\frac{\omega }{v_F}\approx 300\frac{\omega }{c}.$$

(44)

Using the expression from Ref. [42] for $\mathsf{\Pi }$ in the region of propagating waves (42) and in the plasmonic region (43), the transverse dielectric permittivity of graphene (39) in these regions can be written in the same form:

$$\begin{array}{ccc}& & {\epsilon }^{2\mathrm{D},\mathrm{Tr}}(\mathit{q},\omega )-1=i\pi \alpha q\frac{c}{2{\omega }^2}\sqrt{{\omega }^2-v_F^2{q}^2}\hfill \\ & & -\frac{8\alpha c^2}{v_F^{2}q}\left\{\underset{0}{\overset{u^{(-)}}{\int }}\frac{du}{e^{\beta u}+1}\left[1-\frac{1}{2{\omega }^2}\sqrt{{\omega }^2-v_F^2{q}^2}\sum _{\lambda =\pm 1}B(2cu+\lambda \omega )\right]\right.\hfill \\ & & \left.+\underset{u^{(-)}}{\overset{\infty }{\int }}\frac{du}{e^{\beta u}+1}\left[1-\frac{1}{2{\omega }^2}\sqrt{{\omega }^2-v_F^2{q}^2}\sum _{\lambda =\pm 1}\lambda B(2cu+\lambda \omega )\right]\right\}.\hfill \end{array}$$

(45)

Here and below, it is assumed that $\omega >0$ and the following notations are introduced:

$$u^{(-)}=\frac{1}{2c}(\omega -v_{F}q),\beta =\frac{\hslash c}{k_{B}T},B\left(x\right)=\frac{x^2}{\sqrt{x^2-v_F^2{q}^2}},$$

(46)

$\alpha$ denotes the fine structure constant, $k_B$ is the Boltzmann constant and T is the temperature of a graphene sheet.

In the region (44), using the corresponding expression for $\mathsf{\Pi }$ [42], one obtains another expression for the transverse dielectric permittivity of graphene

$$\begin{array}{ccc}& & {\epsilon }^{2\mathrm{D},\mathrm{Tr}}(\mathit{q},\omega )-1=-\pi \alpha q\frac{c}{2{\omega }^2}\sqrt{v_F^2{q}^2-{\omega }^2}\hfill \\ & & -\frac{4\alpha c}{v_F^{2}q}\sqrt{v_F^2{q}^2-{\omega }^2}\underset{0}{\overset{\infty }{\int }}\frac{dw}{e^{Dw}+1}\left[1-\frac{1}{2}\sum _{\lambda =\pm 1}\frac{{(\sqrt{v_F^2{q}^2-{\omega }^2}w+\lambda \omega )}^2}{{\omega }^2\sqrt{1-w^2-\frac{2\lambda \omega w}{\sqrt{v_F^2{q}^2-{\omega }^2}}}}\right],\hfill \end{array}$$

(47)

where $D=\hslash \sqrt{v_F^2{q}^2-{\omega }^2}/\left(2k_{B}T\right)$.

In a similar way, using the expression from Ref. [42] for ${\mathsf{\Pi }}_{00}$ in the region of propagating (42) and plasmonic (43) wave vectors and frequencies, one finds the explicit form of the longitudinal dielectric permittivity of graphene (41) in these regions:

$$\begin{array}{ccc}& & {\epsilon }^{2\mathrm{D},\mathrm{L}}(\mathit{q},\omega )-1=i\pi \alpha cq\frac{1}{2\sqrt{{\omega }^2-v_F^2{q}^2}}\hfill \\ & & +\frac{8\alpha c^2}{v_F^{2}q}\left\{\underset{0}{\overset{u^{(-)}}{\int }}\frac{du}{e^{\beta u}+1}\left[1-\frac{1}{2\sqrt{{\omega }^2-v_F^2{q}^2}}\sum _{\lambda =\pm 1}F(2cu+\lambda \omega )\right]\right.\hfill \\ & & \left.+\underset{u^{(-)}}{\overset{\infty }{\int }}\frac{du}{e^{\beta u}+1}\left[1-\frac{1}{2\sqrt{{\omega }^2-v_F^2{q}^2}}\sum _{\lambda =\pm 1}\lambda F(2cu+\lambda \omega )\right]\right\},\hfill \end{array}$$

(48)

where

$$F\left(x\right)=\sqrt{x^2-v_F^2{q}^2}.$$

(49)

Using the expression of ${\mathsf{\Pi }}_{00}$ [42] in the region (44), for the longitudinal permittivity of graphene (41) in this region, one obtains

$$\begin{array}{ccc}& & {\epsilon }^{2\mathrm{D},\mathrm{L}}(\mathit{q},\omega )-1=\pi \alpha cq\frac{1}{2\sqrt{v_F^2{q}^2-{\omega }^2}}\hfill \\ & & +\frac{4\alpha c}{v_F^{2}q}\sqrt{v_F^2{q}^2-{\omega }^2}\underset{0}{\overset{\infty }{\int }}\frac{dw}{e^{Dw}+1}\left[1-\frac{1}{2}\sum _{\lambda =\pm 1}\sqrt{1-w^2-\frac{2\lambda \omega w}{\sqrt{v_F^2{q}^2-{\omega }^2}}}\right].\hfill \end{array}$$

(50)

Thus, both the transverse and longitudinal dielectric permittivities of graphene are obtained in all ranges of the wave vectors and frequencies (42)–(44). We emphasize that the first lines of Equations (45), (47), (48) and (50) represent the corresponding dielectric permittivity at zero temperature. The terms in the next lines of these equations define the thermal correction to it found in the framework of the Dirac model. These terms make a profound effect on the reflectivity [42] and conductivity [60] properties of graphene, and also on the Casimir interaction between graphene sheets [44,46,47,48]. By construction from the polarization tensor, the obtained permittivities satisfy the Kramers–Kronig relations. The specific form of these relations was investigated in the spatially local limit $q\to 0$ [61] and at zero temperature [62].

For the calculation of the Casimir force in graphene systems, it is helpful to use the reflection coefficients (24) and (38), as well as the dielectric permittivities of graphene, written at the pure imaginary Matsubara frequencies $\omega =i{\xi }_l=2\pi i{k}_{B}Tl/\hslash$, where $l=0,1,2,\dots$.

Substituting $\omega =i{\xi }_l$ into Equations (24) and (38), one obtains, respectively,

$$\begin{array}{ccc}& & r_{\mathrm{TE}}(\mathit{q},i{\xi }_l)=-\frac{{\xi }_l^2[{\epsilon }^{2\mathrm{D},\mathrm{Tr}}(\mathit{q},i{\xi }_l)-1]}{c^{2}q\sqrt{q^2+\frac{{\xi }_l^2}{c^2}}+{\xi }_l^2[{\epsilon }^{2\mathrm{D},\mathrm{Tr}}(\mathit{q},i{\xi }_l)-1]},\hfill \\ & & r_{\mathrm{TM}}(\mathit{q},i{\xi }_l)=\frac{[{\epsilon }^{2\mathrm{D},\mathrm{Tr}}(\mathit{q},i{\xi }_l)-1]\sqrt{q^2+\frac{{\xi }_l^2}{c^2}}}{q+[{\epsilon }^{2\mathrm{D},\mathrm{Tr}}(\mathit{q},i{\xi }_l)-1]\sqrt{q^2+\frac{{\xi }_l^2}{c^2}}}.\hfill \end{array}$$

(51)

These are the Fresnel reflection coefficients in two dimensions calculated at the pure imaginary Matsubara frequencies. The same expressions are obtained if one substitutes Equations (39) and (41) into the reflection coefficients derived in Refs. [40,41] directly in terms of the polarization tensor.

The spatially nonlocal dielectric permittivities of graphene along the imaginary frequency axis are immediately obtainable from Equations (47) and (50) valid in the interval (44) by putting $\omega =i{\xi }_l$. The results are

$$\begin{array}{ccc}& & {\epsilon }^{2\mathrm{D},\mathrm{Tr}}(\mathit{q},i{\xi }_l)-1=\pi \alpha q\frac{c}{2{\xi }_l^2}\sqrt{v_F^2{q}^2+{\xi }_l^2}\hfill \\ & & -\frac{4\alpha c}{v_F^{2}q}\sqrt{v_F^2{q}^2+{\xi }_l^2}\underset{0}{\overset{\infty }{\int }}\frac{dw}{e^{D_{l}w}+1}\left[1+\frac{1}{2}\sum _{\lambda =\pm 1}\frac{{(\sqrt{v_F^2{q}^2+{\xi }_l^2}w+i\lambda {\xi }_l)}^2}{{\xi }_l^2\sqrt{1-w^2-\frac{2i\lambda {\xi }_{l}w}{\sqrt{v_F^2{q}^2+{\xi }_l^2}}}}\right],\hfill \\ & & {\epsilon }^{2\mathrm{D},\mathrm{L}}(\mathit{q},i{\xi }_l)-1=\pi \alpha cq\frac{1}{2\sqrt{v_F^2{q}^2+{\xi }_l^2}}\hfill \\ & & +\frac{4\alpha c}{v_F^{2}q}\sqrt{v_F^2{q}^2+{\xi }_l^2}\underset{0}{\overset{\infty }{\int }}\frac{dw}{e^{D_{l}w}+1}\left[1-\frac{1}{2}\sum _{\lambda =\pm 1}\sqrt{1-w^2-\frac{2i\lambda {\xi }_{l}w}{\sqrt{v_F^2{q}^2+{\xi }_l^2}}}\right],\hfill \end{array}$$

(52)

where now $D_l=\hslash \sqrt{v_F^2{q}^2+{\xi }_l^2}/\left(2k_{B}T\right)$.

These expressions are indeed real as it should be. The same dielectric permittivities are obtained at once from Equations (39) and (41) written at $\omega =i{\xi }_l$ when substituting expressions for $\mathsf{\Pi }(\mathit{q},i{\xi }_l)$ and ${\mathsf{\Pi }}_{00}(\mathit{q},i{\xi }_l)$ derived directly along the imaginary frequency axis [63] rather than analytically continued from the real frequency axis as it was made above.

4. Contribution of Different Polarizations and the Role of Evanescent Waves in the Casimir Pressure between Two Graphene Sheets

The Casimir pressure between two parallel graphene sheets at temperature T separated by distance a is given by the Lifshitz formula, which can be presented in terms of either pure imaginary Matsubara or real frequencies [8,12]. In both cases, the total pressure is the sum of contributions from the electromagnetic waves of TM and TE polarizations.

We begin from the representation in terms of the Matsubara frequencies

$$P(a,T)=P_{\mathrm{TM}}(a,T)+P_{\mathrm{TE}}(a,T),$$

(53)

where

$$P_{\mathrm{TM},\mathrm{TE}}(a,T)=-\frac{k_{B}T}{\pi }\sum _{l=0}^{\infty }{\phantom{\sum }}^{\prime }\underset{0}{\overset{\infty }{\int }}dqq\sqrt{q^2+\frac{{\xi }_l^2}{c^2}}{\left[r_{\mathrm{TM},\mathrm{TE}}^{-2}(\mathit{q},i{\xi }_l)e^{2a\sqrt{q^2+\frac{{\xi }_l^2}{c^2}}}-1\right]}^{-1}.$$

(54)

Here, the prime on the summation sign adds the factor 1/2 to the term with $l=0$, and the reflection coefficients on a graphene sheet for both polarizations are defined in Equation (51) with the dielectric permittivities of graphene presented in Equation (52).

We performed computations of both $P_{\mathrm{TM}}$ and $P_{\mathrm{TE}}$ in the application region of the Dirac model, i.e., under a condition that the characteristic energy of the Casimir force $\hslash {\omega }_c=\hslash c/\left(2a\right)$ should be less than 3 eV [57]. This condition is well satisfied at $a⩾200$ nm, where $\hslash {\omega }_c⩽0.5$ eV.

The computational results for the magnitudes of $P_{\mathrm{TM}}$ and $P_{\mathrm{TE}}$ at $T=300$ K are presented in Figure 4 in the logarithmic scale by the upper and lower lines, respectively, as the function of separation between the graphene sheets. Both $P_{\mathrm{TM}}$ and $P_{\mathrm{TE}}$ are negative, i.e., they contribute to the Casimir attraction.

As seen in Figure 4, the major contribution to the Casimir pressure at $a⩾200$ nm is given by the transverse magnetic polarization, whereas the transverse electric one makes only a negligible small contribution. Thus, at $a=200$ nm, we have $P_{\mathrm{TM}}/P_{\mathrm{TE}}=1530$, i.e., $P_{\mathrm{TM}}/P=0.99935$. The role of the TM polarization only increases with increasing separation. As two more examples, at $a=2$ and $4\mathsf{\mu }$m one finds that $P_{\mathrm{TM}}/P_{\mathrm{TE}}=1.92\times {10}^4$ and $1.49\times {10}^5$, respectively. This results in the following respective fractions of $P_{\mathrm{TM}}$ in the total Casimir pressure: $P_{\mathrm{TM}}/P=0.99995$ and 0.999993.

In Section 3, devoted to the nonlocal dielectric permittivities of graphene, they were considered in the region of propagating (42) and evanescent (43) and (44) waves. In so doing, within the region of propagating (42) and in the plasmonic subregion (43) of evanescent waves, these permittivities have a common analytic form. Nevertheless, keeping in mind an especially important role of the propagating waves (42), which are on the mass shell in free space, it is appropriate to consider their contribution to the Casimir pressure separately. Then, the contribution of the evanescent waves is computed as a sum of two terms using two different forms of the dielectric functions depending on whether the condition (43) or (44) is satisfied. Such a separation into the propagating and evanescent waves is also dictated by the form of the Lifshitz formula written in terms of real frequencies (see below).

The representation mathematically equivalent to Equations (53) and (54) of the Lifshitz formula in terms of real frequencies can be written in the form

$$P(a,T)=P_{\mathrm{TM}}^{\mathrm{prop}}(a,T)+P_{\mathrm{TE}}^{\mathrm{prop}}(a,T)+P_{\mathrm{TM}}^{\mathrm{evan}}(a,T)+P_{\mathrm{TE}}^{\mathrm{evan}}(a,T).$$

(55)

Here, the contributions of the propagating waves with different polarizations are given by [8,12]

$$\begin{array}{ccc}& & P_{\mathrm{TM},\mathrm{TE}}^{\mathrm{prop}}(a,T)=-\frac{\hslash }{2{\pi }^2}\underset{0}{\overset{\infty }{\int }}d\omega coth\frac{\hslash \omega }{2k_{B}T}\underset{0}{\overset{\omega /c}{\int }}qdq\hfill \\ & & \times \mathrm{Im}\left\{\sqrt{q^2-\frac{{\omega }^2}{c^2}}{\left[r_{\mathrm{TM},\mathrm{TE}}^{-2}(\mathit{q},\omega )e^{2a\sqrt{q^2-\frac{{\omega }^2}{c^2}}}-1\right]}^{-1}\right\},\hfill \end{array}$$

(56)

where the reflection coefficients are defined in Equations (24) and (38) and the dielectric permittivities in the region (42) are given by Equations (45) and (48).

The contributions of evanescent waves to Equation (55) with different polarizations take the form [8,12]

$$\begin{array}{ccc}& & P_{\mathrm{TM},\mathrm{TE}}^{\mathrm{evan}}(a,T)=-\frac{\hslash }{2{\pi }^2}\underset{0}{\overset{\infty }{\int }}d\omega coth\frac{\hslash \omega }{2k_{B}T}\underset{\omega /c}{\overset{\infty }{\int }}qdq\sqrt{q^2-\frac{{\omega }^2}{c^2}}\hfill \\ & & \times \mathrm{Im}{\left[r_{\mathrm{TM},\mathrm{TE}}^{-2}(\mathit{q},\omega )e^{2a\sqrt{q^2-\frac{{\omega }^2}{c^2}}}-1\right]}^{-1},\hfill \end{array}$$

(57)

where the reflection coefficients are again defined in Equations (24) and (38). As to the dielectric permittivities entering these reflection coefficients, in the region (43), they are given by Equations (45) and (48), but in the region (44), by Equations (47) and (50).

Equations (56) and (57) are not as convenient for computations as Equation (54). This is most pronounced in $P_{\mathrm{TM},\mathrm{TE}}^{\mathrm{prop}}$ defined in Equation (56), which contains the quickly oscillating functions due to the pure imaginary power in the exponential factor. As to Equation (57), the power of the exponent remains real.

Taking into account that in the application region of the Dirac model nearly the total Casimir pressure is determined by the TM polarized waves, we compute the quantity $P_{\mathrm{TM}}^{\mathrm{evan}}$ by Equations (57) and (38) using the dielectric permittivities defined in Equations (48) and (50). As to the quantity $P_{\mathrm{TM}}^{\mathrm{prop}}$, it is more convenient to not compute it directly by Equation (56), but determine it as a difference

$$P_{\mathrm{TM}}^{\mathrm{prop}}(a,T)=P_{\mathrm{TM}}(a,T)-P_{\mathrm{TM}}^{\mathrm{evan}}(a,T),$$

(58)

where $P_{\mathrm{TM}}$ is already computed by the Lifshitz formula (54) written in terms of the Matsubara frequencies.

The computational results for $P_{\mathrm{TM}}$, $P_{\mathrm{TM}}^{\mathrm{evan}}$, and $P_{\mathrm{TM}}^{\mathrm{prop}}$ at $T=300$ K normalized to the Casimir pressure between two ideal metal plates in the classical limit [12]

$$P_{\mathrm{IM}}(a,T)=-\frac{k_{B}T}{4\pi a^3}\zeta \left(3\right),$$

(59)

where $\zeta \left(z\right)$ is the Riemann zeta function, are presented in Figure 5 as the function of separation by the solid, long-dashed, and short-dashed lines, respectively.

According to Figure 5, at separations of 200–400 nm both the evanescent and propagating transverse magnetic waves contribute significantly to the Casimir pressure. At larger separations, the dominant contribution is given by the evanescent waves. In doing so, the contribution of evanescent waves is attractive at all separation distances. Calculations show, however, that this attraction is combined from the attractive part caused by the plasmonic region (43) and the repulsive part caused by the region (44). The contribution of the TM propagating waves to the Casimir pressure between two graphene sheets changes its character from attraction to repulsion and vice versa with increasing separation.

5. Discussion: Whether Graphene Helps to Solve the Problem Arising for Real Metals

The main distinctive feature of the Casimir pressure in the configuration of two graphene sheets considered above is that in the framework of the Dirac model, the spatially nonlocal dielectric permittivities of graphene are found precisely starting from the first principles of thermal quantum field theory. As to the dielectric permittivities of metals used in computations by means of the Lifshitz formula, they contain phenomenological parameters, such as the relaxation parameter of the Drude model, and have not been tested experimentally within all frequency regions essential for the Casimir effect (i.e., in the region of transverse electric evanescent waves).

The formalism of the Lifshitz theory for two graphene sheets presented in Section 2 and Section 3 is in perfect analogy with that commonly used for two metallic plates. The Lifshitz formula for the Casimir pressure remains unchanged, and only the 3D Fresnel reflection coefficients are replaced with their 2D analogues as it should be done when considering the Casimir interaction of plane structures. Taking into account the fundamental character of the Lifshitz theory, we obtain the conclusion that only some drawback in the used response functions of metals to the electromagnetic field could cause a disagreement of the theoretical predictions with measurements of the Casimir interaction between Au surfaces.

As shown in Section 4, for two graphene sheets, the total Casimir pressure is determined by the contribution of only the transverse magnetic waves. This is because in the application region of the Dirac model at $a⩾200$ nm the Casimir force between graphene sheets is already in the classical limit where the contributions of the TE polarized propagating and evanescent waves cancel each other. The same occurs for the Casimir force between metallic plates described by the Drude model at separations exceeding the thermal length [64], i.e., larger than $7.6\mathsf{\mu }$m at room temperature. At so large separations, however, there are no reliable measurement data available. As to the experimental separations between metallic plates, both the TM and TE polarizations contribute to the Casimir pressure irrespective of whether the experimentally consistent plasma model or the Drude model excluded by the measurement data is used [18].

By and large, the case of graphene suggests to us that when calculating the Casimir force using the Lifshitz theory, it is important to adequately describe the response of boundary materials to both the propagating and evanescent waves with the transverse magnetic and transverse electric polarizations and take proper account of the effects of spatial dispersion.

6. Conclusions

In the foregoing, we considered the Casimir pressure between two graphene sheets using the Lifshitz theory in the form that is most frequently used for a description of the Casimir effect between conventional 3D materials. For this purpose, we presented the detailed derivation of the 2D Fresnel reflection coefficients on a graphene sheet with due account of the spatial dispersion. As a result, the reflection coefficients for two independent polarizations of the electromagnetic field were expressed via the transverse and longitudinal dielectric permittivities of graphene, which depend on the 2D wave vector, frequency, and temperature. These reflection coefficients are equivalent to those expressed directly via the polarization tensor of graphene.

Next, we presented the explicit expressions for the transverse and longitudinal dielectric permittivities of graphene along the real frequency axis in the regions of both the propagating and evanescent waves and also at the pure imaginary Matsubara frequencies. This was made using the polarization tensor of graphene, which was found earlier in the framework of the Dirac model.

Using the Lifshitz formula written in terms of the Matsubara frequencies, we demonstrated that the total Casimir pressure between two graphene sheets at separations exceeding 200 nm is fully determined by the TM polarized electromagnetic field. By applying the Lifshitz formula along the real frequency axis, the contributions of the TM polarized propagating and evanescent waves to the total pressure were found.

Finally, the above results obtained for graphene sheets were confronted with the corresponding results valid for two metallic plates. This confrontation points the way for bringing the Lifshitz theory in agreement with the measurement data by using the more accurate dielectric functions of metallic test bodies. In the future, it is planned to consider different contributions to the Casimir force between two real graphene sheets possessing the nonzero chemical potential, which prevents from reaching the classical limit at the experimental separations.