How many-body effects modify the van der Waals interaction between graphene sheets

Undoped graphene (Gr) sheets at low temperatures are known, via Random Phase Approximation (RPA) calculations, to exhibit unusual van der Waals (vdW) forces. Here we show that graphene is the first known system where effects beyond the RPA make qualitative changes to the vdW force. For large separations, $D \gtrsim 10$nm where only the $\pi_z$ vdW forces remain, we find the Gr-Gr vdW interaction is substantially reduced from the RPA prediction. Its $D$ dependence is very sensitive to the form of the long-wavelength many-body enhancement of the velocity of the massless Dirac fermions, and may provide independent confirmation of the latter via direct force measurements.


I. INTRODUCTION
It is well known that a zero-gap conical π z electronic Bloch band structure of undoped graphene, supporting massless Dirac fermions propagating with speed v, should give this system a number of unusual properties. [1] One such property relates to the low-temperature dispersion (van der Waals) interaction energy per unit area E vdW /A between undoped parallel graphene sheets separated by a large distance D. Commonly used theories such as the summation of pairwise atomic contributions, E vdW ≈ i =j C ij R −6 ij , and other popular and largely successful approaches [2][3][4][5][6], predict the energy for this case (and any case with parallel 2D sheet geometry) to be a power law where C 4 is a system-dependent constant (see Sections 4 and 8 of Ref. 7 for further discussion). By contrast, more microscopic/collective approaches, such as the RPA correlation energy based on the graphene π z -π * z electronic response, yield the result [8][9][10][11][12] The constant C 3 is easily calculated within the random phase approximation (RPA), which treats the electrons in each layer as essentially non-interacting, but subjected to their own time-dependent classical electrostatic field. Indeed, if one makes the simplest ("Casimir-Polder") approximation, in which the interlayer Coulomb interaction is treated by second-order perturbation theory, one finds, in the limit of large separation where α = e 2 v is the effective fine structure constant of graphene, related to the velocity v of the massless Dirac fermions as the fundamental fine structure constant is to the speed of light, and F (a) is a smoothly varying function, which is plotted in Fig. 1 (the derivation of this result, as well as the analytic expression for F (a) will be presented below). For v = 10 6 m/s one has α ≃ 2.2, which makes C 3 very close to the "universal" value C uni 3 = e 2 32π . This value is weakly dependent on small variations of v, and remarkably, the interaction −C uni 3 /D 3 is not changed by the inclusion of electromagnetic retardation, even at large D values. (see Eq. 36 of [13], [14], and Eq. 8b of [8] corrected for a spurious factor 2) In this context it is important to note that for real graphene Eq. (2) is valid only at large separations: at shorter distances, gapped transitions other than the π z → π * z contribute a vdW energy of the conventional form (1). It is only for D 10nm that numerical work within the RPA [15] suggests the D −3 falloff overtakes the conventional D −4 contribution. It is therefore in this regime of larger separations corresponding to small wavevectors that one should check for any many-body effects beyond the RPA due to the π z − π * z graphene response. Furthermore at these separations all electrostatic and metallic overlap forces have long vanished.

II. RENORMALIZATION OF THE VELOCITY
Such corrections are worth investigating because the π z electrons in a graphene sheet are obviously not "weakly interacting": the coupling constant α is larger than one [16][17][18][19][20]. What is believed to be true is that the effective long-wavelength Hamiltonian, generated by a renormalization group (RG) flow, is non-interacting. [16,17,[21][22][23] However, the reason for this simplification is not that the electric charge renormalizes to zero but that the fermion velocity grows to infinity. More precisely, if one introduces an effective Fermion velocity v q , which describes the system at length scales larger than 1/q, then v q → ∞ for q → 0, and, accordingly, the running coupling constant α q = e 2 vq tends to zero. This long-wavelength renormalization is quite distinct from beyond-RPA corrections at short wavelength, arising from modified forms of Adiabatic Local Density Functional Theory (ALDA), which are also sometimes also described as "renormalized" [24] Our long-wavelength corrections can change the asymptotic power law for the van der Waals interaction, whereas the short-ranged ALDAbased ones, unsurprisingly, have no major qualitative effect on long-wavelength vdW phenomena [25].
The stronger the original interaction is at the microscopic scale, the larger the renormalized velocity becomes at any given length scale 1/q. The exact form of the divergence of v q for q → 0 is, of course, unknown, since the many-body problem has not been solved. RG calculations based on first-order perturbation theory [16,23] suggest a weak logarithmic divergence of the form where Λ is a microscopic cutoff of the order of the inverse of the lattice constant (Å). More sophisticated calculations at the "two-loop" level and in the large-N limit, [17,23,[26][27][28][29][30] N being the number of Fermion flavors, predict a stronger power-law divergence of the form where β = 8 N π 2 , (N = 4 for graphene). Recent experiments performed by a variety of techniques have at least partially confirmed these theoretical predictions [18][19][20]31], showing that the effective coupling constant is reduced and the Dirac cones are strongly compressed [19] near the Dirac point.
The many-body enhancement of the fermion velocity has also been shown to affect various many-body phenomena, such as the plasmon dispersion [32,33], the optical Drude weight, [32] and the electronic screening of external charges. [34] To date [24,25], beyond-RPA effects were believed to alter at most the prefactor, not the power exponent, of vdW decay with distance. In this Letter we expose a striking case where beyond-RPA many-body renormalization affects the essential character of a vdW interaction, namely that between graphene sheets. Our main result is that the long-wavelength enhancement of the fermion velocity causes the vdW interaction to decrease asymptotically faster than in Eq. (2) but still slower than in the conventional Eq. (1). This result can be expressed in an intuitively appealing way by saying that the bare α in Eq. (3) must be replaced by the running coupling constant α q evaluated at q = 1/D. Thus we have and since F (x) ≃ π 2 x for x ≪ 1 the asymptotic behavior of the vdW interaction is reduced, relative to the RPA, precisely by the factor α 1/D , which vanishes in the limit D → ∞. The renormalized interaction will therefore decrease as [D 3 ln(DΛ)] −1 or as D −(3+β) , depending on which of the two scenarios, (4) or (5), is realized. Additional many-body effects contained in the so-called vertex corrections turn out to be irrelevant at sufficiently large distance, even though, of course, they can quantitatively change the result at intermediate distances. Throughout the analysis, we assume that the distances are not so large that v q=1/(2D) becomes comparable to the speed of light, at which point electromagnetic retardation effects should be taken into account. Using this criterion, retardation becomes dominant only for D of order 10 220 m (!) for the logarithmic renormalization case [Eq. (4)] or O(10 1 m) for the power law case [Eq. (4)]. Thus retardation here is unimportant in practice, as for the pure RPA theory.

III. VAN DER WAALS CALCULATIONS
We consider the interaction energy between two parallel freely suspended graphene sheets in vacuo, separated by a distance D that is much larger than the twodimensional lattice constant a, as well as the thickness T of each sheet (see Appendix A for other scenarios). Our starting point is the Casimir-Polder (CP) formula obtained by doing straightforward second-order perturbation theory in the inter-layer electron-electron interaction potential where q is the two-dimensional wave vector of density fluctuations in each layer. The result is Here χ(q, iu) is the electronic density-density response function of a single isolated layer, evaluated at wave vec-tor q and imaginary frequency iu. The exponentially decaying factor exp(−2qD) ensures that small wave vectors q ≈ 1/D dominate Eq. (8) for large separations. The response function χ(q, iu) is expressed in terms of the proper response functionχ(q, iu) and the intra-layer interaction potential V intra (q) = 2πe 2 q according to the well-known formula [35] In the RPA one approximatesχ(q, iu) ≈ χ 0 (q, iu), where is the non-interacting zero-temperature response function for the conical π z −bands, [8,16,36] v is the bare velocity, and x ≡ u qv . Intra-layer electron-electron interactions modifyχ in two ways: via self-energy insertions and via vertex corrections. For example, in a recent first-order perturbative calculation coupled with RG arguments, Sodemann and Fogler (SF) find [34] where v q is given by Eq. (4) q ) and α q is the corresponding coupling constant. Here the self-energy insertion has caused the bare velocity v to be replaced by the renormalized, scale-dependent velocity v q . The second term, J, is a dimensionless function representing the combined effects of self-energy and vertex corrections beyond the simple rescaling of v. In the notation of SF J is given by J(x) ≡ I a (ix) + I b (ix) where the functions I a and I b , defined in Ref. 34, are analytically continued here to the imaginary x axis. It is essential to our subsequent arguments that J(x) is a smooth function of x, varying monotonically between J(0) ≈ 0.497 and J(x → ∞) ≈ 0.013 x ≡ C∞ x (see Eq. (12) and Fig. 2, where a numerical fit to J(x) is provided). This feature is expected to persist beyond the first-order approximation, e.g., even in the strong coupling limit, where the renormalized velocity is likely given by v where a = 0.285 and y = √ 1 + x 2 . This is shown in Fig. 2. For imaginary x, this numerical fit provides a good match to the results of Ref. 34.
Substituting Eq. (11) into Eq. (9) we obtain the full interacting response function in the form where H(x) ≡ J(x) √ 1 + x 2 . Plugging into Eq. (8) and changing integration variable from q toq = 2qD we find, after simple manipulations where we have defined the functioñ If now self-energy and vertex corrections are neglected by setting H = 0, we see thatF simplifies to which can be evaluated analytically, yielding and F (a) = 2/3 for a = 1. This is precisely the function F that was introduced in Eq. (3) and was plotted in Fig. (1). Since in RPA α q is constant (α q = α) and ∞ 0q 2 e −q dq = 2, Eq. (3) is seen to follow immediately from Eq. (14).
To go beyond the RPA we must reinstate the selfenergy and the vertex corrections. However, we observe that in the large-D limit the coupling constant tends to zero and therefore theF function reduces to the F function, which in turn can be replaced by its small-a expansion F (a) ≈ π 2 a. Thus, Eq. (14) becomes The asymptotic behavior of the vdW interaction depends solely on the behavior of the running coupling constant α q (or, equivalently, the renormalized velocity) in the q → 0 limit. For the logarithmic renormalization case of Eq. (4), we evaluate (18) by freezing the slowly varying α (1) q/2D , evaluating it at the maximizing value q 0 = 2 of the rapidly varying integrand, getting .
This shows a modest, logarithmic reduction of the vdW interaction relative to the RPA result. If, on the other hand, the strong-coupling model of Eq. (5) is adopted, (18) gives an altogether different power-law behavior: where Γ(x) is the gamma function. Notice that, since β < 1, this is still larger than the D −4 dependence expected for insulating 2D layers, and therefore dominates at large separations in real graphene where gapped insulator type transitions also contribute to the response. The above eqs. (18)(19)(20) are valid at asymptotically large separations. At finite separations the more accurate Eqs (14)-(15) must be used. These now depend not only on the fermion velocity -a measurable quantitybut also on the form of the function J(x), which is not directly accessible to experiment and must be calculated by many-body theory (Notice, however, that the imaginary part of the density response function for real frequency is related to the optical absorption spectrum, which is, in principle, measurable, and could be used to calculate the function J). Making use of J(x) calculated in Ref. 34 and fitted as shown in Fig. 2 we find that is an excellent approximation (relative error under 1%) to the integral of Eq. (15). Using this in Eq. (14), together with Eq. (4) or (5) for the velocity, we find that the distance dependence in the intermediate regime is basically D −3 with only a modest further dependence on distance via F π 2 α q=1/D orF π 2 α q=1/D . [37] In Fig. 3 we plot the force and its local exponent [in F ∝ D −p(D) ] for the RPA energy E RPA /A = −e 2 /(32πD 3 )F (α) and the vdW interaction calculated through (14) using (4) and (5). Using the stretched graphite vdW energy formula of Refs. 11 and 38 yields remarkably similar results for the equivalent in bulk graphite (up to a constant). Naturally, the situation differs between the weak-coupling (4) and strongcoupling (5) models of the renormalized velocity. The interaction energy and force is qualitatively different from RPA at large separations, and shows moderate quantitative differences at intermediate separations, less than 25% for D < 100Å. Observation of such deviations from the RPA will provide additional evidence of the many-body renormalization of the fermion velocity. We note that for intermediate values of D (D 20nm), further modifications are required even at the RPA level to account for departures of the graphene electronic bandstructure from a perfect infinite cone [11], or anisotropy [12].
There have been some proposals ( [23] and references therein) that strong coupling could bring excitonic effects leading to a gap. In that case the vdW energy might show the insulating D −4 behavior as in Eq (1). There seems to be little experimental evidence so far that graphene can be an excitonic insulator, however.

IV. EXPERIMENTAL OPTIONS
While our paper in the main concerns the theory of the renormalization of the vdW interaction, it serves as one of several motivating factors for a renewed experimental effort for direct measurement of the vdW forces on high-quality graphene flakes. These measurements would also be needed to resolve existing controversies about graphitic cohesion in general. In this section we briefly discuss the current prospects for such measurements, and we hope thereby to stimulate experimentalists in this direction.
While there have been a number of experiments [39][40][41][42] that have indirectly determined the binding energy of graphene planes, these have produced a wide range of results, and most have relied on questionable theoretical assumptions for their analysis, so that the whole field is somewhat controversial [43]. Purely theoretical estimates have also varied widely [3,13,15,[44][45][46][47][48][49][50], though recent very large and relatively high-level (RPA [15], DMC [48]) calculations are starting to show consistency. These same high-level types of theory also predict the forces between nanostructures such as graphene planes as a function of distance, showing very different results from popular pairwise-additive-force theories, in the asymptotic region of large separations [see Eqs (2), (6)].
The force between graphenes at such asymptotic distances does not appear to have been measured at all so far. In view of this, as well as the aforementioned binding-energy controversies, direct measurement of such forces at all distances would be very desirable. As the distant forces are relatively small, Atomic Force Microscopy or a related NEMS oscillator approach would be the preferred route. For the basic RPA analysis of two cold undoped graphene sheets of lateral dimension 1 micron, separated by 10 nm, the predicted force is a few nanoNewton, well within AFM capabilities, with larger forces for larger flake areas.
In the ideal case analyzed in the main text there are two undoped, freestanding graphene sheets at a temperature below 10K. Single sheets of high-quality graphene have certainly been subjected to various measurements [51] and vdW forces due to a single supported graphene sheet has been seen [52], but force experiments with two freestanding sheets are rare or absent, and will require some effort to ensure a reasonable degree of parallelism.
Perhaps a better short-term prospect is therefore to measure the force between a single freestanding sheet and its own "vdW image" in the surface of a bulk metallic substrate. A very recent related experiment, using a metal grating instead of graphene, and an O(100µm)radius gold sphere as the substrate to avoid parallellism issues, have been highly successful [53]. This experiment also demonstrated the use of co-located fiber optics for precise control of the separation D down to < 200nm. The theory of this geometry will be analyzed in detail elsewhere, but is expected [see equation (A7)] to yield similar predictions to the ones in the main text.
The unavoidable corrugation of the freestanding graphene sheets should not be a problem, as it is known experimentally not to affect electronic properties substantially, and would only contribute a very few percent uncertainty in the separation D between the sheets, at the separations of a few tens to hundreds of nanometers that are relevant to the main text. This uncertainty would not affect the analysis of the force power law proposed in the text, as one would aim for measurements at D values spanning an order of magnitude.
A more subtle problem is the likely existence of metallic n-and p-type "puddles" on the undoped sheets [54]. Provided that the sheets are of sufficient quality that the puddles are disconnected objects of typical spatial extent λ, we predict that they contribute an "insulatorlike" vdW interaction energy varying with distance D as (const)D −4 when D > λ, clearly distinguishable from the lower powers predicted in the text, arising from the un-doped non-puddle areas. Recent experiments [54] suggest that λ = O(1nm), so that force experiments at D = O(10-200nm) would still exhibit undoped-graphene properties.
Another experimental route would be to avoid separating the sheets in their perpendicular direction, but rather to slide them off one another in a surface-parallel direction. This may be easier experimentally because of the recent centrifugation-based preparation of high-quality micron-sized stacks of 2-5 graphene monolayers [55]. These are spatially staggered like a slipped deck of cards, potentially allowing attachment of an AFM tip to the projecting edges of individual sheets. The force during the entire lateral sliding process, out to wide separation into disjunct coplanar sheets, would be measured. The non-contact part of this force can be expected to show effects from the coupling of long-wavelength electronic charge fluctuations and hence renormalization effects in the graphene polarizability, just as for the case of parallel sheets separated by distance D measured perpendicular to the sheets. From the theory point view, the analysis of this geometry is more difficult and has not yet been attempted in detail.
The above considerations suggest that it will be possible to achieve a much improved understanding of graphenic cohesive forces by direct measurement.

V. CONCLUSION
In conclusion, we have shown that the low-temperature dispersion (vdW) interaction between two infinite parallel non-doped graphene sheets is significantly modified by many-body effects beyond the RPA. Not only is the interaction quantitatively reduced, but also its qualitative asymptotic behavior is modified. The main source of the effect is the many-body renormalization of the velocity of the massless Dirac Fermions. This renormalization has been the subject of many recent investigations [16][17][18][19][20]. It is experimentally observed as a deformation of the Dirac cones near the point of contact. Our findings demonstrate that the same renormalization manifests itself in the long-distance behavior of the dispersion forces.
Direct measurement of the asymptotic graphenegraphene vdW interaction could therefore distinguish between much-debated theory models of electrons in graphene -weak renormalization [Eq. (19)], strong renormalization [Eq. (20)] and excitonic insulator [Eq. (1)]. Such experiments in the asymptotic vdW region will be demanding but we estimate (e.g.) a measurable force of order nN (see inset of Figure 3) between micron-sized graphene sheets separated by O(10nm). Observation of the vdW image force in a metal substrate could avoid the need for two graphene sheets, and there are other possibilities too (see Ref. [55] and second last paragraph of Section IV above). Very recent experiments [53] support the general feasibility of our proposals. We estimate that complications due to graphene wrinkling and puddling [54] will not destroy our effect. Indeed the time is ripe for direct force measurements to clarify this and other recent controversies [43] over graphenic cohesion in general.

ACKNOWLEDGMENTS
TG and JFD were supported by Australian Research Council Grant DP1096240. GV was supported by NSF Grant DMR-1104788. JFD appreciates hospitality at the University of Missouri, the Donostia International Physics Centre, and the European Theoretical Spectroscopy Facility.
Appendix A: van der Waals force between graphene and a metal bulk Here we investigate the van der Waals force between a single graphene layer and a metal bulk, which may be more appropriate for likely successful experimental arrangements (see Sec. IV below). We predict that it will obey the same overall power law as the response between two graphene layers, with at most a logarithmic correction (as a function of D). Similarly the dependence on the many-body beyond-RPA effects is expected to be maintained. Here we offer evidence that this is indeed the case.
We make use of the relationship [56] for the van der Waals potential between semi-infinite bulks and layers interacting across a single surface. In the case of a layer of graphene interacting with a bulk metal where χ Gr (q, iu) is the interacting response of the graphene layer and χ Metal is the bulk response of the metal. Here C Metal = ǫ(iu)−1 ǫ(iu)+1 where ǫ ≈ 1 + ω 2 p u 2 is the dielectric function of the metal.
From equations (9) and (13) we see that C Gr ≡ −2πe 2 q χ Gr (q, iu) can be written as C Gr (x, α q ) ≈ παq/2 √ 1+x 2 +παq/2 where x = u/(v q q) and α q = e 2 /( v q ) varies slowly with q. Thus where we also usedq = 2qD. Clearly the form of α q will have a substantial effect on the asymptotic behaviour of the vdW potential, similar to the bigraphene case. Using C Gr < π/2α x+π/2α and setting αq 2D ≈ α 1/D gives where D 0 ∝ e 2 ωp = O(1Å) i.e. the asymptotic interaction has an extra logarithmic term compared to the bigraphene case kα 1/D /D 3 . The prefactor α 1/D ensures that the difference between different renormalization scenarios is maintained.
One can, of course, perform the integral (A5) [or an equivalent expression for (A1)] numerically for a given graphene velocity v q and dielectric frequency ω p . However, the important asymptotic physics are clearer in the bigraphene case tested in the paper, and the deviation caused by the bulk metal is expected to be small at the O(10-100nm) distances we propose investigating.