Surface van der Waals forces in a nutshell Free

Van der Waals interactions are at the heart of chemical physics. Yet, the standard textbook answer on their essential characteristic is the well known inverse sixth power dependence on the distance. This is a largely biased statement towards the London picture of molecular interactions, which treats intermolecular forces as a result of classical electrodynamic fluctuations. At distances of just a few nanometers, molecular interactions develop a different, faster decay that results from purely quantum electrodynamic fluctuations and was first described in the seminal work of Casimir (cf. Refs. 1 and 2 for a historical perspective).

An acclaimed unification of these two complementary views on the range of molecular interactions has been known for a long time.³ The Dzyaloshinskii-Lifshitz-Pitaevskii (DLP) theory for surface forces across a medium is the solution of the full thermal quantum electrodynamic field theory for the forces between two plates across a dielectric continuum. As a result, it generalizes in one single equation the Hamaker theory of additive dispersive forces, the Debye induction potential between polar and polarizable media, and the Keesom interactions between Boltzmann distributed dipoles. It describes the crossover from purely non-retarded dispersive to retarded long–range interactions, it reduces to the Casimir-Polder formula for the retarded force between metallic plates in vacuum, and it provides the famous result of London for the dispersion interaction between two atoms at short distances.

Not surprisingly, the theory has had a profound impact on fundamental physics;² it has motivated a large number of historical experiments^4–9 and retains its theoretical influence in promising new studies up to date.^10,11

Unfortunately, with some exceptions, the general theory of van der Waals interactions is largely ignored in favor of the Hamaker picture of additive dispersive forces.¹² The reason is simple to infer, as the power and generality of this approach comes at the cost of a complicated and lengthy formula that can hardly be interpreted qualitatively, except in special limiting cases after lengthy manipulations and detailed knowledge on materials properties.^13–17 

Here, we show that the theory of surface van der Waals forces can be formulated as one single simple equation that embodies simultaneously the known low and high temperature limits, the crossover from retarded to non retarded interactions, and the far less appreciated long range exponential suppression of retardation effects. The result allows us to interpret easily the main qualitative features of van der Waals forces for arbitrary distances and provides a convenient means to compute quantitatively the results in analytic form.

The surface free energy between two semi-infinite macroscopic bodies, 1 and 2, separated by a third body, m, of thickness L, may be written in terms of the effective Hamaker function, A(L), as $g(L)=−A(L)12πL2$⁠. In practice, A(L) is a constant at small separations only and develops a complicated L dependence that is given in DLP theory as

where $R1m2M(x,n)=Δ1mM(x,n)Δ2mM(x,n)$⁠, $R1m2E(x,n)=Δ1mE(x,n)Δ2mE(x,n)$⁠, while

and $xi2=x2+(ϵi−ϵm)(2ωnL/c)2$⁠. In these equations, the prime next to the sum implies that the n = 0 term has an extra factor of 1/2. The dielectric function, ϵ_i, is evaluated at the imaginary angular frequencies iω_n, where ω_n = ω_Tn are integer multiples of the thermal Matsubara frequency $ωT=2πkBTℏ$⁠. The magnetic permitivities have been assumed equal to unity. The lower integration limit is $rn=2ϵm1/2ωnL/c$⁠, while k_B and ℏ, c are the usual fundamental constants. Equation (1) provides the leading order result of the full DLP theory and is the starting point of most analytical approximations for the Hamaker function.^{4,12,13,15,16,18,19}

It is conventional to split the sum into n = 0 and n > 0 contributions, such that A(L) = A_ω=0 + A_ω>0(L). For the first term n = 0, it is possible to integrate over x exactly, yielding the well known approximation

For the remaining contribution A_ω>0, both $R1m2M(x,n)$ and $R1m2E(x,n)$ remain non trivial functions of x, and the integral cannot be evaluated exactly by analytical means. However, inspired by Parsegian’s insightful monograph,¹⁴ we perform the integration over x using a generalized one-point Gauss-Laguerre quadrature rule with weight xe^−x in the interval [r_n, ∞]. This yields (see the supplementary material)

where we have written $R=R1m2M+R1m2E$ for short, and R(x, n) is evaluated at the optimized value

Next, we transform the sum of Eq. (4) into an integral. Using the Euler-McLaurin formula, we quantify the first order correction to this transformation as $≈34kBT(1+νTL)e−νTL$⁠, which is a small fraction of the full integral in most practical situations. Therefore, introducing the auxiliary variable ν_n, such that r_n = ν_nL, and defining the constant $νT=2ϵm1/2ωT/c$⁠, we find

where $R̃(x*,ν)=ϵm−1/2jm−1(ν)R(x*,ν)$ and $jm=1+12d⁡lnϵmd⁡ln⁡ω$⁠. The factor j_m is close to unity in most of the frequency interval and becomes strictly equal to one for interactions between two media across vacuum.

Armed with this result, we can now read-off the essence of the crossover behavior between retarded and non retarded interactions. In the limit of large L, the integrand is dominated by the exponential decay, $R̃(x*,ν)$ remains essentially constant, and the integral can be approximated to $R̃(x*,νT)∫(1+νL)e−νLdν$⁠. In the opposite limit, L → 0, the exponential function decays very slowly and the integrand becomes dominated by the algebraic decay of $R̃(x*,ν)$ that takes place at frequencies larger than a characteristic frequency of the material, ν_e. Whence, the integral is now given by $(1+νTL)e−νTL∫R̃(x*,ν)dν$⁠. This analysis highlights the physical origin of the crossover behavior and illustrates the mathematical complexity of the problem. Integrals with a crossover from algebraic to exponential decay are non-elementary functions that cannot be possibly expressed as a finite number of ordinary algebraic, exponential, and logarithmic functions.

In order to circumvent this pessimistic mathematical statement, we introduce an auxiliary exponential function, $e−ν/ν∞$⁠, made to mimic the algebraic decay of $R̃(x*,νT)$⁠. In this mapping, ν_∞ is an effective material parameter that dictates the range where such decay becomes effective. With the help of this function, we then write the trivial identity

The function inside the square brackets shares broadly the properties of the exact integrand in Eq. (6) and remains convergent for all L. Whence, we can use it as a reference weight function and approximate the full integral using again a one-point generalized Gauss-Laguerre quadrature. This leads readily to the expression (see the supplementary material)

where $R̃ξ*(L)=R̃(x*,ν*)eξ$⁠, ν^* = ν_T + ν_∞ξ, and ξ is an adimensional factor given by

The results of Eqs. (8) and (9), together with Eq. (5), provide an analytical approximation that describes the qualitative behavior of A_ω>0(L) in the full range from L = 0 to L → ∞. We call this the Weighted Quadrature Approximation (WQA). Our analysis allows us to identify two different inverse length scales, ν_∞ and ν_T, which dictate the L dependence of the Hamaker function. In the ensuing discussion, we describe the qualitative behavior that follows from the WQA, merely by assuming that the wave numbers ν_∞ and ν_T are sufficiently separated.

• For $L≪ν∞−1$⁠, we retain only the leading order term of Eq. (8) in the limit L → 0, and we find

A_ω>0 is a constant independent of L, with

Imposing the unknown parameter $ω∞=ν∞c/2ϵm1/2$ such that Eq. (10) recovers the known value of A_ω>0 at L = 0, our result becomes exact in this limit by construction.

• For $ν∞−1≪L≪νT−1,Aω>0$ recovers the expected inverse power dependence on L^3,5

$R*$ is now a complicated function of the relative permitivities, while ω_∞ ≫ ω^* ≈ c/L ≫ ω_T.

• For $L≫νT−1$⁠, the retarded interactions become strongly suppressed due to the exp(−ν_TL) factor

with

The expression in Eq. (13) is the exact result for large L at finite temperature.²⁰ In this limit, all the dielectric functions are calculated at iω_T and amount barely to $ϵi(iωT)=ni2$⁠, where n_i may be identified in simple materials with the refractive index in the visible region.^12,14 Notice that A_ω>0 vanishes for $L>νT−1$ and only the static contribution to the van der Waals forces [Eq. (3)] remains.^14,16 Often, in analytical calculations, the low temperature limit ν_T → 0 has been considered so that this exponential suppression does not take place.^3,14,21 At ambient temperature, however, this can become a serious error for L in the micrometer range.

In the above paragraph, it has been the aim to emphasize the crossover of A_ω>0 as L increases. Therefore, only the leading order athermal contributions have been retained. However, by considering also next to leading terms, it is possible to trace non-retarded contributions that operate within the retardation dominated regime. Particularly, for distances Lν_∞ ≫ 1, A_ω>0 features a non-retarded term of order k_BT that adds up to the A_ω=0 contribution of the full Hamaker function. Furthermore, from the analysis it follows that the thermal contribution to the retarded interactions is a small fraction of order ν_T/ν_∞ for Lν_T ≪ 1, but grows steadily and becomes the only remaining source of retardation for Lν_T ≫ 1. This qualifies analytically the observation made from numerical results on the significance of thermal contribution to van der Waals interactions.^{14,20,22–24}

For specific applications, it is required to consider explicitly the frequency dependence of the material’s dielectric response. In view of the limited information that is usually available, it is customary to describe the full dielectric function by a single classical damped oscillator such that

where n_α is the refractive index of the material α in the visible and ω_e is a characteristic absorption frequency with order of magnitude similar to the material’s ionization potential.^{12,14,18,19,25} This expression for the dielectric functions is quite convenient because it provides well known analytical results for the Hamaker constant (L = 0) under the approximation of Eq. (6).^{12,14,18,19,25} The corresponding results may be used to gauge the unknown parameter ν_∞ required for quantitative validation of Eqs. (8) and (17).

We use this model to describe the paradigmatic system of two mica plates interacting across vacuum,^4,5,23,26 with the high frequency oscillator parameters as given by Bergstrom.²⁵ Matching Eq. (8) to the Tabor-Winterton approximation for the Hamaker constant readily provides $ν∞=1.1868(n12+1)1/2ωe/c$ for the cut-off frequency (see the supplementary material). This is all that is required to plot A(L) for arbitrary values of L. A comparison of the exact first order Hamaker function obtained by integration of Eq. (1) shows very good agreement with the WQ approximation (Fig. 1). Particularly, WQA predicts the Hamaker function at large distances almost exactly, but yields a decay rate that is somewhat too slow. In practice, under the approximation of Eq. (6), we notice that $R̃ξ(x*,ν)$ is a bounded function of L, with a quadratic low order expansion in L. Therefore, the auxiliary function exp(−ν/ν_∞) should also exhibit a similar L decay. This can be achieved assuming $ν∞(L)=(1+k(ν∞0L)2)/(1+(ν∞0L)2)ν∞0$⁠, with k determined such that Eq. (8) yields the exact leading order correction to the Hamaker function as given by Eq. (4). Whereas forcing this requirement yields an exceedingly complicated algebraic expression, we illustrate this point by assuming k = 37/20 empirically and call this the WQ-k approximation (Fig. 1). With this device, we obtain almost exact agreement with the Hamaker function for all L. Notice the exact result from Lifshitz theory, as well as the models, is slightly not in line with experimental results that are accompanied in the figure for comparison. However, the early measurements with a prototype surface force apparatus are presently considered as order of magnitude estimates and likely suffer from a number of technical difficulties.^23,26

Unfortunately, even with ν_∞ assumed constant, the WQ approximation may be somewhat too cumbersome for some practical applications. Fortunately, we can obtain a simpler expression for A_ω>0(L) by taking into account that $R̃(x*,ν)$ is expected to be a monotonically decaying function in most cases (an interesting exception is the system made of ice/water/air). Accordingly, applying the second mean value theorem of definite integrals to Eq. (6), one can write

The above result is an exact quadrature for Eq. (6), provided one chooses a suitable L dependent high frequency cutoff for the wave-number, ν_∞. This frequency acts effectively as a natural ultra-violet cutoff for the integral. Performing now a trivial integration yields

where the factor 3/2 inside the first round parenthesis follows by inclusion of Euler-McLaurin corrections to lowest order. This accounts for the transformation of the sum [Eq. (4)] into an integral.

For practical matters, assuming a constant value for ν_∞ is sufficient. Indeed, by matching ν_∞ such that A_ω>0(L) yields the exact quadrature of Eq. (6) for L = 0, we find

where $n¯im=(ni2+nm2)1/2$ are root mean square indexes of refraction and $t=π/2$ is a numerical factor (see the supplementary material). We call this the Frequency Cutoff approximation (FCA).

The model is less accurate than the WQ approximation, but yields a qualitatively good agreement in the full range of relevant distances by using exactly the same parameters that are required to describe the Hamaker constant in usual applications.¹² This contrasts with the few empirical approximations that have been previously suggested, which only provide heuristic estimates of the first crossover length-scale and completely neglect the second one.^6,21,27 For ultimate simplicity, we can use Eq. (17) with the smooth function $R̃ξ=0$ replaced by a constant $R̃ξ=0(L=0)/t$ and choose t = 4 to improve the decay rate of A_ω>0 (see the supplementary material). In this way, the full Hamaker function may be readily given as

Comparison of this very simple procedure (FC-t approximation) yields again an overall very good description of the full Hamaker function, at the cost of somewhat deteriorating the large L behavior. We stress however, that the factor t = 4 is model independent so that this approach can now be applied generally for any system with overall very good accuracy. This is illustrated for interactions between two mica plates across water and for that of octane adsorbed on water (Fig. 2), with dielectric relaxation parameters taken from Israelachvili.¹² Particularly, our approach is able to capture the sign reversal of the Hamaker function for the water/octane/air system, where a description based on the Hamaker constant alone would be unable to predict the stabilization of thick octane films at the air/water interface.

As a caveat, notice that in either the WQ and FC approximations, the tractability of this approach relies on a one to one mapping of the effective cutoff frequency, ν_∞ with the resonance frequency ω_e of Eq. (15). Whence, the method is somewhat limited to this simple optical description. For a more general sum of states, either Eq. (8) or Eq. (17) will still provide a correct qualitative description by mapping ν_∞ to the highest energy oscillator, provided the optical response function R(x, n) decays monotonously.

Aside from the quantitative description discussed above for selected systems, our approach illustrates qualitatively a number of very relevant issues that despite efforts (cf. Parsegian,¹⁴ French et al.¹⁷) are generally not recognized. First, the crossover from retarded to non-retarded interactions sets in for distances of the order c/ω_e. For a large number of materials, including inorganic substrates and hydrocarbons, this corresponds to distances of about 10 nm,^{12,14,18,23,28,29} an order of magnitude less than it is often assumed.²⁷ Second, the rather simple crossover function $A≈A(0)(1+ν∞L)−1$ that is most often used in the literature yields a Hamaker constant that decays as L⁻¹ for large distances.²⁷ In practice, the length scales ν_∞ and ν_T are not sufficiently well separated at ambient temperature, and the decay of the van der Waals interactions does never really attain the L⁻³ power law expected from the Gregory equation. This can be illustrated by representing the effective exponent n = d ln g/d ln L.^5,23 Clearly, a regime of n = 3 constant never really sets in at ambient temperature (Fig. 1). Rather, n reaches a maximum value that is close, but smaller than the ideal value of 3 expected from the standard retardation regime. This observation is not specific to a particular choice of systems (Fig. 2). It is dictated by the separation between the two constants ν_∞ and ν_T, which is just two orders of magnitude at ambient temperature for most substances. In fact, using the FC approximation, it is possible to show that the Hamaker function decreases to half the Hamaker constant at distances $L1/2=4ν∞−1$⁠, whence, only a few times larger than the onset of retardation. Using this estimate for L_1/2 in either Eq. (8) or in Eq. (17) shows the Hamaker function obeys a scaling form A_ω>0(L) = A_ω>0(0)u(L/L_1/2), with u(x) a universal function, as observed empirically by Cheng and Cole.²¹ As a final comment, we remark that the Gaussian Quadrature, which is a well known tool for numerical integration, can be exploited as an exceptional theoretical method in physics provided one chooses suitable problem adapted weight functions.

In summary, it was shown that a simple analytical crossover function [Eq. (19)] can be used to describe the Hamaker function in the range from the angstrom to the micrometer using exactly the same few material parameters that are employed conventionally to calculate the Hamaker constant. It is hoped this will allow for a much better quantification of intermolecular forces in a wide range of applications at the nanometric scale.^17,30–33

See supplementary material for a detailed explanation of the problem adapted one point Gauss quadrature rules, detailed derivation of Eqs. (1), (4), (8), and (17), calculation of the high-frequency cutoff, ν_∞, and parameters for the model systems.

I would like to thank Andrew Archer and David Sibley for their hospitality and stimulating ambience at the Department of Mathematical Sciences at the University of Loughborough where part of this work was performed. I also acknowledge funding from the “Programa Estatal de Promoción del Talento y su Empleabilidad en I+D+i,” by the Spanish Ministerio de Educación, Cultura y Deporte (Plan Estatal de Investigación Científica y Técnica y de Innovación 2013-2016) for the generous funding of this visit under Grant Ref. No. PRX17/00460 and to Agencia Estatal de Investigación and Fondo Europeo de Desarrollo Regional (FEDER) for funding under Grant No. FIS2017-89361-C3-2-P.