Quantum friction between oscillating crystal slabs: Graphene monolayers on dielectric substrates

We present a theoretical description of energy transfer processes between two noncontact quasi- twodimensional crystals separated by distance a, oscillating with frequency omega0 and amplitude rho0 , and compare it with the case of two quasi-twodimensional crystals in uniform parallel motion. We apply the theory to calculate van der Waals energy and dissipated energy in two oscillating slabs where each slab consists of a graphene monolayer deposited on SiO2 substrate. The graphene dielectric response is determined from first principles, and SiO2 surface response is described using empirical local dielectric function. We studied the modification of vdW attraction as function of the driving frequency and graphene doping. We propose the idea of controlling the sticking and unsticking of slabs by tuning the graphene dopings EF i and driving frequency omega0 . We found simple rho02 dependence of vdW and dissipated energy. As the Dirac plasmons are the dominant channels through which the energy between slabs can be transferred, the dissipated power in equally doped EF1 = EF2 = 0 graphenes shows strong omega0 = 2omegap peak. This peak is substantially reduceed when graphenes are deposited on SiO2 substrate. If only one graphene is pristine (EFi = 0) the 2omegap peak disappears. For larger separations a the phononic losses also become important and the doping causes shifts, appearance and disappearance of many peaks originating from resonant coupling between hybridized electronic-phononic excitations in graphene-substrate slabs.

We present a theoretical description of energy transfer processes between two noncontact quasitwodimensional crystals separated by distance a, oscillating with frequency ω0 and amplitude ρ0, and compare it with the case of two quasi-twodimensional crystals in uniform parallel motion. We apply the theory to calculate van der Waals energy and dissipated energy in two oscillating slabs where each slab consists of a graphene monolayer deposited on SiO2 substrate. The graphene dielectric response is determined from first principles, and SiO2 surface response is described using empirical local dielectric function. We studied the modification of vdW attraction as function of the driving frequency and graphene doping. We propose the idea of controlling the 'sticking' and 'unsticking' of slabs by tuning the graphene dopings EF i and driving frequency ω0. We found simple ρ 2 0 dependence of vdW and dissipated energy. As the Dirac plasmons are the dominant channels through which the energy between slabs can be transferred, the dissipated power in equally doped EF 1 = EF 2 = 0 graphenes shows strong ω0 = 2ωp peak. This peak is substantially reduceed when graphenes are deposited on SiO2 substrate. If only one graphene is pristine (EF i = 0) the 2ωp peak disappears. For larger separations a the phononic losses also become important and the doping causes shifts, appearance and disappearance of many peaks originating from resonant coupling between hybridized electronic/phononic excitations in graphene/substrate slabs.

I. INTRODUCTION
Detailed understanding of non-contact friction and energy transfer processes in nanostructures is of great importance, both from the conceptual and practical viewpoints. Existing theoretical studies, starting with the seminal paper by Pendry [1], mostly consist of calculations of friction coefficients, i.e. friction force between two parallel dielectric plates (e.g. supported graphenes) in uniform relative motion which is experimentally not easily measured (e.g. current drag in one graphene caused by current flow in another one) [2][3][4][5][6][7][8].
While the experiments with two slabs in parallel relative motion with constant velocity are difficult to perform, we suggest here that for the same systems experiments with slabs in relative oscillatory motion with fixed or variable frequency might be easier to perform, and could lead to new and interesting observations. Recently a similar approach has been realized experimentally [9][10][11][12]. In these experiments the system (usually an AFM tip above the surface) oscillates at some characteristic frequency. These oscillations are then, because of various dissipation mechanisms (which includes quantum friction), damped. Our model is based on a slightly different concept; one of the slabs, e.g. the AFM tip, is driven with variable frequency. This means that the friction can be deduced from the energy dissipated in one oscillating cycle. In this paper we provide a general the-oretical description of such processes, expecting that this method might become a useful tool to study dynamical properties of low-dimensional systems [13].
The main objective of this paper is therefore a theoretical description of these phenomena in systems consisting of two non-touching polarizable media, specifically conservative (van der Waals or Casimir) and dissipative forces (quantum friction) between two quasitwodimensional (q2D crystals) in relative parallel and oscillatory motion. While the case of slabs in parallel uniform motion has been extensively studied [1,5,[14][15][16][17][18], here we develop an analogous theory describing interaction of atomically thick slabs (q2D crystals) in oscillatory motion.
In Sec.II the expressions for van der Waals and dissipative energies and forces are derived for such a q2D system in a very general case, for variable slab temperatures and dynamical properties characterized by their surface response functions D 1 and D 2 , and for variable oscillating frequencies and amplitudes. We assume 2D translational invariance and neglect retardation for the slab distances in consideration. For the sake of clarity and comparison, in AppendixA we derive analogous results for the case of parallel uniform motion, recovering but also generalizing some earlier results [19,20].
In Sec.III we derive general expressions for surface response functions D i for multilayer slabs, later to be specified for monolayers of a substance like graphene or silicene adsorbed on dielectric substrates. Surface response functions D 1 and D 2 will be the key ingredients in the expressions describing dissipative and reactive processes in Sec.II and Sec.III. In Sec.III we also show how to calculate surface response functions D i for a specific case of q2D crystals on a dielectric substrate The expression for the surface excitation propagator of a system of two coupled slabs is also derived.
In Sec.IV we present the models used to describe the q2D crystal and substrate dynamical response. We study the specific case of a graphene monolayer on a dielectric substrate, which is chosen to be ionic crystal SiO 2 . The substrate is considered as a homogenous semiinfinite ionic crystal SiO 2 with the appropriate dielectric function in the longwavelength limit. Graphene monolayer dynamical response is determined from first principles. Also some computational details are specified.
In Sec.V general expressions of previous sections are applied to the system of two slabs, where each slab represents a graphene(E F i )/SiO 2 system, and where graphene doping is characterized by Fermi energy E F relative to the Dirac point.
In Sec.V A we demonstrate how the spectra of electronic excitations in one slab and in two coupled slabs depend on graphene doping E F .
The form of these coupled ecitations is responsible for the behaviour of the atractive forces and dissipation. We first discuss in Sec.V B the modification of van der Waals force for oscillating in comparison with the static slabs. Van der Waals energies depend on two factors. They increase with the increased graphene doping, but are reduced for the asymmetric doping when excitations in two slabs are off-resonance. Dynamical vdW energy shows unusual behavior: it starts as plateau, and then decreases. This is, because the fast Dirac plasmon in one slab for low driving frequencies ω 0 < ω p , still perfectly follows Doppler shifted charge density fluctuations in another slab. For larger driving frequencies this is not the case and vdW energies decrease. Finally, for small or zero doping the π → π * and π → σ excitations cause linear weakening of the dynamical vdW energy.
In Sec.V C we calculate and discuss how dissipated power depends on various parameters: driving amplitude ρ 0 and frequency ω 0 , on the separations between slabs a and on the substrate. We find simple ρ 2 0 dependence, while the ω 0 dependence is determined by the intensity of resonant coupling between hybridized Dirac plasmons and substrate TO phonons. We found that in realistic grahenes (in comparison with Drude model when excitation of undamped Dirac plasmons provides unrealistically strong 2ω p peak in the dissipated power) the dissipation power peak is strongly reduced and red shifted. We also explain why the substrate substantially reduces dissipated power peak. For larger separations a additional peaks appear in dissipated power originating from the excitations of hybridized substrate phonons.
In Sec.V D we explore how the dissipated power depends on graphene dopings. We show that if one graphene is pristine (E F = 0) it causes the disappearance a ω 0 In Appendix A 1 we have derived van der Waals energy and force between two slabs in uniform relative motion in some detail because it will help us to treat a similar problem of two oscillating slabs.
We shall later assume that the slabs consist of graphene monolayers with variable doping, deposited on dielectric slabs of thickness ∆ described by local dielectric functions (ω), as shown in Fig.1. The left slab mechanically oscillates with frequency ω 0 and amplitude ρ 0 relative to the right slab. Again we calculate the diagram in Fig.8 as in the A 1, but now the slab parallel coordinates change in time as so that instead of (A.3) we have If we use where J m are Bessel functions, after Fourier transformation in ω space, using expressions (A.5-A.7), (A.9) and integration over z coordinates we obtain Here we have also used the fact that ImD 2 (Q, ω) is an antisymmetric function of ω and does not contribute to integration. We see that the energy oscillates in time with frequencies (m − m )ω 0 . If we assume to measure energies on a time scale ∆t > T , where T = 2π ω0 is the maximal duration of one cycle, then we can average over T and find the result independent of time: where the expression in curly brackets is fully analogous to the one in (A.11), but now ω → ω m = ω+mω 0 . Inclusion of higher order processes follows the same procedure as for the parallel motion in A 1. After integration over the coupling constant, we obtain the result analogous to (A.16) where A is given by (A.17) and (A.18), with ω m = ω + mω 0 . Again, the limiting cases can be obtained from Sec.A 1. For ω 0 = 0 (ω = ω) and ρ 0 = 0 we find the well known result for van der Waals interaction when the slabs are at rest [21,22]: For finite frequency ω 0 and D 1 = D 2 = D we find: We notice that the frequency integrals are the same as in (A. 16-A.20). Also, the attractive van der Waals force between two oscillating slabs is given by where the function B is given by

B. Dissipated power
We can perform the calculation of the dissipated power for two slabs oscillating parallel to each other with amplitude ρ 0 and frquency ω 0 in analogy with the previous treatment of two slabs in uniform relative motion in Sec.A 2. Again, we have to transform the parallel coordinates in the left slabs as in (1). Then (A.29), after integration over t 1 becomes We see that the energy transfer rate is time dependent and oscillates with frequency (m − m)ω 0 . Again, from (2) we see that for time intervals large with respect to the oscillation period T the terms m = m do not contribute and the energy transfer rate is If we now use (A.5), the definitions (A.6) and (A.7) of the surface correlation function and the surface excitation propagator, respectively, and the connection (A.9) between the surface correlation function and the imaginary part of surface excitation propagator, equation (6) can be written as Evaluating (7) we have used the fact that the real part of the function under summation and integration is odd and the imaginary part is an even function of n and ω. P 12 is the energy transferred from the left to the right slab. Now we have to repeat the discussion in Sec.A 2 and substract the part of this energy which will be reversibly returned to the left slab. The same arguments, leading to (A.37), will give this energy to be Expression (8) represents the energy transferred from the left to right but which will be reversibly returned, as shown in Fig.10b. Therefore the energy which is irreversibly transferred from the left to the right, i.e. the dissipated power, is Analogous calculation would give the energy dissipated in the process where the charge fluctuation in the right slab induces fluctuations in the left slab. We have to exchange 1 and 2 in (9) and replace m → −m. Repeating the steps in (A.40) the final result becomes: This expression is analogous to (A.40). For T = 0 2n(ω) + 1 → sgnω and (10) can be written as Adding higher order terms (A.12,A.13) we obtain the energy dissipated per unit time: where C is given by (A.42). Limiting casses are also obtained from (12). For ω 0 = 0 and/or for ρ 0 = 0 obviously P = 0.
The main quantities which appear in the formula for van der Waals interaction E c or dissipated power P are the surface excitation propagators D 1 (Q, ω) and D 2 (Q, ω) of the left (first) and right (second) slab, respectively. The derivation of D 1 and D 2 is analogous for both slabs, so here we shall derive just one surface excitation propagator D. The structure of the monolayer-substrate composite (e.g. graphene on SiO 2 ) is shown in Fig.2. The slab consists of the graphene monolayer adsorbed at some small distance h (e.g. h = 0.4nm) above the substrate of macroscopic thickness ∆. The dielectric, e.g. the SiO 2 slab is placed in the region −∆ − h ≤ z ≤ −h and the graphene layer occupies z = 0 plane. The same model system is used in Refs. [23,24] where the authors explore plasmon-phonon hybridization, stopping power and wake effect produced by the proton moving parallel to the composite. The unit cell for such huge nanostructure would consist of hundreds of atoms, so it is impossible to perform full ab initio ground state and structure optimization calculation. Moreover, an ab initio calculation of the response function would be even more demanding so we need an approximation for the response function calculation. The easiest (and probably the best) approximation is to treat the SiO 2 slab as a homogeneous dielectric described by some local dielectric function S (ω) and to consider graphene as a purely 2D system described by the response function R(Q, ω), as sketched in Fig.2. In order to derive the surface excitation propagator D(Q, ω) we start from its definition: which connects the surface excitation propagator with the screened Coulomb interaction W (Q, ω, z = 0, z = 0) at z = z = 0 surface. Here R(Q, ω, z, z ) represents the nonlocal dielectric function of graphene/dielectric composite which we assume occupies the region z, z ≤ 0. It is well known [25][26][27][28] that physical properties of a graphene monolayer in the low (Q, ω) region can be described to a very good approximation assuming the monolayer to be strictly twodimensional, so that the nonlocal independent electron response function can be written as where we assume that the graphene lies in the z = 0 plane and the response function R 0 (Q, ω) can be derived from first principles, as decribed in Sec.IV. Dynamically screened response function R(Q, ω) in RPA is given as a series of terms .
If we assume for the moment that there is no dielectric in the system (e.g. S (ω) = 1) then the screened Coulomb interaction is simply given by Using the definition (13) the surface excitation propagator becomes When the dielectric slab is introduced, the external charges and charge density fluctuations in the graphene layer do not interact via the bare Coulomb interaction v Q but via the Columob interaction modified by the presence of the dielectric slab [29] v where the substrate surface excitation propagator is and represents the surface excitation propagator of a semiinfinite (∆ → ∞, h = 0) dielectric. This causes that the screened Colulomb interaction (16) becomes the function ofṽ Q (ω) where, because charge density fluctuations inside graphene also interact viaṽ Q (ω), the screened response function is modified as Finally, after inserting (21) into (13) we obtain the surface excitation propagator in the presence of the dielectric which can be rewritten in a more transparent form as The spectrum of coupled excitations in a single slab can be calculated from For the coupled slabs described by their surface excitations propagators D 1 and D 2 , separated by the distance a, in a similar way we can derive the propagatorD for the coupled system and the excitation spectrum of this system is

IV. DESCRIPTION OF SUBSTRATE AND GRAPHENE DYNAMICAL RESPONSE
The results in Sec.III are quite general and can be applied to a monolayer of any material on any dielectric substrate. Now we shall specify the dielectric substrate to be the homogenous layer of ionic crystal SiO 2 .
Dielectric properties (or dynamical response) of bulk ionic crystals in the long-wavelength limit can be described in terms of their optical phonons at the Γ point. More complex polar crystals such as SiO 2 possess a multitude of different optical phonons of different symmetries and polarizations. However, here we suppose that SiO 2 posses two well-defined, non-dispersing transverse optical (TO) phonon modes at the frequencies ω T O1 and ω T O2 with the corresponding damping rates γ T O1 and γ T O2 , giving rise to a dielectric function of the form [23,24] where 0 , i , and ∞ represent the dielectric constant for SiO 2 at the zero, intermediate, and very large frequencies. This dielectric function will be inserted in the expression (19) for the substrate surface excitation propagator D S (Q, ω). The graphene response function R(Q, ω) is given by (22) in terms of the noninteracting response function where the 3D Fourier transform of independent electron response function is given by [30] R 0 where Q is the momentum transfer vector parallel to the x − y plane, G = (G , G z ) are 3D reciprocal lattice vectors and r = (ρ, z) is a 3D position vector. Integration in (31) is performed over the normalization volume Ω = S × L, where S is the normalization surface and L is the superlattice constant in z direction (separation between graphene layers is superlattice arrangement). Plane wave expansion of the wave function has the form where the coefficients C nK are obtained by solving the Local Density Approximation-Kohn Sham (LDA-KS) equations selfconsistently as will be discussed below. However, this straightforward calculation of graphene response functions R(Q, ω) is not sufficient if we want to investigate the hybridization between the Dirac plasmon and Fuchs-Kliewer (FK) phonons at dielectric surfaces. Namely, due to the very low energy of FK phonons (∼ 50meV) the crossing of their dispersion relations with Dirac plasmon occurs for very small wave vectors (Q < 0.001a.u.). On the other hand even for very dense K-point mesh sampling, as for example 601×601×1 used in this calculation, the minimum transfer wave vector Q which can be reached (e.g. Q = 0.0026a.u. −1 in this calculation) is still bigger than FK phonon-Dirac plasmon crossing wave vector. Therefore we have to find the way how to calculate R(Q, ω) for a denser Q-point mesh in the optical Q ≈ 0 limit. One possible way is that instead of calculating response function R 0 (Q, ω) we calculate the optical (Q = 0) conductivity σ(ω). The optical conductivity in graphene can be written as [27] σ where is intraband or Drude conductivity and where represents the effective number of charge carriers. The interband conductivity is where the current vertices are given by and In the optical Q ≈ 0 limit the independent electron response function can be written in terms of optical conductivities (32) as [31] Finally, the RPA or screened response function R(Q, ω) can be obtained from (38) using (15).
In the calculation of Sec.V we shall assume the graphene response to be isotropic in the small (Q, ω) limit. This means that the graphene response functions and the corresponding surface excitation functions are functions of Q and not of Q.

A. Computational details
The first part of the calculation consists of determining the KS ground state of the single layer graphene and the corresponding wave functions φ nK (ρ, z) and energies E n (K). For graphene unit cell constant we use the experimental value of a = 4.651 a.u. [32], and superlattice unit cell constant (separation of graphene layers) is L = 5a. For calculating KS wave functions and energies we use a plane-wave self-consistent field DFT code (PWSCF) within the QUANTUM ESPRESSO (QE) package [33]. The core-electron interaction was approximated by the norm-conserving pseudopotentials [34], and the exchange correlation (XC) potential by the Perdew-Zunger local density approximation (LDA) [35]. To calculate the ground state electronic density we use 21 × 21 × 1 Monkhorst-Pack K-point mesh [36] of the first Brillouin zone (BZ) and for the plane-wave cut-off energy we choose 50 Ry. The second part of calculation consists of determining the independent electron response function (30) and conductivity (32)(33)(34)(35). In order to achieve better resolution in the long wavelength (Q ≈ 0) and low energy (ω ≈ 0) limit the response function (30,31) and conductivity (32)(33)(34)(35)(36)(37) are evaluated from the wave functions φ nK (r) and energies E n (K) calculated for the 601 × 601 × 1 Monkhorst-Pack K-point mesh which coresponds to 361801 K-points in the first Brillouin zone (1BZ). Band summations (n, m) in (30), (34) and (35) are performed over 30 bands. In the calculation we use two kinds of damping parameters: η intra = 10meV for transitions within the same bands (n ↔ n), and η inter = 50meV for transitions between different bands (n ↔ m). For bulk SiO 2 dielectric function given by (28) we use the following parameters: 0 = 3.9, i = 3.05, ∞ = 2.5, ω T O1 = 55.6 meV, ω T O2 = 138.1 meV, γ T O1 = 5.368 meV and γ T O2 = 8.947 meV taken from Ref. [37]. For the gap between graphene and the SiO 2 surface, we take h = 4Å[7.55 a.u.] [38].

V. RESULTS FOR GRAPHENE MONOLAYERS ON SIO2 SUBSTRATES
Theoretical expressions derived in Sec.II (and in Appendix A) are quite general, i.e. are valid for any pair of crystal slabs described by their response functions, while the corresponding surface excitation functions derived in Sec.III are valid for any 2D adsorbed monolayer on any dielectric substrate. In this section we shall apply these results to calculate reactive and dissipative response of various combinations of slabs consisting of graphene monolayers with variable doping on SiO 2 substrate, using the dynamical surface response functions of these materials given in Sec.III. Before proceeding with detailed calculations a few general comments are in order. Though the derived expressions for van der Waals and dissipated power (3) and (12), respectively, include temperature dependence, in the systems studied here inclusion of finite temperature leads to practically no effects, therefore all results will be reported for T = 0. The dependence of these two physical properties on the two parameters, the distance between the slabs a and the oscillation amplitude ρ 0 , can be analyzed if we recognize in the expressions (3) and (12) the function which is possible because of the assumed isotropy of graphene response. The function f m (x) is shown in Fig.3 for first four m's, where x = Qρ 0 . Another important factor in (3) and (12) is e −2Qa which defines the cutoff wave vector Q c , depending on the slab separation a. The separations we shall consider in this calculation are a = 10 − 50nm which defines the cutoff wave vector Q c ≈ 0.05a.u.. On the other hand, the ampitudes which will be considered are ρ 0 ≈ 0.1 − 1nm. This finally provides the maximum argument Q of the functions (39) which is x cut ≈ 1. From Fig.3 is obvious that up to x cut only the m = 0 and m = 1 terms will contribute. Moreover, for x < x cut the Bessels functions can be approximated as J 0 ≈ 1 − x 2 4 and J m (x) ≈ x m 2 m m! ; m > 1 and therefore In Fig.3 we see that approximation (40) is valid almost up to x cut .

A. Spectra of coupled modes
In this section we shall first discuss the spectra of coupled plasmon/phonon excitations in one and two graphene/SiO 2 slabs separated by distance a in order to understand the dominant dissipation mechanisms. Fig.4(a) shows the spectrum of surface excitations S(Q, ω) = −ImD(Q, ω) in graphene(200meV)/SiO 2 slab (as shown in Fig.2) and Fig.4(b) in the system which consists of two graphene/SiO 2 slabs (as shown in Fig.1) separated by distance a = 5nm. In the lonwavelength limit the SiO 2 surface suports two surface polar (FK) TO phonons with flat dispersions and the doped graphene contains a Dirac plasmon with square root dispersion. Coupling between these modes results in three branches, as shown in Fig.4(a). For larger Q the first and second flat branches are phononlike, i.e. their induced electrical fields mostly come from polarization modes on the dielectric surface. On the other hand, the third square root branch is plasmon-like, i.e. its induced electrical field mostly comes from charge density oscillations localised in the graphene layer. However, in the Q → 0 limit the strong hybridization (avoided crossings) between these modes occur and they possess mixed plasmon-phonon character. When another slab is brought in the vicinity the three modes in each slab interact which results in the mode splitting and formation of six coupled modes as shown in Fig.4(b). Figure 4(c) shows the spectrum of surface excitations in the graphene(0meV)/SiO 2 slab. Because the pristine graphene does not support Dirac plasmon the spectrum consist just of two weak phonon branches ω T O1 and ω T O2 damped by π → π * excitations. The spectrum of surface excitations in two equal graphene(0meV)/SiO 2 slabs separated by 5nm (not shown here) is very similar to the one shown in Fig.4(c) which indicates weak interaction between phonons in the two slabs. This could be the consequence of strong screening of FK phonons by graphene adlayers which reduces the range of their induced electrical field. Figure 4(d) shows the spectrum in the system which consists of two different slabs, graphene(0meV)/SiO 2 and graphene(200meV)/SiO 2 , separated by 5nm. One can notice interesting hybridization between the Dirac plasmon and two phonons in one slab and two phonons in another slab giving five branches. In the next section we shall explore how particular plasmon-phonon modes contribute to the dissipated power in two oscillating slabs.

B. Modification of van der Waals force
Van der Waals energy and attractive force are usually calculated and measured for static objects. Here we show how their relative oscillating motion can reduce this attraction, which can be relevant not only from the theoretical standpoint but also in some experimental situations and applications. This phenomenon is present also in the case of parallel motion, as shown in the Appendix A, but this situation would be more difficult to realize in practice.
Making use of the approximation (40) for the lowest order terms of the functions f 0 and f 1 given by (39) we can rewrite the expression (3) for the van der Waals energy as where A is given by (A.17) and (A.18). In the T → 0 limit and neglecting higher order terms A reduces to We see that for ρ 0 → 0 the van der Waals energy reduces to the standard result for the static case, and for ρ 0 = 0 and ω 0 = 0 the lowest order corrections scale with ρ 2 0 . From (41), and also from (42), we see that the slab separation a (because of exponential factor e −2Qa ) reduces the wave vector range to Q < 1/2a. Fig.5 shows van der Waals energies E c of two variously doped, unsupported full conductivity (32-35) graphenes as functions of the driving frequency ω 0 . The driving amplitude is ρ 0 = 20nm and separation between slabs is a = 10nm. For the case of two heavily and equally doped graphenes 1 − 1eV (thick black solid line) the 'static' (ω 0 = 0) van der Waals energy is the largest in comparison with other doping combinations. This is reasonable considering that then except of π and π + σ plasmons (and corresponding electron-hole excitations) the graphenes support strong Dirac plasmons which are all in resonance. Therefore, the charge density fluctuation in one slab ImD 1 (ω) resonantly induces electrical field in another slab ReD 2 (ω) to which it couples, and vice versa. As the driving frequency ω 0 increases the fluctuation and the induced field do not match any more, i.e. ImD 1 (ω) and ReD 2 (ω + nω 0 ) become Doppler shifted and vdW energy is expected to decrease. However, the vdW energy first exhibits a wide plateau until ω 0 < 50THz. We performed a separate vdW energy calculation for two unsupported Drude (32,33) graphenes (not shown here) and noticed that it shows the same features as presented in Fig.5. This suggests that Dirac plasmons are responsible for all characteristic features in vdW energy (for larger dopings). Therefore, the plateau arises probably because the Dirac plasmon fluctuation in one slab, e.g. at ω p , can be efficiently screened by induced plasmon field in another slab which is not necessarily at the same frequency ω p . Moreover, graphene, regardless of doping, exhibits perfect screening ReD(Q ≈ 0, ω ≈ 0) ≈ −1 [39] causing that the static point charge feels image potential. This causes that E c shows almost identical plateau for the case of differently doped graphenes 1 − 0.2eV (black solid line) and 1 − 0eV (thin black solid line). As the doping difference increases plateau energy decreases which is reasonable because of plasmon resonance breakdown. For larger ω 0 > 50THz the Dirac plasmon in one slab does not match any more the perfect screening regime in another one, resulting in a rapid decrease or weakening of vdW energy. In the case of weakly doped graphenes, such as the combinations 0.2 − 0.2eV (red dashed line) and 0.2 − 0eV (thin red dashed lines), the 'static' ω 0 ≈ 0 van der Waals energy reduces in comparison with the heavy doping (combinations with 1eV) cases. This is rea-sonable considering that Dirac plasmon spectral weight decreases with doping. Additionally, it can be noted that for lower doping the vdW plateau shifts to ω 0 < 25THz. This is because the perfect screening frequency region can be roughly estimated as ReD(ω < ω p ) ≈ −1, so, as the plasmon energy decreases the frequency interval whithin which fluctuations are perfectly screened becomes narrower. It is interesting to notice that for some frequencies (e.g. ω 0 > 100THz) the resonant but low doping vdW energy (e.g. 0.2 − 0.2eV case) overcomes the heavily doped but off resonance vdW energy (such as the cases 1 − 0.2eV and 1 − 0eV). The static ω 0 = 0 vdW energy of pristine graphenes 0 − 0eV (blue dashed dotted line) is the weakest and shows smooth decreasing, almost linear behaviour. In this case there are no Dirac plasmons in the graphenes spectra. Therefore, only resonant coupling between π → π * electron-hole excitations, π and π +σ plasmons contribute to the vdW energy. As the frequency ω 0 increases the overlap between these electronic excitations decreases causing smooth and linear vdW energy weakening. The same linear behaviour (for ω 0 > 50THz) can be noticed for doping combinations 0.2−0.2eV and 0.2−0eV which proves that for lower dopings the dominant vdW energy weakening mechanism becomes off-resonant coupling between π → π * electron-hole excitations, π and π + σ plasmons.
It should be noted here that such designed (graphene based) slabs might enable modification of attraction between slabs, e.g. controlled 'sticking' and 'un-sticking' of two slabs. For example, two heavily doped graphenes (1 − 1eV case in Fig.5) are strongly bound, however binding energy between pristine graphenes (0 − 0eV case achieved, e.g. simply by electrostatic gating) is reduced more than twice. Moreover, for larger ω 0 (and fixed doping) the dynamical binding energy is substantially reduced, leading to 'un-sticking' of two slabs, and vice versa, their 're-sticking' by reducing the driving frequency.

C. Dissipated power -substrate dependence
In this section we shall explore how the dissipation power in two oscillating slabs depends on the conductivity model we use to describe graphene and how substrate influences the dissipation power.
In order to facilitate the analysis of the results we shall again use the approximation (40). The lowest order term which contributes in (12) is f 1 , and from Fig.3 it is obvious that, for x < x cut , the higher order terms (m = 2, 3, ...) do not contribute and f 1 can be freely approximated by (40) (red dotted line). Furthermore, because the higher order processes (see Fig.9) included in (12) weakly influence the power P it can be calculated using equation (11) which includes only the lowest order process. Therefore the formula for the dissipated power can be rewritten as This suggests that the dissipated power, within the parameter space used in this investigation (for x < x cut ), behaves as P ∼ ρ 2 0 . Also Eq.42 suggests that the resonant condition (maximum in P ) will occur when the driving frequencies satisfy the condition ω 0 = n 1 ω i + n 2 ω j ; n 1 , n 2 = 1, 2, 3, ...
where ω i = ω p , ω T O1 and ω T O2 are the frequencies of hybridized Dirac plasmons and TO phonons, respectively. Figs 6 show the dissipated power P (ω 0 ) for two oscillating graphene monolayers, calculated in several approximations: unsupported graphene (no substrate) using Drude expression (32)(33) for the conductivity (blue thin line), and using full conductivity (32)(33)35) (red dashed line), as well as for graphenes on semiinfinite (∆ → ∞) SiO 2 substrates with full expression for conductivity (black solid line). Both graphene monolayers are doped so that E F 1 = E F 2 = 200meV. In Fig.6(a) the separation between slabs and oscillation amplitude are a = 10nm and ρ 0 = 0.1nm, respectively. We see that in the Drude model P shows a strong peak which comes from the excitation of undamped Dirac plasmons.
In the full conductivity model plasmon peak is strongly suppressed and interband π → π * excitations become the dominant dissipation mechanism. The fingerprint of π → π * excitations in Fig.6(a) is linear P (ω 0 ) behaviour starting at about 200THz, where we also added cyan dashed lines to guide the eye. It can also be noted that the plasmon peak is red shifted which is reasonable considering that π → π * transitions push Dirac plasmon dispersion toward lower energies.
In the presence of the substrate dissipation is additionally reduced by almost a factor of three. This is because for smaller separations (a = 10nm) the modes with higher wave vectors (e.g. Q ≈ 0.01a.u.), which is in this case only the Dirac plasmon, dominantly contribute to P . In this wavevector region the Dirac plasmon already has high enough frequency (ω ≈ 60THz) that the dynamical part of the substrate screening in not active and the substrate dielectric function can be approximated by S (ω) ≈ ∞ . This causes the reduction of substrate screened Coulomb interactionṽ Q (ω) = 2 1+ S (ω) v Q (see Eq.18) and then (considering Eq.23) reduction of the plasmon intensity, which finally causes the reduction of P . Reduction of the screened Coulomb interaction (18) also causes the reduction of the plasmon frequency which can also be noted. Fig.6(b) shows the dissipated power P (ω 0 ) for the same set of parameters as in Fig.6(a) except that the separation between slabs is increased to a = 50nm. As expected, from the discussion in Sec.V B, P is reduced by about four orders of magnitude and plasmon peaks are shifted toward lower frequencies. The latter is also expected considering that for larger separations the modes with smaller Q contribute, and here the Dirac plasmon has lower energy. We can notice qualitative difference between P in Figs.6(a) and (b) for the case when substrate is present (black lines). In Fig.6(b) P possesses additional structures (two additional peaks at ω T O1 + ω T O2 and 2ω T O2 ) which are not present in Fig.6(a). This is because for larger a the modes with smaller wave vectors (e.g. Q ≈ 0.002a.u.) start contributing to P , and this is exactly the region where plasmon/phonon hybridization occurs (as ilustrated in Fig.4(a)), so the additional peaks at ω T O1 + ω T O2 and 2ω T O2 represent the resonant dissipation to two phonon modes.
Figs.6(c) and (d) show the dissipated power P for the same parameters as in Figs.6(a) and (b), respectively, except that the oscillation amplitude is increased to ρ 0 = 1nm. P in Figs.6(c) and (d) are qualitatively the same and exactly hundred times larger than P in Figs.6(a) and (b). This confirms P ∼ ρ 2 0 behaviour of the dissipated power with amplitude as predicted by Eq.42.

D. Dissipated power -graphene doping and distance dependence
In this section we shall explore the dissipated power for two oscillationg slabs for different graphene dopings. If both graphenes are doped P shows the plasmon peak at about 2ω p = 100THz, and starting at about 200THz it increases linearly, which is the consequence of interband π → π * excitations, as already observed in Fig.6. However, if one doped graphene is replaced by pristine graphene (E F = 0), which does not support the Dirac plasmon (as shown in Fig.4(c)), the Dirac plasmon in doped graphene can no longer resonantly transfer energy to the Dirac plasmon in another graphene and P loses the plasmon peak at 2ω p . However, the visible step remains (at about 75THz) which is the consequence of energy transfer between Dirac plasmon in the doped graphene and π → π * excitations in the undoped one. In this case (small a and larger Q) phonons are still very weak and do not represent important dissipation channel. When both graphenes are pristine the only dissipation comes from the resonant energy transfer between π → π * excitations in different graphenes, resulting in the strictly linear behaviour of P . Fig.7(b) shows the dissipated power P for the same parameters as in Fig.7(a) except that the separation between slabs is increased to a = 50nm. As we have already discussed in Fig.6(a), in this case the modes with smaller wave vectors Q contribute and the dissipated power P gets additional structures coming from resonant phonon excitations. For the case E F 1 − E F 2 = 200 − 200meV (black solid line) (coupling between modes in Fig.4a) the dissipated power shows three peaks at ω T O1 + ω T O2 ≈ 40THz, 2ω T O2 ≈ 60THz and 2ω p ≈ 75THz. For there is a possibility for resonant coupling between two phonons in the slab with pristine graphene and three hybridized plasmon/phonon modes in the slab with doped graphene (coupling beteen modes in Fig.4a and modes in Fig.4c). The three peaks correspond to resonant couplings at ω T O1 + ω T O2 , 2ω T O2 and ω T O2 + ω p , as denoted in Fig.7(b). When both graphenes are pristine, i.e. E F 1 − E F 2 = 0 − 0meV (thin solid blue line) the dominant dissipation channels become the resonant coupling between phonons in both slabs (coupling between modes in Fig.4(c)). The three peaks correspond to resonant couplings at 2ω T O1 , ω T O1 + ω T O2 and 2ω T O2 , as denoted in Fig.7(b). Of course, in all three cases P shows linear behaviour for larger ω 0 coming from the resonant π → π * excitations in both slabs. Figs.7(c) and (d) show the same as Figs.7 (a) and (b), except that the oscillation amplitude is increased to ρ 0 = 1nm. As in Figs.6, P is qualitatively similar and exactly hundred times larger than P in Figs.7(a) and (b). This again confirms the P ∼ ρ 2 0 behaviour. This strong dependence of dissipated power on graphene doping suggests many opportunities for applications.

VI. CONCLUSIONS
In this paper we have provided a complete theoretical description of van der Waals and friction forces for two slabs in relative oscillatory motion which includes variable temperatures in two slabs, their dynamical properties, and dependence on characteristic oscillation amplitude and frequency. In Appendix we also provide, for comparison, analogous expressions for the slabs in parallel uniform motion.
We applied this formulation to explore van der Waals and friction forces between two oscillating slabs, each consisting of atomically thick crystal (e.g. graphene) adsorbed on a dielectric substrate (SiO 2 ). We explore dependence of these forces on osillator characteristics such as driving amplitude ρ 0 and frequeny ω 0 , but also on slab separation a, on graphene doping E F and on substrate properties. We show how the spectra of coupled electronic/phononic excitations in slabs determine the energy transfer processes in this system.
We show that, in general, as the driving frequency ω 0 increases the vdW energy first shows an unusual plateau, and then decreases. We propose the idea of controlling the 'sticking' and 'un-sticking' of slabs by tuning the graphene dopings E F i and driving frequency ω 0 .
We also found a simple ρ 2 0 dependence of both the vdW force and dissipated power. The dissipated power between Drude model graphenes, as function of ω 0 , shows unrealistically strong 2ω p peak. However, in a realistic graphene (whose dielectric properties are calculated from first principles) this peak is strongly reduced and red shifted. We also explain why the substrate substantially reduces dissipated power peak 2ω p . For larger separations a additional peaks appears in dissipation power originating from the excitations of hybridized substrate phonons.
We showed that if one graphene is pristine (E F = 0) it causes the disappearance of the strong 2ω p peak in the dissipated power. Moreover, for larger separations a the doping causes shifts, appearance and disappearance of many peaks originating from resonant coupling between hybridised electronic/phononic excitations in   graphene/substrate slabs.
Obviously, when present, the Dirac plasmons are the dominant channels through which the energy between slabs can be transferred, so the studied model system strongly supports the possibility to control the energy or heat transfer between the slabs by tuning the graphene doping, e.g. by electrostatic gating. More radically, for zero doping E F = 0 the energy transfer can be locked, and vice versa.
In conclusion, it is expected that studies of energy transfer processes in the case of osillating slabs, based on our complete theoretical description, will provide supplementary and more practical approach as compared to those in parallel uniform motion. We shall first derive the van der Waals potential and force between two inequivalent slabs, described by their response functions R 1 and R 2 , moving with relative parallel velocity v and separated by a, as can be seen in Fig.8. In the following we shall briefly summarize the derivation presented in Ref. [19], modified to describe the most general case, i.e. for the slabs with different response functions R 1 = R 2 and different temperatures T 1 = T 2 , including the case of graphene monolayers deposited on dielectric substrates. In the diagram in Fig Here V is the Coulomb potential, S 1 is the correlation function of the left slab and R 2 is the response function of the right slab. We assume that the slab 1 is moving with velocity v so that the parallel coordinates in S 1 are transformed as If we use translational invariance in time and in the parallel direction and perform the Fourier transform in parallel coordinates we find The Fourier transform in time gives: where we have introduced ω = ω + Qv. Because the charge densities in slabs 1 and 2 do not overlap, z integrals in (A.4) contribute only for z 3 > z and z 2 > z 1 , so that we can write where v Q = 2πe 2 Q . Also, if we use the definition of the surface correlation function (A.6) and the definition of the surface excitation propagator [40,41]  Moreover, after we use the connection between the thermal/quantum mechanical charge density fluctuations and the dissipation in the left slab: where n 1 (ω) = 1/(e β1 ω −1) represents the Bose-Einstein distribution, β = k B T 1 and T 1 is the temperture of slab 1, the expression (A.8) becomes ImD 1 (Q, ω)ReD 2 (Q, ω ). (A.10) Here we have used the fact that ImD 2 (Q, ω) is an odd function of ω and does not contribute in (A.8). To this we have to add the contribution from the process in which the charge density fluctuation is created in the slab 2.
Because then slab 2 moves with parallel velocity v relative to slab 1 this contribution can be obtained from (A.10) by exchanging v → −v and 1 ↔ 2, and the result for the van der Waals energy is: E c given by (A.11) includes only the lowest order processes shown in Fig.1a. If we want to include higher order processes shown in Fig.9, we have to replace the interaction v Q which appears in D 1 : (A.12) and the one which appears in (A.13) and integrate over the coupling constant λ to find Notice that (A.14) does not change for v → −v. In order to do the λ integration we transform this expression into: which finally gives the van der Waals energy in the case of unequal slabs and finite velocity: where A(Q, ω, ω ) = [2n 1 (ω) + 1]A 12 (Q, ω, ω ) + (1 ↔ 2) (A.17) and From the van der Waals potential E c (a) we can derive the perpendicular attractive force F ⊥ (a) between two moving slabs: We note that for v = 0 our results agree with the previous ones, but for v = 0 they differ from those in Ref. [1,14]. The functions A and B will also appear in the same form in the expressions for van der Waals potential and force between the oscillating slabs, but with the different choice for ω . We can verify, using spectral representations for ReD's, that our results correspond exactly to the well known result for the van der Waals attraction between two moving or oscillating objects in the lowest order [21], e.g. for T = 0: where ω = ω + ∆ω and ∆ω = Qv for uniform motion or ∆ω = nω 0 for an oscillator. The frequency integral can be rewritten as: which is exactly the lowest order term in (A.16).

Dissipated power and friction force
Now we shall calculate the energy dissipated by the two slabs in parallel uniform motion following the derivation in Ref. [19].
Suppose that the left slab is moving parallel to the right one with relative velocity v and that a charge density fluctuation is spontaneously created in the left slab at the moment t 1 (Fig.1). Propagating in time between t 1 and t it induces charge density fluctuations in the right slab with which it can subsequently interact. In such a process the left slab can be considered as a source which is transferring energy to the right slab, and in analogy with Eqs.3 and 4 of Ref. [19], the energy loss rate operator in this process can be written aŝ where D 2 is the retarded response function of the right slab andρ(r, t) andρ(r, t) are density operators which represent quantum mechanical charge density fluctuations created and annihilated at points (r 1 , t 1 ) and (r, t), respectively. Energy transfer rate from the left to the right slab can be obtained by taking the ground state matrix element of Eq.(A.28) where S 1 (r, r 1 , t, t 1 ) = ρ(r, t)ρ(r 1 , t 1 ) + ρ(r 1 , t 1 )ρ(r, t) (A.30) is the correlation function of the left slab which represents real charge density fluctuation. Eq. (A.29) can be illustrated by the the Feynman diagram in Fig.10. We note that in the inertial system of the right slab the charge density in the left slab, apart from the fluctuations, has an additional parallel component of motion, so all parallel coordinates in the left slab have to be transformed as in (A.2). Explicitly, the correlation function (A.30) becomes S 1 (r, r 1 , t, t 1 ) = S 1 (z, z 1 , ρ − vt, ρ 1 − vt 1 , t, t 1 ). (A.31) After inserting (A.31) into (A.29) and the Fourier transformation in parallel coordinates and in time we get the formula for energy transfer rate per unit surface area from the left to the right slab  Fig.10a. We see that if the charge fluctuation is created with the energy ω it can create excitations in the right slab with the energy ω = ω + vQ. This is expected, namely, ω is the energy in the inertial system of the left slab, but in the inertial system of the right slab it is Doppler shifted by vQ. In (A.35) we have calculated energy transferred from the left to the right slab. However, the part of this energy belongs to the quantum mechanical fluctuation which will be reversibly returned back to the left slab. We can calculate this part of energy which fluctuates between the slabs by going to the inertial system of the left slab and forgeting for the moment the right one. Sitting in the inertial system of the left slab we know that it is in the quantumechanical (and thermodynamical) equilibrium with the environment (in this case with the right slab). So, the energy just fluctuates between the left slab and the environment, i.e. the energy which is given to the environment is exactly equal to the energy which is received from the environment. The energy given to the environment, i.e. to the right slab, can be calculated using exactly the same ideas as before, except that now the right slab is moving with the velocity −v and the left one is at rest. Therefore, following the same procedure (A. 28-A.35) with the response functions of the right slab transformed as D 2 (r, r 1 , t, t 1 ) = D 2 (z, z 1 , ρ + vt, ρ 1 + vt 1 , t, t 1 ) (A. 36) we obtain the energy that is reversibly given to the right slab (A.37) ImD 1 (Q, ω)ImD 2 (Q, ω ).
This means that the energy which is irreversibly given to the right slab or dissipated energy can be obtained by substracting the reversible contribution P 12 from the total energy transfer P 12 [2n 1 (ω) + 1] ImD 1 (Q, ω)ImD 2 (Q, ω ).
(A.40) As in the case of van der Waals energy in Sec.A 1 the higher order terms can be included by replacing v Q 's in D i 's in (A.40) by an infinite series (A.12,A.13), as also shown in Fig.9, so that we get  Obviously, for v → 0 both P and F vanish. This result agrees with Pendry's alternative derivation [1,17]. The above derivation repeats and generalizes some previously well known results [1, 14? ? -19]. We should note that this derivation takes into account not the local but the full microscopically calculated nonlocal response functions R i (Q, ω, z, z ); i = 1, 2. However, its main purpose is to facilitate the derivation of analogous results for the oscillating slabs in Sec.II.