Impact of nuclear vibrations on van der Waals and Casimir interactions at zero and finite temperature

Retarded many-body framework

The RMB expression for the vdW interaction free energy at temperature T among a collection of N_mol disjoint molecules in an arbitrary macroscopic environment is (12)F=kBT∑n=0∞′ln (det (T∞T−1))︸Φ(iξn)(1)whose integrand Φ(iξ), evaluated at the Matsubara frequencies ξn=2πkBTnℏ (with a half-weight on the n = 0 term in the sum, denoted by a prime), depends on the inverse of the scattering T operator T−1=∑lVl−1−Genv, which encodes EM scattering to all orders among the molecules l of susceptibilities V_l mediated by the macroscopic environment whose Maxwell Green’s function is Genv, and T∞−1=∏lTl∞−1, which encodes the scattering properties of each molecule isolated in vacuum via Tl∞−1=Vl−1−G0 in terms of the vacuum Maxwell Green’s function G0; note that all of these quantities depend on iξ, but this is suppressed for notational brevity.

For the insulating and weakly metallic atomistic systems we consider, we describe the molecular susceptibilities via (15, 16)Vl=−ξ2c2∑p,i,q,jαpi,qj∣fpi〉〈fqj∣(2)in terms of localized basis functions ∣f_pi〉 representing the EM response of valence electrons via dipolar ground-state wave functions of effective polarizabilities α_pi,qj. The polarizability matrix α_pi,qj is constructed from a partitioning of the ground-state electron density of the microscopic system, computed from ab initio density functional theory (DFT) methods including short-range quantum exchange and correlation effects, into atomic fragments, which are then mapped onto a set of coupled valence electronic and nuclear oscillators for each atom (Fig. 1) (see Methods for more details). Our inclusion of nuclear oscillations contrasts with past iterations of the RMB framework, which only captured the high-frequency behavior of electronic oscillations, and leads to the emergence of collective vibrations (phonons) at infrared and lower frequencies, in turn delocalizing the susceptibility. The impact of this nonlocality on EM interactions, especially at short distances, is captured through our use of atom-centered Gaussian basis functions (3, 15)fpi(x)=(2πσp)−3exp[−(x−rp)22σp2]ei(3)centered at the equilibrium atomic positions r_p, polarized along the Cartesian direction e_i, and normalized so that ∣ ∫ f_pi(x) d³x ∣ = 1. The isotropic widths σ_p of these basis functions are chosen at each frequency asσp(iξ)=14π(∣αp(iξ)∣3)1/3(4)such that a dipolar oscillator of isotropic polarizability α_p has the same self-interaction energy −ξ23c2∑i〈fpi∣G0fpi〉 as that of a Gaussian dipole distribution in vacuum. The effective atomic polarizabilities α_p are constructed from the overall molecular polarizability matrix α_pi,qj viaαp(iξ)=13∑q,jαpj,qj(iξ)(5)so that they and, in turn, the Gaussian basis functions capture the spatial extent of nonlocality in the molecular response due to phonons at low frequencies.

• Download high-res image
• Open in new tab
• Download Powerpoint

We demonstrate that the combination of nonlocal response due to phonons, leading to the smearing and renormalization of divergent short-range EM interactions, and the interactions of those phonons with long-range EM fields to yield phonon polaritons strongly modifies vdW free energies at zero and room temperatures (T = 0 and 300 K, respectively), as captured in the power law dependence of the free energy on the separation of each microscopic system considered from the macroscopic gold surface, as well as the ratios of free energies between zero and room temperatures. In particular, we define the power law as followsPower law=∂ln (∣F∣)∂ln(z)(6)

While we primarily consider low-dimensional carbon allotropes, we find that these modifications depend strongly on molecular geometry and dissipation. For the cases of fullerene and carbyne, we compare the results in our RMB framework with phonons to those in the RMB framework neglecting phonons, which changes the form of the susceptibility to be local, as encoded in the polarizability matrix α_pi,qj being diagonal in the atom-centered basis; the corresponding free energy is denoted Fe. In addition, for the case of carbyne, we compare with the Casimir-Polder (CP) approximate free energyFCP=−kBT∑n=0∞′Tr[α↔CP(iξn)⋅G↔sca(iξn,r,r)](7)defined in terms of the scattering Green’s function Gsca=Genv−G0 and an effective polarizability tensor αijCP=∑p,q〈fpi∣Tl∞fqj〉 contracting the molecule to a point dipole such that it includes long-range interactions among molecular degrees of freedom in vacuum to infinite order in the scattering, but only to lowest order with the surface. Last, for the case of graphene, we compare to predictions obtained from common macroscopic (continuum) models.

Zero-dimensional fullerene

We begin with the case of a zero-dimensional fullerene above the gold surface (Fig. 2). An isolated fullerene will not have the same dissipation mechanisms as a fullerene in solution (6), so we neglect the dissipation altogether; however, we constrain the center of mass of this isolated molecule by fixing the positions of two nuclei on opposite sides of the fullerene. At large z, the finite size of the fullerene is negligible, so it may be treated like a point dipole with respect to scattering from the plate. Without phonons, the susceptibility is characterized by a single frequency scale arising from the electronic response (Fig. 2A), so when the cutoff frequency c/z falls below that as z increases (Fig. 2B), the power law monotonically approaches the standard retarded dipolar limit of −4 at zero temperature (Fig. 2C). At room temperature, the power law will eventually increase to −3 only when z becomes comparable to the thermal wavelength, ℏc/(k_BT) ≈ 7.6 μm, because once the cutoff frequency c/z falls below the first Matsubara frequency 2πk_BT/ℏ, only the zero Matsubara frequency will contribute. These predictions are similar to predictions from macroscopic formulations of Casimir physics (6); namely, in the absence of phonons, the free energies at zero versus room temperatures are essentially identical for z ≤ 1 μm, like in typical macroscopic situations.

• Download high-res image
• Open in new tab
• Download Powerpoint

Matters change drastically when phonons are considered, in which case the susceptibility is characterized by two frequency scales due to the vastly different nuclear and electronic masses. Even for z ≤ 1 μm, as z increases, these frequency scales compete with the cutoff frequency c/z to produce pronounced nonmonotonic interaction power laws for both zero and room temperatures; the onset of this deviation based on temperature occurs at separations far smaller than one would expect from common macroscopic predictions. This substantial sensitivity to temperature at small z ≈ 100 nm illustrates that even for a small, compact molecule like fullerene, the interplay between phononic response and long-range EM fields can make the interaction power laws deviate substantially from typical macroscopic predictions. The increased relative importance of phononic response at large z also leads to strong deviations of the free energy ratios F(T)/F_e(0) from 1 for large z at both temperatures (Fig. 2C). At smaller separations, the finite size and curved spherical shape of the fullerene dominate the interaction power law. Because of the small size of the fullerene, at larger frequencies ∼ c/z, phonons neither substantially delocalize the molecular response nor strongly couple to EM fields, consistent with the fact that Gaussian widths throughout the molecule are smaller than the smallest value of z = 1 nm considered here. Hence, all three power laws converge upon each other for z ∈ [1 nm,10 nm], and the free energy ratios also converge to 1 in this limit. Note that at z ≈ 100 nm, where the power laws exhibit nonmonotonic behavior, the vertical vdW force on the fullerene by the surface is 1.38 × 10⁻¹⁸ N at room temperature with phonons, which is far smaller than currently measurable in state-of-the-art vdW or Casimir force experiments (14, 17, 18); likewise, the magnitude of the interaction free energy at that separation is less than 10⁻²⁵ J [see the Supplementary Materials for more details, which refers to (19–31)]. That said, if these forces were measurable, the force at that separation would be 15% larger than the prediction of 1.20 × 10⁻¹⁸ N for the room temperature force without phonons.

One-dimensional carbyne

We now turn to the case of a one-dimensional carbyne wire (Fig. 3). As with the fullerene, we assume dissipation to be negligible when the molecule is in isolation and constrain the nuclei at each end of the wire to remain fixed. While the qualitative behaviors of the interaction power laws ∂ ln (F)/∂ ln (z) for a wire at large z above a gold plate are similar to those of the fullerene, the elongated shape of the wire allows it to support longer-wavelength phonons, which couple much more strongly to low-frequency and infrared EM fields, leading to a richer dependence on separation and temperature even for z < 100 nm. To better understand the physics of these interactions, we compare the RMB predictions to those from the CP approximation of (7) for a 500-atom-long carbyne wire above a gold plate.

• Download high-res image
• Open in new tab
• Download Powerpoint

For a 500-atom-long carbyne wire, delocalization in the response due to phonons has the strongest effect at low frequencies, leading to static Gaussian basis function widths σ₀(ω → 0) ≈ 3.3 nm at the middle of the wire. Figure 3A plots α₀ ∼ σ₀(iξ)³ as a function of ξ, showing the existence of two characteristic frequency scales arising from the much stronger phonon-induced electronic delocalization at infrared wavelengths. In contrast, the response in the absence of phonons exhibits only a single frequency scale, and the Gaussian widths never exceed 1 Å. These enlarged Gaussian smearing widths lead to nonmonotonicity in the RMB integrand (Fig. 3B) with respect to ξ for z < σ₀(0), as the Gaussian basis functions overlap with the response of the gold plate (i.e., interactions with image basis functions in the perfectly conducting limit). This nonmonotonicity cannot be observed in the absence of phonons, as the bare susceptibility is essentially local in that case. However, as this nonmonotonicity occurs for very small ξ in the integrand, the behavior of the RMB power law at zero temperature with phonons with respect to z is less sensitive to the nonmonotonic integrand, simply approaching the power law without phonons as z decreases (Fig. 3C). By contrast, at room temperature, the RMB power law with phonons shows substantial deviations from that at zero temperature even up to z < 20 nm; essentially, the sampling of Matsubara frequencies at room temperature makes the free energy disproportionately sensitive to the response at static and infrared ξ. In addition, the sensitivity of these vdW interactions at room temperature to the response at infrared and smaller ξ leads to a free energy ratio F(T)/F(0) that exceeds 2 even for z < 20 nm, and that energy ratio is nonmonotonic because the zero and room temperature RMB free energy power laws with phonons cross each other. At large z, the magnitude of the RMB free energy ratio F(300 K)/F(0) is consistent with recent work by Maghrebi et al. (32), which shows the sensitivity of this energy ratio to geometry for continuum objects even in the absence of nonlocal response. Note that at z = 4 nm, the vertical force on the wire at room temperature is 1.11 × 10⁻¹¹ N with phonons, and as can be observed from the power laws and energy ratios, the ratio of room to zero temperature forces is approximately 1.2; in addition, the room temperature force without phonons is 8.3 × 10⁻¹² N, which is 25% smaller than the prediction with phonons. Both the forces themselves and their differences with respect to temperature should therefore be measurable and resolvable in state-of-the-art Casimir experiments (33), although it should be pointed out that long free-standing carbyne wires have not been stably fabricable, and carbyne has only been found in solution or confined to supramolecular structures like carbon nanotubes (34, 35).

To better understand these phenomena, we compare these results to those obtained by the CP approximation, where the phonon polaritonic response is contracted into a point dipolar polarizability. We find that if the response is associated with a Dirac delta basis function characteristic of an ideal point dipole, the integrand monotonically and markedly drops as ξ increases, so the resulting vdW interaction power laws for z < 20 nm are monotonic at zero and finite temperature, and the energy ratios asymptotically approach a constant larger than 1 as z → 0. However, if the response is associated with a Gaussian basis function whose width σ₀(iξ) is artificially chosen to match that used in the RMB framework, the overlap of the basis function with its reflected image for z < σ₀(0) will lead to nonmonotonic integrands, and thus nonmonotonic power laws at room temperature, although such a “smeared CP approximation” still leads to quantitative differences compared to RMB, as it neglects explicit consideration of the finite molecular size and geometry. In any case, for z > σ₀(0), the CP approximate power laws converge to each other irrespective of basis function choice, and both converge to the corresponding RMB power laws at each respective temperature at much larger z.

To conclude the discussion of carbyne, we point out that nonmonotonic power laws are not merely an artifact of the overlap of enlarged Gaussian basis functions with the gold plate but can be seen even at zero temperature for parallel carbyne wires in vacuum. In addition, we point out that the nonmonotonic behavior in the room temperature power law can be seen as a bump in the corresponding free energy in that range of separations, with the magnitude of the interaction free energy at room temperature exceeding 10⁻²⁰ J there. Further details concerning both of these points are provided in the Supplementary Materials.

Two-dimensional graphene or BN

Last, we consider two-dimensional atomically thin materials above the gold plate, starting with an infinite graphene sheet. Unlike fullerene molecules or carbyne wires, atomically thin graphene sheets can be exfoliated and suspended in vacuum (36, 37); therefore, compared to the other aforementioned carbon allotropes, there is substantially more theoretical and experimental work characterizing the mechanical and vibrational (22, 38), electronic (39), and thermal (40) properties of graphene. In the particular context of vdW interactions, several theoretical macroscopic models for the response of graphene have been used. Given this, we specifically compare our RMB predictions to those obtained from a Lifshitz formula of the Casimir interaction between a graphene and a plate, based on a tight-binding model of the electronic band structure of graphene. Such a model is consistent with the random phase approximation (RPA) (41) and includes spatial dispersion and the possibility of doping but does not consider contributions from phonons or other dissipation mechanisms. In contrast, our RMB model of graphene explicitly includes phonons, and because extended graphene has many more channels for damping of the response, we include explicit nonzero dissipation terms in the susceptibility (see Methods for more details), unlike our neglect of dissipation in fullerene and carbyne.

Casimir and vdW interactions involve sums over all frequencies (1), in which case the infrared response (including both temporal and spatial dispersion) is expected to be relevant. We therefore begin by demonstrating the differences in the graphene nonlocal susceptibility χ(iξ, k) (see Methods for more details) in our RMB framework (Fig. 4A) compared to the RPA model (Fig. 4B) in the absence of doping and at zero temperature. In the RMB model, at small ξ, and especially at small k, corresponding to nearly spatially uniform incident fields, contributions from long-wavelength acoustic phonons tend to dominate, delocalizing the electronic response over several primitive unit cells and drastically increasing ∣χ∣ to more than 100; as ∣k∣ increases, so too do the contributions of localized optical phonons, decreasing the magnitude and nonlocality of the response. It is only for ξ > 10¹⁴ rad/s that phonons no longer contribute appreciably to the response, so ∣χ∣ is essentially independent of k and its magnitude is less than 10, consistent with RPA and decaying quickly as ξ increases further. In contrast to RMB, the electronic orbitals in the RPA model are tightly bound to fixed nuclei, in which case the resulting response vanishes as k → 0 for nonzero frequency, and hence χ increases with increasing ∣k∣, in qualitative opposition to the RMB prediction. In particular, in the short-wavelength regime, χ approaches a constant asymptote qe2g/(32ℏϵ0vF)≈4 at a rate controlled by ξ, which is far smaller than the RMB prediction at low ξ. Given that the π orbitals within RPA are tightly bound to nearby nuclei, one would expect that phonons should have a substantial impact upon the valence electronic and static response of graphene. The low-frequency behavior of metals is particularly important in the context of vdW interactions, as exemplified by recent discrepancies in the response of Drude versus plasma models of gold (33).

• Download high-res image
• Open in new tab
• Download Powerpoint

We now compare the predictions of both RMB and RPA models for the vdW interaction free energy per unit area of graphene above a gold plate, at both zero and room temperatures. In particular, Fig. 4C shows the free energy power law and corresponding energy ratios, with respect to the gap separation z. The RPA model is analyzed in the absence of electronic doping. Our RMB predictions for the power laws are qualitatively similar to those of the fullerene in that the delocalization in the response due to phonons leads to power laws at both temperatures that are nonmonotonic, although the higher dimensionality of graphene compared to fullerene or carbyne causes the onset of nonmonotonic behavior to arise at smaller z, around 20 nm, than for the fullerene or carbyne. Likewise, as the impact of phonons is more pronounced at larger temperatures, the energy ratio starts to deviate noticeably from 1 at z, not much smaller than 20 nm. Note that, at z = 100 nm, the ratio of the vertical forces at room to zero temperatures is approximately 1.3 and the corresponding pressure at zero temperature is approximately 1 N/m², so these differences should be measurable in state-of-the-art experiments involving graphene (42, 43); in addition, the free energy per unit area exceeds 10⁻⁸ J/m² in magnitude at that separation, with more details in the Supplementary Materials. Compared to our previous predictions in carbyne systems, the assumption of much larger dissipation in graphene substantially dampens and smears the impact of long-wavelength acoustic phonons, limiting the phonon mean free path and hence the spatial extent of delocalization in our Gaussian widths σ to a much greater degree. Consequently, we observe energy ratios much closer to 1 and monotonic power laws at small separations z < 20 nm. In contrast, the RPA power law is a constant −3 at zero temperature over all z, while it increases monotonically toward −2 at finite temperatures. In addition, since this RPA model does not include phonons, there is no nonmonotonic behavior in the power laws, and the RPA power laws are consistently above −3 in the range of separations z that we consider, unlike the RMB power laws that drop below −3 for z around 20 nm. However, the RPA free energy ratio is 2, rather than 1, even at z = 1 nm, and increases monotonically to 3.2 by z = 100 nm. This is because, in the absence of doping, electrons fill the band structure exactly up to the Dirac point, so thermal fluctuations at room temperature can have a much larger relative effect on the response. In contrast, the material model entering the RMB framework for all materials we consider, including both DFT calculations and atom-centered partitioning of the density, is formulated at zero temperature, and thermal effects are only considered via the Matsubara summation. A more consistent model of atomistic response would account for the effects of finite temperature on the mechanical structure and phonon dispersion of graphene itself.

Developments before the RMB framework, particularly the many-body dispersion (MBD) framework, have shown strong agreement with experimental measurements and state-of-the-art theoretical predictions of static binding energies, equilibrium structures, and elastic moduli in single and multiple layers of graphene, BN, and other two-dimensional materials (44, 45). This suggests that the RMB framework can accurately capture the contributions of phonons and EM retardation to long-range vdW interaction forces involving large molecules and spatially extended low-dimensional materials, including two-dimensional materials, due to our ab initio consideration of both effects. However, the interplay between phonons and inherently delocalized electrons, like that observed in the semimetallic Dirac cone electronic structure in pristine graphene, may not be fully captured in the material model of single electronic oscillators associated with each atom, and this could have further effects [Dobson type C nonadditivity (46)] on the interaction power laws at nanometric and larger separations, which go beyond the regimes typically considered by the MBD framework for binding energies and structural properties. Having examined semimetallic pristine graphene, we now consider atomically thin hexagonal BN (Fig. 5), which behaves as a polar dielectric and which we expect, in the absence of other experimental evidence, may be more accurately described by the material model of single electronic and nuclear oscillators associated with each atom. Numerical difficulties preclude consideration of an infinite BN sheet with periodic boundary conditions as we did with graphene, so we use supercells of approximate dimensions 11 nm × 10 nm in each case, with atoms at two of the four opposite corners constrained, and only consider z ≤ 10 nm to mitigate sampling of finite size effects. The RMB calculations of vdW interactions show similar power laws for F_e (without phonons) for graphene and BN in the range z ∈ [1 nm,10 nm], which suggests the need for more work to distinguish electronic structure effects between graphene and BN even at low frequency; that said, phonons dominate the response, so when their contributions are included, it is less unexpected that the similar nuclear masses, geometries, and resulting internuclear couplings K_I of the two systems lead to negligible qualitative distinctions between the two systems. In particular, at z = 2 nm, where the nonmonotonic behavior of the room temperature power law is most prominent for graphene as well as BN, the room temperature force per unit area is 2.80 × 10⁶ N/m² with phonons versus 2.57 × 10⁶ N/m² without phonons for BN, compared to 2.86 × 10⁶ N/m² with phonons versus 2.62 × 10⁶ N/m² without phonons for graphene. We emphasize the sensitivity of the RMB predictions to the confluence of material effects, dimensionality, and boundary conditions: Finite size leads not only to quantitative deviations when the separation approaches the supercell size but also to qualitative deviations as the room temperature power law for the graphene supercell exhibits nonmonotonic behavior for z < 10 nm, which is similar to the behavior of the compact carbyne wire but different from infinite pristine graphene.

• Download high-res image
• Open in new tab
• Download Powerpoint

At this initial stage, the RMB framework may do a better job handling insulating two-dimensional materials like BN or graphene with defects (whether in the form of vacancies or foreign atoms or molecules) (47) compared to semimetallic pristine graphene. That said, we expect that further development of the electronic response model used in the RMB framework will allow it to accurately model vdW interactions along the range of insulating to semimetallic low-dimensional materials as well as in electronically doped systems like graphene (which we discuss in the Supplementary Materials), which become more metallic due to the fully delocalized response of free mobile charge carriers. This ultimately requires changing the form of the susceptibility V and the corresponding basis functions, without changing the general formula for the vdW energy in terms of EM Green’s functions and scattering operators, thereby clarifying that the improvements must come from the side of material modeling rather than from modified consideration of long-range EM interactions.