On the spectrum of fluctuations of a liquid surface: from the molecular scale to the macroscopic scale

On the spectrum of fluctuations of a liquid surface: from the molecular scale to the macroscopic scale

Blokhuis, E.M.

Citation

Blokhuis, E. M. (2009). On the spectrum of fluctuations of a liquid surface: from the molecular scale to the macroscopic scale. Journal Of Chemical Physics, 130(1), 014706.

doi:10.1063/1.3054346

Version: Not Applicable (or Unknown)

License: Leiden University Non-exclusive license Downloaded from: https://hdl.handle.net/1887/62397

Note: To cite this publication please use the final published version (if applicable).

On the spectrum of fluctuations of a liquid surface: From the molecular scale to the macroscopic scale

Edgar M. Blokhuis

Citation: The Journal of Chemical Physics 130, 014706 (2009); doi: 10.1063/1.3054346 View online: https://doi.org/10.1063/1.3054346

View Table of Contents: http://aip.scitation.org/toc/jcp/130/1 Published by the American Institute of Physics

Articles you may be interested in

Triezenberg-Zwanzig expression for the surface tension of a liquid drop The Journal of Chemical Physics 138, 194711 (2013); 10.1063/1.4805001

The intrinsic structure of the water surface

The Journal of Chemical Physics 125, 014709 (2006); 10.1063/1.2209681

Capillary waves at the liquid-vapor interface and the surface tension of water The Journal of Chemical Physics 125, 014702 (2006); 10.1063/1.2209240

Thermodynamic expressions for the Tolman length

The Journal of Chemical Physics 124, 074701 (2006); 10.1063/1.2167642

Test-area simulation method for the direct determination of the interfacial tension of systems with continuous or discontinuous potentials

The Journal of Chemical Physics 123, 134703 (2005); 10.1063/1.2038827

The Statistical Mechanical Theory of Surface Tension

The Journal of Chemical Physics 17, 338 (1949); 10.1063/1.1747248

On the spectrum of fluctuations of a liquid surface: From the molecular scale to the macroscopic scale

Edgar M. Blokhuis^a兲

Colloid and Interface Science, Leiden Institute of Chemistry, Gorlaeus Laboratories, P.O. Box 9502, 2300 RA Leiden, The Netherlands

共Received 1 August 2008; accepted 1 December 2008; published online 7 January 2009兲

We show that to account for the full spectrum of surface fluctuations from low scattering vector qdⰆ1 共classical capillary wave theory兲 to high qdⲏ1 共bulklike fluctuations兲, one must take account of the interface’s bending rigidity at intermediate scattering vector qdⱗ1, where d is the molecular diameter. A molecular model is presented to describe the bending correction to the capillary wave model for short-ranged and long-ranged interactions between molecules. We find that the bending rigidity is negative when the Gibbs equimolar surface is used to define the location of the fluctuating interface and that on approach to the critical point it vanishes proportionally to the interfacial tension. Both features are in agreement with Monte Carlo simulations of a phase-separated colloid-polymer system. © 2009 American Institute of Physics.关DOI:10.1063/1.3054346兴

I. INTRODUCTION

The description of the spectrum of surface fluctuations of a liquid from the macroscopic scale down to the molecular scale remains a challenging experimental and theoretical problem. Using grazing incidence light scattering experi- ments, Daillant and co-workers¹were able, for the first time, to determine the full spectrum of surface fluctuations, where in previous experiments共ellipsometry and reflectivity兲 only certain aspects of the spectrum could be determined. At the same time, the spectrum can now be analyzed in computer simulations with ever increasing accuracy.^2–4

Theoretical insight into the structure of a simple liquid surface is provided by density functional theories共DFTs兲 on the one hand^5–7and the capillary wave 共CW兲 model on the other hand.^8–10 DFTs provide a description of the interface on a microscopic level. The prototype of such theories, the van der Waals squared-gradient model, was very successful in describing, for the first term, the density profile and sur- face tension in terms of molecular parameters.⁷ It, however, fails to capture the subtle role of long wavelength interfacial fluctuations described by the CW model.^8–10

The CW model introduced in 1965共Ref.8兲 describes the spectrum of fluctuations in terms of a height function h共rជ储兲 with the surface tension␴and gravity g acting as the domi- nant restoring forces. The length scale involved in describing CWs is the capillary length, Lc⬅

冑

␴/共m⌬␳g兲, which may be as large as a tenth of a millimeter. The theoretical challenge is to incorporate both theories and to describe the spectrum of fluctuations of a liquid surface, as determined from light scattering experiments and computer simulations, from the molecular scale to the scale of CWs.

An important ingredient in “bridging the gap” between CWs and the molecular scale is an extension of the CW model that incorporates the energy associated with bending the interface.^11–13Bending is important when the wavelength

of the height fluctuations is approximately

冑

k_BT/␴^{, which is} typically of the order of a few times the molecular diameter, i.e., close to the scale where the molecular structure becomes important and the density fluctuations are more bulklike. The natural question that arises is whether it is possible to de- scribe the full spectrum of surface fluctuations by the CW model at long wavelengths and bulklike fluctuations at the molecular scale. Is it then necessary to include the leading order correction to the CW model from bending or are even higher order terms, relevant at even smaller length scales, required?

This article addresses these questions in two parts 共a condensed version has appeared in Ref.14兲. In the first part, we analyze the spectrum of fluctuations recently obtained by Vink et al.³ in computer simulations of a phase-separated polymer-colloid system^15–18 in which the interactions are strictly short ranged 共SR兲. It is shown that the simulation data are very accurately described by the combination of the CW model extended to include a bending correction, with the bending rigidity as an adjustable parameter, and bulklike fluctuations.

In the second part, a molecular basis for the bending correction to the CW model is offered and the results are compared with the simulations. The theoretical framework used for the comparison is mean-field DFT in which the interactions are described by a nonlocal integral term.^5–7,19 The advantage of this approach is that it features the full shape of the interaction potential enabling the analysis of different forms and ranges of the interaction potential. We consider both SR interactions and long-ranged共LR兲 interac- tions, which fall of as U共r兲⬀1/r⁶ at large intermolecular separations.

An important ingredient in our theoretical analysis is the modification of the density profile, described by␳1共z兲, due to the local bending of the interface.²⁰ The determination of

␳1共z兲 requires one to formulate precisely the thermodynamic conditions used to vary the interfacial curvature. Several ap-

a兲Electronic mail: e.blokhuis@chem.leidenuniv.nl.

共2009兲

0021-9606/2009/130共1兲/014706/15/$25.00 130, 014706-1 © 2009 American Institute of Physics

proaches for the determination of ␳1共z兲 have appeared in literature.^19–22 They differ in the form of the external field used to set the curvature to a specific value; in the equilib- rium approach²⁰ the external field is uniform throughout the system, whereas in the approach by Parry and Boulter²¹ and Blokhuis et al.²² it is infinitely sharp-peaked 共Vext⬀␦共z兲兲 at the interface. In this article we suggest to add an external field acting in the interfacial region only with a peak width of the order of the thickness of the interfacial region. The ad- vantage of this approach is that the bulk regions are unaf- fected by the addition of the external field and the resulting

␳1共z兲 is a continuous function.

Our paper is organized as follows: In Sec. II, the general form of the surface structure factor to describe the spectrum of interfacial fluctuations is derived as the combination of the CW model extended to include a bending correction and bulklike fluctuations. This form is then compared in Sec. III to the Monte Carlo共MC兲 simulation results by Vink et al.³ for the phase-separated polymer-colloid system. In Sec. IV, the mean-field DFT used to provide a molecular basis for the bending extension to the CW model is presented. Explicit results are obtained for SR interactions 共Sec. V兲 and LR interactions共Sec. VI兲. We end with a discussion of results.

II. THE FLUCTUATING LIQUID SURFACE

In the classical CW model, the fluctuating interface is described by a two-dimensional surface height function h共rជ储兲, where rជ储=共x,y兲 is the direction parallel to the surface.^8–10 The fluctuating density profile can then be written in terms of an “intrinsic density profile” shifted over a distance h共rជ储兲 as follows:

␳共rជ兲 =␳0共z − h共rជ储兲兲, 共1兲

where␳0共z兲 is the intrinsic density profile. Often, fluctuations are assumed to be small so that an expansion in h can be made, neglecting terms ofO共h²兲,

␳共rជ兲 =␳0共z兲 −␳0⬘共z兲h共rជ储兲 + ¯ . 共2兲 An important consequence of the above linearization is that one may now identify the intrinsic density profile as the av- erage density profile,␳0共z兲=具␳共rជ兲典, in view of the fact that 具h共rជ储兲典=0. It is convenient to locate the z=0 plane such that it coincides with the Gibbs equimolar surface,^7,23i.e.,

冕

−⬁

⬁

dz关具␳共rជ兲典 −␳step共z兲兴 =

冕

−⬁

⬁

dz关␳0共z兲 −␳step共z兲兴 = 0, 共3兲 where ␳step共z兲=␳_ᐉ⌰共−z兲+␳v⌰共z兲 with ⌰共z兲 the Heaviside function and ␳_ᐉ,v the bulk density in the liquid and vapor regions, respectively.

In the above model for␳共rជ兲, the density correlations are essentially given by the correlations of h共rជ储兲, which are de- scribed by the height-height correlation function

Shh共r储兲 ⬅ 具h共rជ1,储兲h共rជ2,储兲典, 共4兲 where rជ储⬅rជ2,储− rជ1,储 and r储⬅兩rជ储兩.

To determine the height-height correlation function, one should examine the change in free energy, ⌬⍀, associated with a fluctuation of the interface. In the CW model it is described by considering the change in free energy associ- ated with a distortion of the surface against gravity and sur- face area extension,⁸

⌬⍀ =1

2

冕

^{drជ储关m⌬␳gh共rជ储兲2+␴兩ⵜជh共rជ储兲兩2兴. 共5兲}

It is convenient to express⌬⍀ in terms of the Fourier Trans- form of h共rជ储兲, h共qជ兲=兰drជ储e^{−iqជ·rជ储}h共rជ储兲,

⌬⍀ =1

2

冕

共2^dq␲^ជ兲²关m⌬␳^{g +}␴^q2兴h共qជ兲h共− qជ兲. 共6兲 In the CW model the height-height correlation function is determined by a full statistical mechanical analysis^9,10 in which the above expression for the change in free energy is interpreted as the so-called CW Hamiltonian, ⌬⍀

=HCW关h共rជ储兲兴. In general, one has Shh共r储兲 =1

Z

冕

^{Dhh共rជ1,储兲h共rជ2,储兲e−HCW关h兴/kBT, 共7兲}

where Z is the partition function associated withHCW关h兴, kB

is Boltzmann’s constant, and T is the temperature. It can be shown that^9,10

Shh共q兲 =

冕

^{drជ储e−iqជ·rជ储Shh共r储兲 =}m⌬␳^{kg +BT}␴^q2= k_BT

␴共Lc−2+ q²兲. 共8兲 For simplicity, we ignore gravity effects in the following and set Lc=⬁ 共g=0兲.

A. Extended capillary wave model

In the derivation of the classical CW model, one as- sumes an expansion in gradients of h共rជ储兲, 兩ⵜជh兩Ⰶ1. In the extended capillary wave共ECW兲 model, one wishes to extend the expansion by including higher derivatives of h共rជ储兲. To leading order one may then write the fluctuating density as^19,22

␳共rជ兲 =␳0共z兲 −␳0⬘共z兲h共rជ储兲 −␳1共z兲

2 ⌬h共rជ储兲 + ¯ . 共9兲 The function ␳1共z兲 is identified as the correction to the den- sity profile due to the curvature of the interface, ⌬h共rជ储兲⬇

−1/R1− 1/R2, with R₁ and R₂ the共principal兲 radii of curva- ture. The prefactor of −1/2 is chosen such that the notation is consistent with an analysis in which the curvature does not result from a fluctuation of the planar interface, but is due to the fact that one considers a spherical liquid droplet 共R1

= R₂= R兲 in 共metastable兲 equilibrium with a bulk vapor phase.^20,24,25 An expansion in the curvature of the density profile␳s共r兲 then gives

␳s共r兲 =␳0共r兲 +␳1共r兲

R + ¯ , 共10兲

which parallels the expansion in Eq.共9兲.

The inclusion of curvature corrections in the ECW model leads to higher order terms in an expansion in q², terms beyond␴^q2, in the expression for⌬⍀ in Eq.共6兲. It is customary to capture these higher order terms by introducing a wave vector dependent surface tension␴共q兲 as follows:¹²

⌬⍀ =1

2

冕

共2^dq␲^ជ兲²␴共q兲q²h共qជ兲h共− qជ兲, 共11兲 which gives for the height-height correlation function

S_hh共q兲 = k_BT

␴共q兲q². 共12兲

The precise form of␴共q兲 depends sensitively on the behavior of the interaction potential at large distances.¹⁹ When the interaction potential is sufficiently SR, the expansion of␴共q兲 in q²is regular and the leading correction is of the form

␴共q兲 =␴^{+ kq2+}O共q⁴兲 共SR兲. 共13兲 The coefficient k is identified as the bending rigidity.^11–13 This is because the form for⌬⍀ in Eq.共11兲, with␴共q兲 given by Eq.共13兲, can also be derived from the Helfrich free en- ergy expression,¹¹which reads for a fluctuating interface

⌬⍀ =1

2

冕

^{drជ储关␴兩ⵜជh共rជ储兲兩2+ k共⌬h共rជ储兲兲2兴. 共14兲}

When the interaction potential is LR, specifically when it falls of as U共r兲⬀1/r⁶ at large intermolecular distances, which is the case for regular fluids due to London-dispersion forces, one finds that the leading correction to␴共q兲 picks up a logarithmic contribution¹⁹

␴共q兲 =␴^{+ k}sq²ln共qᐉk兲 + O共q⁴兲 共LR兲, 共15兲 with ksandᐉkparameters independent of q. The coefficient ks depends on the asymptotic behavior of U共r兲 but is other- wise a universal constant.¹⁹The bending lengthᐉkdepends, like the bending rigidity k, on the microscopic parameters of the model. In principal, all the parameters␴^{, k, k}s, andᐉkcan be expressed in terms of the density profiles␳0共z兲 and␳1共z兲 by inserting the fluctuating density as given in Eq.共9兲into a microscopic model for the free energy and comparing the result with Eq.共11兲.

It is important to realize that the ECW model assumes a curvature expansion in Eq. 共9兲, which translates into an ex- pansion in q² in Eq. 共11兲 that is valid only up to O共q⁴兲.

Higher order terms are not systematically included. The re- sult is that one should limit the expansion of␴共q兲 in Eq.共13兲 or Eq.共15兲to the order in q indicated.

B. Definition of the height profile

An important subtlety in the preceding analysis is the fact that the location of the interface, i.e., the value of the height function h共rជ储兲, cannot be defined unambiguously.²³A certain procedure must always be formulated to determine h共rជ储兲. It turns out that the choice for h共rជ储兲 influences the density profile␳1共z兲, which, in turn, determines the value of the bending parameters k andᐉk.

We explicitly consider two canonical choices for the de- termination of h共rជ储兲; the crossing constraint 共cc兲 and the in-

tegral constraint共ic兲.^21,22Other choices are certainly possible and equally legitimate as long as they lead to a location of the dividing surface that is “sensibly coincident” with the interfacial region.²³ In this context we like to mention the work by Tarazona et al.,⁴ who define the location of the interface based on the distribution of molecules rather than the molecular density alone.

In the cc, h共rជ储兲 is defined as the height where the fluctu- ating density equals some fixed value of the density that lies in between the limiting bulk densities, say,␳共rជ兲=␳0共z=0兲,

␳共rជ储,z = h共rជ储兲兲 =␳0共0兲 共cc兲. 共16兲 Using this condition in Eq.共9兲, one finds the following con- straint for␳1共z兲=␳1cc共z兲:

␳1cc共0兲 = 0. 共17兲

In the ic, h共rជ储兲 is defined by the integral over the fluctuating density²³

h共rជ储兲 = 1

⌬␳

冕

−⬁

⬁

dz关␳共rជ兲 −␳step共z兲兴 共ic兲. 共18兲

With this condition inserted into Eq.共9兲, one now finds that

␳1共z兲=␳1ic共z兲 is subject to the following constraint:

冕

−⬁

⬁

dz␳1

ic共z兲 = 0. 共19兲

We show in Sec. IV that the ambiguity in locating the divid- ing surface translates into the density profile␳1共z兲 being de- termined up to an additive factor proportional to␳₀⬘共z兲.²² In particular, ␳₁^ic共z兲 and␳₁^cc共z兲 are related to

␳1 ic共z兲 =␳1

cc共z兲 +␣␳₀⬘共z兲. 共20兲

The value of the constant␣can be determined by integrating both sides of the above equation over z,

␣⁼_⌬¹

␳

冕

_−⬁^{⬁ dz␳1cc共z兲. 共21兲}

One may further show that the ambiguity in the determina- tion of ␳1共z兲 is of influence to the value of the bending pa- rameters k andᐉk. In Sec. IV we show that because␳1

ic共z兲 and

␳1

cc共z兲 are related to Eq. 共20兲, we have for the bending parameters²²

k^ic= k^cc−␣␴^,

共22兲 ksln共ᐉk

ic兲 = ksln共ᐉk cc兲 −␣␴^.

Naturally, all experimentally measurable quantities cannot depend on the choice made for the location of the height function h共rជ储兲. The implication is that it is necessary to for- mulate precisely the quantity that is determined experimen- tally and verify that its value is independent of the choice for h共rជ储兲. This is explicitly shown next.

The quantity studied in experiments and simulations is the共surface兲 density-density correlation function. It is an in- tegral into the bulk region to a certain depth L of the density- density correlation function

S共r储兲 ⬅ 1 共⌬␳兲²

冕

−L

L

dz₁

冕

−L L

dz₂具关␳共rជ1兲 −␳step共z1兲兴

⫻关␳共rជ2兲 −␳step共z2兲兴典. 共23兲 When we insert the general expression for␳共rជ兲 as given by Eq.共9兲 into Eq.共23兲, one finds that

S共r储兲 = 具h共rជ1,储兲h共rជ2,储兲典 − 1

⌬␳

冕

−⬁

⬁

dz␳1共z兲具h共rជ1,储兲⌬h共rជ2,储兲典, 共24兲 where we can neglect a term 具⌬h⌬h典 to the order in the curvature expansion considered. Furthermore, we have as- sumed that L is sufficiently large so that we can approximate

冕

−L L

dz␳0⬘共z兲 ⬇ − ⌬␳^,

冕

−L 共25兲

L

dz␳1共z兲 ⬇

冕

−⬁

⬁

dz␳1共z兲.

Rather than S共r储兲, we consider its Fourier transform, S共q兲, which we shall term the surface structure factor,

S共q兲 =

冕

^{drជ储e−iqជ·rជ储S共r储兲}

= Shh共q兲 + 1

⌬␳

冕

_−⬁^{⬁ dz␳1共z兲q2Shh共q兲. 共26兲}

We now verify that S共q兲 is independent of the choice for h共rជ储兲 by determining S共q兲 using both the ic and cc. For sim- plicity, we consider the case of SR forces only共the verifica- tion for the case of LR forces follows analogously兲. The surface structure factor using both constraints is given by

S^cc共q兲 = k_BT

␴q²+ k^ccq⁴+¯+ ␣k_BTq²

␴q²+ k^ccq⁴+¯

=k_BT

␴^q2− k_BTk^cc

␴^{2 +}

␣^kBT

␴ ^+O共q2兲,

共27兲 S^ic共q兲 = k_BT

␴^{q2+ kicq4+}¯=k_BT

␴^q2− k_BTk^ic

␴^{2 +O共q2兲,} where we have used the explicit expression for S_hh共q兲 in Eq.

共12兲 together with Eq. 共13兲. On account of the fact that k^ic

= k^cc−␣␴, one finds that S^cc共q兲=S^ic共q兲⬅S共q兲 as required.

This analysis shows that S共q兲 equals the height-height correlation function when the integral constraint is used to define the location of the height profile, i.e.,

S共q兲 = Shhic共q兲. 共28兲

It is therefore convenient, but by no means necessary, to use the ic to define the location of the dividing surface.

Finally, we consider the contribution of “bulklike” fluc- tuations to the fluctuating density profile, which are predomi- nantly present at short wavelengths, qdⲏ1.

C. Bulklike fluctuations

Adding short wavelength, bulklike fluctuations to the fluctuating density, the full picture that emerges for ␳共rជ兲 is that schematically depicted in Fig.1. It can be described as

␳共rជ兲 =␳0共z兲 −␳0⬘共z兲h共rជ储兲 −␳1共z兲

2 ⌬h共rជ储兲 +␦␳b共rជ兲, 共29兲 where ␦␳b共rជ兲 represents the bulklike fluctuations. We shall consider only small fluctuations so that具␦␳b典=0 and assume that there are no correlations between height fluctuations and bulklike fluctuations,具h␦␳b典=0. When we insert the expres- sion for ␳共rជ兲 as given by Eq. 共29兲 into the expression for S共r储兲 in Eq.共23兲, one finds that

S共r储兲 = Shh

ic共r储兲 + 1 共⌬␳兲²

冕

−L

L

dz₁

冕

_−⬁^{⬁ dz12具␦␳b共rជ1兲␦␳b共rជ2兲典.}

共30兲 Here we have made a further approximation by replacing the integration over z₂ from −L to L with an integral over z₁₂ from −⬁ to ⬁. The integral over z1that is left gives rise to a term that increases linearly with L. That means that the bulk- like contributions to S共r储兲 eventually dominate the height fluctuations when L becomes larger. To study surface fluc- tuations via S共r储兲 it is therefore important that on the one hand L is sufficiently large in order to make the approxima- tions in, e.g., Eq.共25兲but on the other hand not so large as to completely dominate the contribution from surface height fluctuations. In Sec. III we show how these two conditions pan out for the circumstances under which the simulation results are obtained.

A further issue is that the bulk density correlation func- tion具␦␳b␦␳b典 differs in either phase 共liquid or vapor兲. When one then considers the integral over z₁, it seems appropriate to approximate 具␦␳b␦␳b典 by the density correlation function in the bulk liquid region

具␦␳b共rជ1兲␦␳b共rជ2兲典 =␳_ᐉ²关g_ᐉ共r兲 − 1兴 +␳_ᐉ␦共rជ12兲, 共31兲 and introduce an L-dependent prefactorNLto account for the integral over z₁. The surface structure factor thus becomes

S共q兲 = Shh

ic共q兲 + NLSb共q兲, 共32兲

with the bulk structure factor Sb共q兲 defined as

0 z

ρ(z)

>

>

>

h

FIG. 1. Sketch of the fluctuating density profile as a function of z; the height h = h共rជ储兲 is the distance over which the intrinsic density profile␳0共z兲 共dashed line兲 is shifted.

S_b共q兲 = 1 +␳_ᐉ

冕

^{drជ12e−iqជ·rជ12关gᐉ共r兲 − 1兴. 共33兲}

This approximation may be justified by arguing that close to the critical point there is no distinction between the two bulk correlation functions, whereas far from the critical point the contribution from the bulk vapor can be neglected since ␳v

⬇0.

The value for the L-dependent prefactorNLmay be de- termined from a fit to the limiting behavior of S共q兲 at q→⬁. For an explicit evaluation of Sb共q兲, we have taken for g_ᐉ共r兲 the Percus–Yevick 共PY兲 solution²⁶for the hard-sphere correlation function, g_ᐉ共r兲=g_hs^PY共r;␳ᐉ兲.

III. COMPARISON WITH MONTE CARLO SIMULATIONS In this section, the surface structure factor in Eq.共32兲is compared to the results from MC simulations by Vink et al.³ The system considered consists of a mixture of colloidal par- ticles with diameter d and polymer particles with diameter 2R_g. The colloid-colloid and colloid-polymer interactions are considered to be hard-sphere-like, whereas polymer-polymer interactions are taken ideal. The presence of polymer induces a depletion attraction between the colloidal particles, which may ultimately lead to phase separation.^15–18 The resulting interface of the demixed colloid-polymer system is studied by Vink et al.³for a number of polymer concentrations and for a polymer-colloid size ratio parameter ␧⬅1+2Rg/d

= 1.8.

To study the interfacial fluctuations, Vink et al.³ intro- duced the local interface position as

zG共rជ储兲 ⬅ 1

⌬␳

冕

−L L

dz关␳共rជ兲 −␳step共z兲兴, 共34兲

where␳共rជ兲 can be taken to be either the colloid or polymer density. The integration limits ⫾L are inside the bulk re- gions, but different values for it are systematically considered.³ One may easily verify that the correlations of the local interface position are exactly described by the sur- face structure factor defined earlier in Eq.共23兲,

具zG共rជ1,储兲zG共rជ2,储兲典 = S共r储兲. 共35兲 In Fig.2, typical results for the Fourier transform of the

surface structure obtained in the MC simulations of Vink et al.³are shown共Fig. 13 of Ref.3兲. In this example the inte- gration limit is varied, L/W=1, 2, 3, 4, where W is some measure of the interfacial thickness. One clearly observes that when L/W is too small, the results do not match the classical CW behavior for small q共dashed line兲, and that the contribution from bulklike fluctuations at high q increases with L/W.

In Fig. 3, we consider the result from Fig. 2 for L/W

= 3. For small q the results asymptotically approach the result of the classical CW model 共dotted line兲 with the value of␴ taken from separate simulations. The dashed line is the com- bination of the CW model with the bulk correlation function

S共q兲 =k_BT

␴^{q2+NLSb共q兲. 共36兲}

The value of NL is chosen such that it matches the q→⬁

limit for S共q兲 in Fig. 3. One finds that Eq. 共36兲 already matches the simulation results quite accurately except at in- termediate values of q, qd⬇1.

As a next step, we investigate whether the inclusion of a bending rigidity is able to describe the simulation results at the following intermediate values:

S共q兲 = k_BT

␴^{q2+ kicq4+}¯+NLS_b共q兲. 共37兲

The bending rigidity describes the leading order correction to the classical CW model in an expansion in q². Its value is therefore obtained from analyzing the behavior of S共q兲 when qdⱗ1. The fact that the simulation results in Fig.3are sys- tematically above the CW prediction in this region indicates that the bending rigidity thus obtained is negative, k^ic⬍0.

Unfortunately, a negative bending rigidity prohibits the use of Eq.共37兲to fit the simulation results in the entire q-range since the denominator becomes zero at a certain value of q.

It is therefore convenient to rewrite the expansion in q² in Eq. 共37兲in the following form:

1 10

q

1 10

S(q)

L/W = 4 L/W = 3 L/W = 2 L/W = 1 CW

FIG. 2. MC results by Vink et al.共Ref.3兲 for the surface structure factor 共in units of d⁴兲 vs q 共in units of 1/d兲 for various values of the integration limit L/W=1, 2, 3, 4. The dashed line is the CW model. In this example, ␧

= 1.8,␩p= 1.0, and the colloidal particles are used to define zG.

1 10

q

1 10

S(q)

Vink CW CW + bulk ECW + bulk

FIG. 3. MC results by Vink et al.共Ref.3兲 for the surface structure factor 共in units of d⁴兲 vs q 共in units of 1/d兲. The dotted line is the CW model, the dashed line is the combination of the CW model and the bulk correlation function, and the drawn line is the combination of the ECW model and the bulk correlation function. In this example,␧=1.8,␩p= 1.0, L/W=3, and the colloidal particles are used to define zG.

S共q兲 =k_BT

␴^q2

冉

^{1 −k}␴^icq2+¯

冊

^{+NLSb共q兲, 共38兲}

which is equivalent to Eq.共37兲to the order in q considered, but which has the advantage of being well behaved in the entire q-range. Other forms to regulate S共q兲, which are equivalent to Eq. 共37兲 to the order in q considered, may certainly be formulated. In analogy with a similar treatment of CWs by Parry et al.²⁷in the context of wetting transitions, one might suggest that the appearance of a negative bending rigidity indicates the missing of a correlation length that would replace Eq.共38兲with an explicit formula valid for all values of q, not just to the order in q considered.

The above form for S共q兲 in Eq.共38兲, with the bending rigidity used as an adjustable parameter 共k=−0.045kBT兲, is plotted in Fig.3as the drawn line. Exceptionally good agree- ment with the MC simulations is obtained. In TableI, we list values of the bending rigidity obtained for a number of poly- mer volume fractions,␩p. These values are the results of fits of S共q兲 from MC simulations for several system sizes and for several values of L/W, with the error estimated from the standard deviation of the various results. For the L/W=1 and L/W=2 curves 共see Fig.2兲, one needs to adjust for the fact that the CW limit is not correctly approached at low q. For the results in Fig. 2, one ultimately obtains for the bending rigidity k = −0.040, −0.040, −0.045, −0.060k_BT, for L/W=1, 2, 3, and 4, respectively.

In TableI, it should be reminded that, rather than the true polymer volume fraction in either phase,␩pshould be inter- preted as the polymer volume fraction of a reservoir fixing the polymer chemical potential.²⁸ Furthermore, the “liquid”

is defined as the phase relatively rich in colloids and the

“vapor” as the phase relatively poor in colloids.

The excellent agreement between Eq.共38兲 and the MC simulations is even more apparent in Fig.4where the results in Fig. 3 are redrawn on a linear scale. In Fig. 4 we also show the simulation results³ and the corresponding fit using the polymer particles to define the location of the interface.

As the polymer-polymer interactions are considered ideal, the bulk structure factor Sb共q兲=1 in this case.

It is important to note that, effectively, the inclusion of a bending rigidity in the CW model results in the presence of an additive factor in S共q兲 that is adjusted, see Eq. 共38兲. The determination of the value for k^icfrom the behavior of S共q兲 near q = 0, therefore, requires one to take into account the presence of the bulklike fluctuations since they also contrib- ute as an additive constant, NLSb共0兲, near q=0. This means that even though the MC simulation results of Vink et al.³are very accurately described by Eq. 共38兲, the resulting value

obtained for k^icsensitively depends on the theoretical expres- sion used for S_b共0兲. Here we have simply approximated the bulk correlation function by the PY hard-sphere expression in the liquid,²⁶ but one could imagine more sophisticated expressions leading to a somewhat different value for k^ic.

In Secs. IV–VI we investigate whether the values for the bending rigidity obtained from the simulations共TableI兲 can also be described in the context of a molecular theory.

IV. DENSITY FUNCTIONAL THEORY

Our task in this section is straightforward. Using the expression for ␳共rជ兲 given in Eq.共9兲, we determine⌬⍀ and the resulting␴共q兲. To achieve this, we need a model for the free energy and define a procedure to determine the density profiles,␳0共z兲 and␳1共z兲, that are present in the expression for

␳共rជ兲. We choose to perform these tasks in the context of DFT.

In the DFT for an inhomogeneous system that we consider,^5–7,19 the free energy is given by the free energy of the reference hard-sphere system augmented by an integral nonlocal term that considers the attractive part of the inter- action potential, U共r兲=Uhs共r兲+Uatt共r兲,

⍀关␳兴 =

冕

^{drជ1ghs共␳兲 +12}

冕

^drជ1

冕

^{drជ2Uatt共r兲␳共rជ1兲␳共rជ2兲.}

共39兲 For explicit calculations, g_hs共␳兲 is taken to be of the Carnahan–Starling form²⁹

TABLE I. Listed are the simulation results 共Ref.3兲 for the polymer volume fraction␩p, liquid and vapor colloidal volume fractions,␩ᐉand␩v, surface tension␴共in units of kBT/d²兲, bending rigidity k 共in units of kBT;

in parentheses is the estimated error in the last digit兲, and冑−k/␴共in units of d兲.

␩p ␩ᐉ ␩v ␴ k 冑−k/␴

0.9 0.2970 0.0141 0.1532 −0.045共15兲 0.40

1.0 0.3271 0.0062 0.2848 −0.07 共2兲 0.50

1.1 0.3485 0.0030 0.4194 −0.10 共3兲 0.49

1.2 0.3647 0.0018 0.5555 −0.14 共3兲 0.50

0 2 4 6 8 10

q

0 2 4 6 8 10

S(q)

ECW + bulk CW Vink (col) Vink (pol)

FIG. 4. MC results by Vink et al.共Ref.3兲 for the surface structure factor 共in units of d⁴兲 vs q 共in units of 1/d兲 using the colloidal particles 共circles兲 and polymer particles共triangles兲 to define zG. The dotted line is the CW model and the drawn lines are the combination of the ECW model and the bulk correlation function. In this example,␧=1.8,␩p= 1.0, and L/W=3.

g_hs共␳兲 = kBT␳^ln共␳兲 + kBT␳^共4␩^{− 3}␩²兲

共1 −␩兲² −␮␳^, 共40兲 where ␩⬅共␲/6兲␳^d3. In the uniform bulk region, the free energy equals

⍀共␳兲

V ⬅ g共␳兲 = ghs共␳兲 − a␳^2, 共41兲 with the van der Waals parameter a given by⁷

a⬅ −1

2

冕

^{drជ12Uatt共r兲. 共42兲}

The integration over rជ12is restricted to the region r⬎d. This is not explicitly indicated; instead, we adhere to the conven- tion that the attractive part of the interaction potential U_att共r兲=0 when r⬍d. The chemical potential␮ is fixed by the condition of two-phase coexistence,␮⁼␮coex, which im- plies that␮coex,␳v, and␳_ᐉare determined from the following set of equations: g⬘^共␳v兲=0, g⬘^共␳_ᐉ兲=0, and g共␳v兲=g共␳_ᐉ兲=−p.

To determine the change in free energy due to density fluctuations, we insert the expression for ␳共rជ兲 given by Eq.

共9兲into the expression for ⍀ in Eq.共39兲. One then finds for

⌬⍀=⍀−␴^A

⌬⍀ =1

8

冕

^{drជ1兵ghs⬙共␳0兲␳1共z1兲2关⌬h共rជ1,储兲兴2}

^其

+ 1

8

冕

^drជ1

冕

^{drជ12Uatt共r兲兵− 2␳⬘0共z1兲␳0⬘共z2兲关h共rជ2,储兲}

− h共rជ1,储兲兴²+ 4␳1共z1兲␳⬘0共z2兲⌬h共rជ1,储兲关h共rជ2,储兲 − h共rជ1,储兲兴 +␳1共z1兲␳1共z2兲⌬h共rជ1,储兲⌬h共rជ2,储兲其. 共43兲 Even though the derivation is somewhat different, this ex- pression equals that given by Mecke and Dietrich¹⁹ apart from a gravity term that was included in their expression. To cast ⌬⍀ in the form of Eq.共11兲, we take the Fourier trans- form. One then finds for␴共q兲 共Ref.19兲

␴共q兲 =

冕

−⬁

⬁

dz₁

冕

−⬁

⬁

dz₁₂

冋

^{␻共q,z12兲 −q2␻0共z12兲}

册

⫻关␳0⬘共z1兲␳0⬘共z2兲 − q²␳1共z1兲␳0⬘共z2兲兴 +q²

4

冕

−⬁

⬁

dz₁

冕

−⬁

⬁

dz₁₂␻共q,z12兲␳1共z1兲␳1共z2兲

+q²

4

冕

_−⬁^{⬁ dz1ghs⬙共␳0兲␳1共z1兲2. 共44兲}

Here we have defined the共parallel兲 Fourier transform of the interaction potential

␻共q,z12兲 ⬅

冕

^{drជ储e−iqជ·rជ储Uatt共r兲 = 2␲}

冕

0

⬁

dr_储r_储J₀共qr储兲Uatt共r兲.

共45兲 As a first step, we determine the leading contribution to

␴共q兲 given by the surface tension of the planar interface,␴

=␴共q=0兲. Then, one needs to consider the two leading con- tributions in the expansion of␻共q,z12兲 in q²as follows:

␻共q,z12兲 =␻0共z12兲 +␻2共z12兲q²+ ¯ , 共46兲 where

␻0共z12兲 ⬅

冕

^{drជ储Uatt共r兲,}

␻2共z12兲 ⬅ −1

4

冕

^{drជ储r储2Uatt共r兲. 共47兲}

The surface tension thus becomes

␴⁼

冕

_−⬁^{⬁ dz1}

冕

_−⬁^{⬁ dz12␻2共z12兲␳0⬘共z1兲␳0⬘共z2兲. 共48兲}

The 共planar兲 density profile␳0共z兲, featured in the above ex- pression for ␴, is determined from minimizing the free en- ergy functional ⍀关␳兴 in Eq. 共39兲 in planar symmetry. The Euler–Lagrange equation that minimizes ⍀关␳兴 is then given by

g_hs⬘共␳0兲 = −

冕

^{drជ12Uatt共r兲␳0共z2兲 = −}

冕

_−⬁^{⬁ dz12␻0共z12兲␳0共z2兲,}

共49兲 which can be solved explicitly to obtain␳0共z兲 and thus ␴^.

The evaluation of further contributions to␴共q兲 requires one to determine the density profile␳1共z兲. Just like␳0共z兲, one would like to determine the density profile␳1共z兲 from a mini- mization procedure. One then has to determine the energeti- cally most favorable density profile for a given curvature of the surface.²¹ This turns out to be not so straightforward, since one then has to specify in what way the curvature is set to its given value. Several approaches have been suggested, which we shall now discuss.

• Mecke and Dietrich approach. In this approach a cer- tain form for␳1共z兲 is directly hypothesized,¹⁹

␳1

MD共z兲 = −CH

2␲^⌬␳␰^fH共z/␰兲, 共50兲

with ␰ ^{the bulk} correlation length and fH共x兲

⬅x sinh共x/2兲/cosh²共x/2兲. The coefficient CH in this expression can be used as a fit parameter. This practical approach is certainly legitimate, but one would like to also be able to formulate a molecular basis for this ex- pression.

• Equilibrium approach. Rather than the surface being curved by surface fluctuations, in this approach the in- terface is curved by changing the value of the chemical potential to a value off coexistence. One then considers the density profile of a spherically or cylindrically shaped liquid droplet in metastable equilibrium with a bulk vapor.^20,25This approach is equivalent to adding an external field to the free energy

⍀⬘^关␳兴 = ⍀关␳兴 +

冕

^{drជVext共rជ兲␳共rជ兲 = ⍀关␳兴 −}

冕

^{drជ⌬␮␳共rជ兲,}

共51兲

where⌬␮⁼␮⁻␮coex. The downside of the equilibrium approach is that the external field V_ext共rជ兲=−⌬␮ ^{is uni-} form throughout the system and thus also affects the bulk densities far from the interfacial region. This seems inappropriate for the description of the density fluctuations considered here since we have that ␮

=␮coexand the bulk densities are unaltered by the cur- vature of the surface fluctuations.

• Local external field. In this approach, one again adds to the free energy an external field, but, to ensure that the bulk regions are unaffected, one assumes that it is peaked infinitely sharply at z = 0 as follows:^21,22

V_ext共rជ兲 = ␭␦共z兲⌬h共rជ储兲. 共52兲 In this case, the external field only acts as a Lagrange multiplier in the minimization procedure to ensure that the curvature⌬h共rជ储兲 is set to a certain value; it is not included in the expression for the free energy. The downside of this method is that the resulting density profile ␳1共z兲 has a discontinuous first derivative at z

= 0, which is, from a physical point of view, not so appealing.³⁰ Furthermore, the discontinuous nature of

␳1共z兲 prohibits an analytical simplification using a gra- dient expansion.

In the present approach, we suggest to add an external field acting as a Lagrange multiplier that is unequal to zero only in the interfacial region 共the bulk densities are unaf- fected兲, but which is not infinitely sharp peaked. It seems natural to choose a peak width of the order of the thickness of the interfacial region. It thus seems convenient to choose V_ext共rជ兲⬀␳0⬘共z兲 as follows:

V_ext共rជ兲 = ␭␳0⬘共z兲⌬h共rជ储兲. 共53兲 This choice for V_ext共rជ兲 constitutes our fundamental “ansatz”

for the determination of␳1共z兲. The Lagrange multiplier ␭ is not a free parameter but set by the imposed curvature, as demonstrated below.

The addition of an external field to the free energy re- sults in the following Euler–Lagrange equation:

g_hs⬘共␳兲 = −

冕

^{drជ12Uatt共r兲␳共rជ2兲 − Vext共rជ兲. 共54兲}

Using the external field given in Eq.共53兲, we insert the fluc- tuating density given by Eq. 共9兲 into the above Euler–

Lagrange equation. In order for the resulting equation to hold independently of the value of h共rជ储兲 or ⌬h共rជ储兲, one finds, be- sides Eq.共49兲, the following equation to determine␳1共z兲:

g_hs⬙共␳0兲␳1共z1兲 = −

冕

_−⬁^{⬁ dz12␻0共z12兲␳1共z2兲}

+ 2

冕

−⬁

⬁

dz₁₂␻2共z12兲␳⬘0共z2兲 + 2␭␳0⬘共z1兲.

共55兲 The value of the Lagrange multiplier can be determined by multiplying both sides of the above expression by␳0⬘共z1兲 and integrating over z₁as follows:

␭ = −␴

冒 ^{冋 冕}

^{dz关␳0⬘共z兲兴2}

^册

^{. 共56兲}

One may now verify that if ␳1共z兲 is a particular solution of Eq. 共55兲 then␳1共z兲+␣␳0⬘共z兲 is a solution on account of Eq.

共49兲.

It is convenient to use the Euler–Lagrange equation in Eq. 共55兲to remove the explicit appearance of g_hs⬙共␳0兲 in the expression for␴共q兲 in Eq.共44兲. The resulting␴共q兲 is written as the sum of a term that depends only on the density profile

␳0共z兲 and one term that also depends on the density profile

␳1共z兲 as follows:

␴共q兲 =␴0共q兲 + k1q²+O共q⁴兲, 共57兲 with

␴0共q兲 ⬅

冕

−⬁

⬁

dz₁

冕

−⬁

⬁

dz₁₂

冋

^{␻共q,z12兲 −q2␻0共z12兲}

册

^{␳0⬘共z1兲␳0⬘共z2兲,}

k₁⬅ −1

2

冕

_−⬁^{⬁ dz1}

冕

_−⬁^{⬁ dz12␻2共z12兲␳1共z1兲␳0⬘共z2兲}

+␭

2

冕

_−⬁^{⬁ dz1␳1共z1兲␳0⬘共z1兲. 共58兲}

With the above expression for k₁ it is now also possible to verify that when the density profile␳1共z兲 is shifted by a fac- tor␣␳0⬘共z兲, the resulting effect on the bending parameters is that given by Eq.共22兲.

The procedure to determine␴共q兲, and therefore␴^{, k, k}s, andᐉk, is now as follows: Assuming a certain form for the attractive part of the interaction potential, ␳0共z兲 is obtained from solving Eq.共49兲, which is then inserted into Eq.共55兲to solve for ␳1共z兲 explicitly. The two density profiles thus ob- tained are inserted into Eq. 共58兲to yield␴0共q兲 and k1. This procedure is carried out in Secs. V and VI considering SR forces and LR forces 共U共r兲⬀1/r⁶兲. In general the density profiles ␳0共z兲 and ␳1共z兲 need to be determined numerically.

We shall, however, also provide an approximation scheme, based on the gradient expansion, that is exact near the critical point, but which also gives an excellent approximation far from it.

A. Gradient expansion

The gradient approximation⁷is based on the assumption that the spatial variation of the density profile is small, i.e.,