Evidence for the role of fluctuations in the thermodynamics of nanoscale drops and the implications in computations of the surface tension

Evidence for the role of fluctuations in the thermodynamics of nanoscale drops and the implications in computations of the surface tension

José G. Sampayo    Alexandr Malijevský    Erich A. Müller    Enrique de Miguel    George Jackson111Corresponding author: g.jackson@imperial.ac.uk Department of Chemical Engineering, Imperial College London, South Kensington, London SW7 2AZ, UK
Departamento de Física Aplicada, Facultad de Ciencias Experimentales, Universidad de Huelva, 21071 Huelva, Spain
E. Hála Laboratory of Thermodynamics, Institute of Chemical Process Fundamentals of ASCR, 16502 Prague 6, Czech Republic

Test area deformations are used to analyse vapour-liquid interfaces of Lennard-Jones particles by molecular dynamics simulation. For planar vapour-liquid interfaces the change in free energy is captured by the average of the corresponding change in energy, the leading-order contribution. This is consistent with the commonly used mechanical (pressure tensor) route for the surface tension. By contrast for liquid drops one finds a large second-order contribution associated with fluctuations in energy. Both the first- and second-order terms make comparable contributions, invalidating the mechanical relation for the surface tension of small drops. The latter is seen to increase above the planar value for drop radii of particle diameters, followed by an apparent weak maximum and slow decay to the planar limit, consistent with a small negative Tolman length.

Density functional theory, fundamental measure theory, Monte Carlo and molecular dynamics simulation, Lennard-Jones fluid, interfaces, vapour-liquid equilibria, drops, nucleation, finite-size effects, Laplace relation, Kelvin equation, Tolman length.
Valid PACS appear here

It is striking that though almost a century has passed since Gibbs formulated his thermodynamic theory of curved interfaces, there is still widespread controversy about the dependence of the surface tension on the curvature (size of a drop), and the validity of the mechanical route to the surface tension Rowlinson and Widom (1982); Schofield and Henderson (1982); Croxton (1986). The formal approach of Gibbs is intimately connected with the relations of Laplace, , and Tolman, , for drops of radius . Here, is the pressure difference inside (l) and outside (g) the drop, is the interfacial tension associated with the surface of tension , the value for the planar gas-liquid surface, and the Tolman Tolman (1949) length is defined relative to the radius of the equimolar surface as .

There are 3 basic routes to the definition of the tension Rowlinson and Widom (1982): thermodynamic (Gibbs and Tolman), mechanical (Laplace and Young), and statistical mechanical (density functional and related theories). The thermodynamic and mechanical routes are macroscopic theories so there has been much debate about their applicability to small systems such as nanoscale liquid drops or bubbles. One key question is whether the mechanical relations based on the pressure virial (formulated in terms of the appropriate tensorial components) that make use of the concept of the bulk pressure of the coexisting states are appropriate at these length scales for curved surfaces.

While the Laplace equation essentially defines the ratio , the first-order form of Tolman’s theory is appropriate only for sufficiently large drops. One can view as the leading-order correction to the tension of a planar surface. Despite its fundamental role in studies of interfacial properties of curved surfaces and theories of nucleation, there is still much controversy as to even the sign of . Microscopic statistical mechanical approaches including square gradient theories (SGTs) Falls et al. (1981); Guermeur et al. (1985), curvature expansions of the planar interface Blokhuis and Bedeaux (1992, 1992), and density functional theories (DFTs), including local Lee et al. (1986); Oxtoby and Evans (1988); Koga et al. (1998) and non-local Bykov and Zeng (2002); Li and Wu (2008); Malijevský and Jackson (2009) treatments, have led to conflicting views on the magnitude and sign of , as well as the curvature dependence of the surface tension. The widely accepted view from this body of work is that and that there is a small maximum in as the drop radius is decreased, then followed by a sharp decrease. This is supported by studies on the penetrable sphere model Hemingway et al. (1981) (which can be solved exactly at the mean-field level at zero temperature) where one finds a negative Tolman length (), with the molecular diameter.

By contrast, the vast majority of computer simulation studies suggest that . In most simulations of liquid drops, the mechanical route to the interfacial tension is employed, usually involving an integration of the gradient of the normal component of the pressure tensor from the centre of the drop to the bulk vapour phase Thompson et al. (1984); Nijmeijer et al. (1992); Lei et al. (2005); Vrabec et al. (2006). In this case one predicts a monotonous decrease in the surface tension with increasing curvature (decreasing drop radius) from the planar limit (infinite radius); this would correspond to throughout. As was pointed out early on by Schofield and Henderson Schofield and Henderson (1982) there are fundamental problems in employing local pressure tensors and the associated definition of the internal pressure for microscopic (high curvature) drops. This leads to a mismatch in the free energy of the formation of a drop determined via the mechanical and thermodynamic routes as observed in simulation ten Wolde and Frenkel (1998). Macroscopic thermodynamic routes based on a combination of the Laplace and Tolman relations have been employed Thompson et al. (1984) but also suffer from the ill-definition of the internal pressure and density of the liquid. One can also estimate the interfacial tension from the free energy change accompanying a volume deformation of spherical surfaces using a virial-like expression El Bardouni et al. (2000); these results for the surface tension are in disagreement with those obtained from the direct mechanical route. Recent grand canonical simulations Schrader et al. (2009), and a thermodynamic analysis of large drops based on Laplace-Tolman relations van Giessen and Blokhuis (2009) both now appear to suggest a small negative , which is consistent with the findings of DFT.

The aim of this Letter is to use a new method for the calculation of the surface tension of small liquid drops in molecular simulation, highlighting the role played by the fluctuations in the energy of deformation. The method relies on the thermodynamic definition of the surface tension and is thus free from the inconsistencies associated with the application of the mechanical route. A variant of the test-area (TA) method Gloor et al. (2005) is used where small virtual perturbations are made in the box dimensions of systems with interfaces to obtain the change in free energy associated with the corresponding change in surface area. For a fluid drop of radius the change in the Helmholtz free energy is expressed thermodynamically as Rowlinson and Widom (1982)


with the entropy, the vapour and liquid volumes, the temperature, the chemical potential, the number of particles, the interface area, and the conjugate variable for . The surface tension of a drop is given by


where the minimal interfacial tension defines and corresponds to taking . The change in free energy due to a virtual change in area can be expressed as the average of the Boltzmann factor of the corresponding change in configurational energy  Gloor et al. (2005)


The averages are over configurations of the unperturbed reference system. In Eq. (Evidence for the role of fluctuations in the thermodynamics of nanoscale drops and the implications in computations of the surface tension) is expressed as a perturbation series to , where the first-order average of the change in energy is , the second-order energy fluctuation term is , and the third-order contribution is denoted by . The full Boltzmann form, Eq. (3), is employed in, e.g., the test-particle approach for the chemical potential Widom (1963), or the volume perturbation method for the pressure Eppenga (1994) and the pressure tensor deMiguel2006 (2006).

Figure 1: Test-area deformations of a planar liquid-vapour interface of the LJ-TS fluid. MD simulations of particles in a periodic box of dimensions and at over timesteps. The deformations correspond to changes in the box dimensions (particle coordinates) of , , and . a) The contributions and to the change in free energy per unit area (in units of (, ; , ; average, ). The interfacial tension is obtained by extrapolation to . b) The distribution of the change in energy (relative to its average in units of ) scaled at the maximum peak height for different relative deformations . The width (standard deviation, ) is depicted in the inset.

The tension is obtained as the change in free energy per unit area for infinitesimal perturbations to :


The leading term, , corresponds to the mechanical work involved in changing the area of the interface, which can be directly associated with the so-called virial expression for the tension Lekner and Henderson (1977) (expressed in terms of averages of the appropriate components of the virial, , at the Hookian linear-response level). The corresponding entropic contribution due to the deformation is Lekner and Henderson (1977): .

In the case of the a planar interface, it is well known that the interfacial tension can be obtained formally from the virial expression Rowlinson and Widom (1982); Lekner and Henderson (1977), i.e., entirely from the leading-order contribution of Eq. (5). This is exemplified for a planar vapour-liquid interface of Lennard-Jones particles (of diameter and well depth , truncated and shifted TS at ) as shown in Fig. 1. A planar interface is first stabilised during an molecular dynamics (MD) simulation of the inhomogeneous system with a liquid slab in the centre of the box separated by two vapour regions. The change in configurational energy due to small test changes in the dimensions of the box such that the area is increased or decreased at fixed overall volume is then computed to estimate the various contributions in Eq. (5); the limit of infinitesimal deformations is obtained by extrapolation to . From Fig. 1a it is clear that only the leading mechanical term in contributes to the interfacial tension of a planar interface, confirming the validity of the pressure-tensor route in this case. The fluctuation term is very small by comparison and does not contribute to the tension in an appreciable way; this is also true for the third-order term. In Fig. 1b we plot the distribution of the change in configurational energy (relative to ) for different area perturbations; the distribution is well represented by a Gaussian, the width of which () decreases to zero with , consistent with a very small .

The overall physical picture is fundamentally different for a nanosized spherical drop of liquid in contact with its vapour. Once the drop has been stabilised its size can be characterised from the density profile as a function of the distance from its centre by calculating the Gibbs dividing surface , corresponding to an area of . Virtual perturbations from the equilibrium spherical drop geometry are made with test changes in the dimensions of the simulation cube: two of the Cartesian axes are decreased (or increased) in length and the third is increased (or decreased) such that the overall volume remains constant. The perturbed states correspond to ellipsoidal drops of prolate (or oblate) shape which always have larger surface areas than the original drop, . This essentially corresponds to the longest (Legendre polynomial) capillary-wave oscillations possible for the drop Croxton (1986); the capillary-wave surface tension is equivalent to the thermodynamic one at least to leading order in curvature . Averages are then accumulated over very long runs of timesteps, corresponding to microsecond runs for typical molecular parameters. The term is more than two orders of magnitude larger in the case of the drop than for a planar interface systems of comparable size (cf. Fig. 1a and Fig. 2a). The most significant difference is the large contribution from the second-order energy “fluctuation” term for the drop, which was negligible for the planar interface; this term is now comparable in magnitude but of opposite sign to the first-order term. The third-order terms remains essentially negligible. Thus both the first- and second-order terms contribute to the surface tension of the drop. A thermodynamic characteristic of the drop is thus the non-vanishing (and large) fluctuation term, which is clearly an indication of an additional entropic contribution. This can be seen in the distribution of the change in configurational energy for different TA perturbations (Fig. 2b). The data are again well described as Gaussians, but now the width does not vanish in the limit of infinitesimal deformations, i.e., and correspondingly . The fact that for both the planar and curved systems suggests symmetrical Gaussians.

Figure 2: Test-area ellipsoidal deformations of a spherical drop of LJ-TS liquid of radius in coexistence with its vapour. MD simulations of particles in a periodic box of dimensions at over timesteps. The deformations correspond to changes in the box dimensions (particle coordinates) of , , and (where , , denotes any of the Cartesian axes). a) The contributions , , and to the change in free energy per unit area (in units of (prolate, ; oblate, ; average, ). The tension is obtained by extrapolation to . b) The distribution of the change in energy scaled at the maximum peak height for different relative deformations ; the width is depicted in the inset. The Gaussians for the planar interface are shown dotted (note the very small scale in comparison).

Figure 3: The surface tension of spherical drops of LJ-STS fluids with average radii , 6.12, 8.04, 11.7 and 14.6 at from TA ellipsoidal deformations (), compared with the values from the mechanical route Vrabec et al. (2006) (), and the data of Thomson et al. Thompson et al. (1984) (), El Bardouni et al. El Bardouni et al. (2000) (), and Schrader et al. Schrader et al. (2009) (continuous curve); the planar limit is shown dotted. The predictions of FMT Malijevský and Jackson (2009) are depicted in the inset (dashed curve).

The dependence of the surface tension computed from ellipsoidal deformations as a function of the drop size (for systems with to 11,334) is depicted in Fig. 3. Here the tension is computed for rather than , though to . The behaviour obtained with our thermodynamic TA approach does not support the findings obtained from a standard pressure-tensor route (e.g., the recent MD data of Vrabec et al. Vrabec et al. (2006)). This is in line with the concerns of Schofield and Henderson Schofield and Henderson (1982), and others Blokhuis and Bedeaux (1992); ten Wolde and Frenkel (1998); Reiss and Reguera (2004) about the inadequacy of the mechanical route for very small systems. For drops larger than we observe values of the surface tension which appear to be slightly larger than the planar limit, ; because the tension has to converge to when this suggests a non-monotonic behaviour of the tension with increasing curvature, and a corresponding weak maximum. Our values are consistent with the data point reported by El Bardouni et al. El Bardouni et al. (2000) estimated from the surface free energy change, and with the small maximum observed by Schrader et al. Schrader et al. (2009) using a Landau free energy approach in the canonical ensemble (though the authors do not comment explicitly on this point). The calculations of the tension of curved interfaces from curvature corrections, SGT and DFT (which have been brought into question because of their failure to reproduce existing simulation data) are now supported by our data. In the inset of Fig. 3 we compare the TA data for the surface tension of drops with those from a non-local DFT using fundamental measure theory (FMT) Malijevský and Jackson (2009); a maximum is predicted with FMT at .

Three main conclusions can be gleaned from our study. Firstly, there is clearly a large fluctuation contribution to the interfacial tension of nanoscale spherical drops (and most likely other curved surfaces) in addition to the underlying first-order (mechanical) contribution which fully describes the planar interface. Such contributions from fluctuations in the energy are not found in the planar limit to any significant degree. As a consequence, our results do not support the validity of a mechanical (pressure-tensor) route to the interfacial tension for surfaces of high curvature such as small drops. This is in line with the warning of Blokhuis and Bedeaux Blokhuis and Bedeaux (1992) that the use of a mechanical approach in this context “is still a matter of concern” and that it is “advisable not to use the pressure-tensor whenever this can be avoided”. Our data are not consistent with the monotonic dependence of the surface tension with curvature obtained from a mechanical treatment. As well as contributions in , the correct “virial” expression for the surface tension would have to contain terms in averages of the type and , which would involve up to four-body correlations for pair-wise additive potentials. This suggests that there are additional contributions to the change in the entropy due to the deformation of small drops involving quadratic terms in : , , , , , and . As a final point, the rise of the surface tension above that of the planar limit after a certain drop size would be consistent with a negative Tolman length. Our data for the larger drops suggests a value of . Though the statistical uncertainty is large, our finding supports the exact mean-field predictions for the penetrable sphere model Hemingway et al. (1981), and is consistent with the latest accurate value of determined from the Laplace relation for a large particle system van Giessen and Blokhuis (2009).

We are very appreciative to Jim Henderson for useful discussions. Financial support from the Czech, Mexican, Spanish and UK councils is gratefully acknowledged.


  • Rowlinson and Widom (1982) J. S. Rowlinson and B. Widom, Molecular Theory of Capillarity (Oxford University Press, Oxford, 1982)
  • Schofield and Henderson (1982) P. Schofield and J. R. Henderson, Proc. Roy. Soc. A379, 231 (1982).
  • Croxton (1986) J. R. Henderson Statistical Mechanics of Spherical Interfaces in Fluid Interfacial Phenomena, edited by C. A. Croxton (John Wiley & Sons Ltd., New York, 1986)
  • Tolman (1949) R. C. Tolman, J. Chem. Phys. 17, 333 (1949).
  • Falls et al. (1981) A.H. Falls, L.E. Scriven, H.T. Davis, J. Chem. Phys. 75, 3986 (1981).
  • Guermeur et al. (1985) R. Guermeur, F. Biquard, C. Jacolin, J. Chem. Phys. 82, 2040 (1985).
  • Blokhuis and Bedeaux (1992) E. M. Blokhuis, and D. Bedeaux, Physica A184, 42 (1992).
  • Blokhuis and Bedeaux (1992) E. M. Blokhuis, and D. Bedeaux, J. Chem. Phys. 97, 3576 (1992).
  • Lee et al. (1986) D. J. Lee, M. M. Telo da Gama, and K. E. Gubbins, J. Chem. Phys. 85, 490 (1986).
  • Oxtoby and Evans (1988) D.W. Oxtoby, R. Evans, J. Chem. Phys. 89, 7521 (1988).
  • Koga et al. (1998) K. Koga, X. C. Zeng, A. K. Shchekin, J. Chem. Phys. 109, 4063 (1998).
  • Bykov and Zeng (2002) T.V. Bykov, X.C. Zeng, J. Chem. Phys. 117, 1851 (2002).
  • Li and Wu (2008) Z. Li, and J. Wu, Ind. Eng. Chem. Res. 47, 4988 (2008).
  • Malijevský and Jackson (2009) A. Malijevský, and G. Jackson, to be submitted (2010).
  • Hemingway et al. (1981) S. J. Hemingway, J. R. Henderson, and J.S. Rowlinson, Faraday Symp. Chem. Soc. 16, 33 (1981).
  • Thompson et al. (1984) S.M. Thompson et al., J. Chem. Phys. 81, 530 (1984).
  • Nijmeijer et al. (1992) M. J. P. Nijmeijer et al., J. Chem. Phys. 96, 565 (1992).
  • Lei et al. (2005) Y. A. Lei, T. Bykov, S. Yoo, and X. C. Zeng, J. Am. Chem. Soc. 127, 15346 (2005).
  • Vrabec et al. (2006) J. Vrabec, G. K. Kedia, G. Fuchs, and H. Hasse, Mol. Phys. 104, 1509 (2006).
  • ten Wolde and Frenkel (1998) P. R. ten Wolde, and D. Frenkel, J. Chem. Phys. 109, 9901 (1998).
  • El Bardouni et al. (2000) H. El Bardouni, M. Mareschal, R. Lovett, and M. Baus, J. Chem. Phys. 113, 9804 (2000).
  • Schrader et al. (2009) M. Schrader, P. Virnau, K. Binder, Phys. Rev. E 79, 061104 (2009).
  • van Giessen and Blokhuis (2009) A.E. van Giessen, E.M. Blokhuis, J. Chem. Phys. 131, 164705 (2009).
  • Gloor et al. (2005) G. J. Gloor, G. Jackson, F. J. Blas, and E. de Miguel, J. Chem. Phys. 123, 134703 (2005).
  • Widom (1963) B. Widom, J. Chem. Phys. 39, 2808 (1963).
  • Eppenga (1994) R. Eppenga, and D. Frenkel, Mol. Phys. 52, 1303 (1994).
  • deMiguel2006 (2006) E. de Miguel, and G. Jackson, J. Chem. Phys. 125, 164109 (2006).
  • Lekner and Henderson (1977) J. Lekner, and J. R. Henderson, Mol. Phys. 34, 333 (1977).
  • Reiss and Reguera (2004) H. Reiss, and D. Reguera, J. Phys. Chem. B 108, 6555 (2004).
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description