JOURNAL OF CHEMICAL PHYSICSVOLUME 109, NUMBER 10

8 SEPTEMBER 1998

Validity of Tolman’s equation: How large should a droplet be?
Kenichiro Kogaa) and X. C. Zeng
Department of Chemistry and Center for Materials Research and Analysis, University of Nebraska, Lincoln,
Nebraska 68588

A. K. Shchekin
Department of Statistical Physics, St. Petersburg State University, 198904, St. Petersburg, Russia~Received 23 February 1998; accepted 2 June 1998!
Surface tension and the length d ~distance between the Gibbs surface of tension R s and theequimolar surface R e ! of simple liquid droplet ~Lennard-Jones and Yukawa! are computed over awide range of droplet sizes up to about 43106 molecules. The study is based on the Gibbs theoryof capillarity combined with the density–functional approach to gas–liquid nucleation. Since thismethod provides behavior of the surface tension fully consistent with the tension of the planner
surface, the constant in Tolman’s equation d ` can be determined unequivocally from the asymptotic
behavior of s s . Comparison of the tension given by Tolman’s equation against the result of exact
thermodynamic relations reveals that Tolman’s equation is valid only when the droplet holds morethan 106 molecules for the simple fluid systems near their triple points, in contrast to theconventional wisdom that Tolman’s equation may be applicable down to droplets holding a few
hundreds of molecules. © 1998 American Institute of Physics. @S0021-9606~98!50334-0#

I. INTRODUCTION

strong function of R s , i.e., assumption B is invalid for verysmall droplets.13 Size dependence of surface tension andlength d becomes more and more controversial as one considers larger and larger droplets. This is mainly due to thefact that numerical errors in the calculations overwhelm theprecision one needs for large droplets for which dependenceof s s and d on R s is very weak. To our knowledge, noquantitative information has been reported as to a range ofdroplet sizes where we can apply Tolman’s equation for calculating surface tension.
The purpose of this article is to shed some light on towhat extent Tolman’s equation is valid. To this end, we compute the surface tension and the length d for a wide range ofdroplet sizes on the basis of the Gibbs thermodynamic theoryof capillarity combined with the density–functional theory ofgas–liquid nucleation.14,15
We first study behavior of d for sizes of droplet ranging
from 40 to about 4.03106 particles. Two different routes tod, both of which are exact, are taken: a direct route whichcalculates R e and R s separately, and an adsorption route we
propose in Sec. IV. We will show that the direct route gives
rise to unphysical behavior of d when we approach the planner surface limit whereas the adsorption route can lead to
self-consistent and precise results of asymptotic behavior of
the length d. From the asymptotic behavior, the value of
Tolman’s length d ` can be determined accurately. It is thenpossible to examine the range of droplet sizes where Tolman’s equation is valid by comparing results of Tolman’sequation ~using d ` ! with those from rigorous calculations
based on the thermodynamic definition of surface tension.
We also discuss the validity of other thermodynamic
routes ~mostly used in computer simulations! to surface tension, length d, and Tolman’s length d ` . One such route is
the thermodynamic route proposed by Thompson et al.10 in
the molecular dynamics simulation study of Lennard-Jones

One fundamental question on properties of droplets ofliquid is the dependence of surface tension on droplet radius.
Numerous theoretical investigations have been devoted tothis question, including the original thermodynamic consideration by Gibbs,1 later development by Tolman,2 and recent
studies based on statistical mechanical theories3–8 and computer simulations.9–12
When a droplet is sufficiently large, size dependence ofthe surface tension can be expressed by Tolman’s equationof the form2

ss
Rs
5
.
s ` R s 12 d `~1!

Here, s s is the surface tension with respect to the Gibbs
surface of tension and R s denotes the radius of this dividing
surface; s ` is the planner surface tension, and d ` is defined
as limR s →` d , where d is the distance between the equimolar
surface with radius R e and the surface of tension. Tolman
derived Eq. ~1! based on the Gibbs theory of capillarity1 plus
two additional assumptions: d !R s ~assumption A! and d
5 d ` ~assumption B!.
On the other hand, in the case of small droplets, no analytical relation has been found so far to express s s in termsof R s . There have been several numerical calculations of s s
for various sizes of droplets via the gradient theories,3,4,6 thedensity–functional theory,7 and computer simulations.9–12
Most studies reached a common conclusion that very rapid
decrease in surface tension takes place with decrease in droplet size when R s is sufficiently small ~e.g., smaller than about
ten molecular radii3!. Some studies3,7 also showed that d is aa!

Electronic mail: koga@phase2a.unl.edu

0021-9606/98/109(10)/4063/8/$15.00

4063© 1998 American Institute of Physics

4064

J. Chem. Phys., Vol. 109, No. 10, 8 September 1998

Koga, Zeng, and Shchekin~LJ! liquid drops. Recently, Nijmeijer et al.11 adopted another thermodynamic route to d in large scale computer
simulations. Both of the routes seem particularly convenient
for the computer simulation because they require neither theequation of state of the bulk fluid nor the grand potential ofthe droplet–vapor system. We note, however, that theseroutes are based on several approximations to the rigorous
thermodynamic theory and become exact only in the large
droplet limit. We therefore also examined the extent towhich these routes are applicable to the determination of s s
and d.
This paper is organized as follows: Section II gives
Gibbs thermodynamic description of the droplet–vapor system. Section III summarizes the density–functional theory ofgas–liquid nucleation. Results and discussion are given inSec. IV and concluding remarks are given in Sec. V.
II. THERMODYNAMIC DESCRIPTION OF THE
DROPLET–VAPOR SYSTEM

We consider a critical droplet of a single-component liquid in ~unstable! equilibrium with a supersaturated vapor attemperature T. The system is taken to be a conical section~with solid angle v! of a sphere ~with radius R! which includes the droplet at the center with surrounding vapor. The
radius R is chosen to be sufficiently large such that thevapor at the boundary acquires the bulk properties. The
chemical potential m of the system is higher than the chemical potential m eq at vapor–liquid ~stable! equilibrium at T.
The reversible transformations of the system are describedby the Gibbs fundamental relation
dU5T dS2dW1 m dN,~2!

where the symbols U, S and N denote, respectively, the energy, entropy and number of molecules of the system. The
work performed by the system is represented by dW. To
define the surface or excess variables including surface tension, one needs to introduce a mathematical dividing surface
with an arbitrary radius R, which partitions the total volume(V) into two parts: one inside the dividing surface V a andthe other part V b . Then a hypothetical reference system is
introduced, which possesses the same m, T, and V but iscomposed of the uniform liquid and vapor phases ~let us callthem a and b phases, respectively! right up to the dividing
surface; namely, the reference system consists of two parts: aportion (V a ) of a bulk liquid at m and T and a portion (V b )of a bulk gas at m and T. Since m . m eq , phase a is thermodynamically stable, whereas phase b is metastable. Now thesurface variables can be defined as
N s 5N2N a 2N b ,~3!

U s 5U2U a 2U b ,

~4!

a

~5!

b

S 5S2S 2S ,swhere the superscripts a and b indicate the quantities of eachpart of the reference system. The infinitesimal change in theinternal energy of each part of the reference system is described by the fundamental relationdU a 5T dS a 2 p a dV a 1 m dN a ,~6!

dU b 5T dS b 2 p b dV b 1 m dN b .

~7!

a

It is important to emphasize that p is not the pressure p(0)at the center of the droplet but is the pressure of homogeneous phase a in the reference system; in general, these twovalues are different. The change in the surface energy isassumed to be given bydU s 5T dS s 1 s dA1C dR1 m dN s ,~8!

where A is the surface area of the dividing surface and s is
the surface tension defined as a conjugate variable of A; C isthe conjugate variable of R. Combining Eqs. ~6!–~8! results
in a relation
dU5T dS2p a dV a 2 p b dV b 1 s dA1C dR1 m dN.
~9!
Comparing Eqs. ~9! with ~2!, one can find that s dA
1C dR is the difference between the work done on the realsystem (2dW) and the work done on the hypothetical reference system (2p a dV a 2 p b dV b ). Integrating Eq. ~9!
with respect to the solid angle v at fixed T, m and R yields5

V1 p a V a 1 p b V b DV1D pV a
5
,
A
A

~10!

where V[U2TS2 m N is the grand potential, DV[V1 p b V is the work of formation of the critical droplet, andDp stands for p a 2 p b . Equation ~10! holds for surface tension with respect to any dividing surface. Note that s depends not only on the thermodynamic state of the system but
also on the choice of the parameter R. However, if a particular dividing surface such as the surface of tension or theequimolar surface is chosen, the corresponding surface tension turns out to be a state function.
The surface of tension with radius R s is a special dividing surface for which the coefficient C in Eq. ~9! vanishes
and the Laplace equation holds:16
Dp52 s s /R s .

~11!

Recalling that V a 54 p R 3s /3 and A54 p R 2s , one can solve thecoupled equations ~10! and ~11! to obtain

s s5

R s5

and

S

3DV ~ Dp ! 2 1/3
,
16p

D

~12!

S

3DV 1/3.
2 p DpD~13!

These are rigorous thermodynamic expressions for s s andR s . One may question the validity of thermodynamic arguments when a droplet is too small to possess bulk liquid
properties inside its surface. However, as noted by Gibbs,1
his original thermodynamic arguments ~presented above! do
not assume any thermodynamic property inside the droplet;
instead they rely upon the properties of the homogeneous
liquid phase having the same chemical potential as thedroplet–vapor system has. Whether or not there is a homogeneous liquid phase inside the droplet, the thermodynamic
quantities of the reference system are well defined. Therefore, the surface tension and other surface quantities are well

J. Chem. Phys., Vol. 109, No. 10, 8 September 1998

Koga, Zeng, and Shchekin

defined no matter how small the droplet is. In this sense, theGibbs thermodynamic arguments hold for any size of a droplet. On the other hand, any thermodynamic arguments relying upon thermodynamic quantities in the actual droplet~e.g., pressure at the center of the droplet! become less andless reliable when we go to smaller and smaller droplets.
Tolman’s equation is derived as follows. From the
Gibbs–Duhem equation for the surface variables (A d s s
1N s d m 1S s dT50) together with the definition of the surface adsorption (G s [N s /A), we haveG s 52 ~ ]s s / ] m ! T .

~14!

Combining a similar relation for the reference system, D r[ r a 2 r b 5( ] Dp/ ] m ) T , with Eq. ~14! yields a differential
equation:

S D

Gs
]s s
52
.
] Dp TDr~15!

In the planner surface limit, the quantity G s /D r on the righthand side of Eq. ~15! is identical to Tolman’s length d ` .
This is seen from an exact relation ~originally derived byTolman2!

SD

d 1 d2Gs
5 d 11 1.
Dr
R s 3 R 2s

~16!

If we integrate Eq. ~15! with respect to D p from the planar
surface limit (Dp50), assuming that G s /D r is constantover the range, we finds s 2 s ` 52 d ` D p.

~17!

Substitution of Laplace equation ~11! into Eq. ~17! results in
Tolman’s equation ~1!.

III. APPLICATION OF THE DENSITY–FUNCTIONAL
THEORY TO THE DROPLET–VAPOR SYSTEM
OF SIMPLE FLUIDS

Although Eqs. ~12! and ~13! are rigorous thermodynamic
expressions for the surface tension s s and the radius R s ofthe surface of tension, the determination of s s and R s requires values of DV ~the work of droplet formation! and Dp.
Thermodynamics does not provide these data for the
droplet–vapor system; a molecular-based approach ~statistical mechanics or computer simulation! is needed to determine them. The density–functional theory of nucleation issuch an approach.14,15 This approach has the advantage thatmolecular-level detail can be incorporated, and in the limit oflarge droplets the theory goes naturally to the classical nucleation theory. Effects of curvature dependence of the surface
tension arise naturally, rather than being added as ad hocassumptions. Using the density–functional theory, one canobtain the density profile r~r! of a spherical droplet, which isa critical nucleus of the liquid phase, and also the grandpotential V of the droplet–vapor system.14
The formalism of the density–functional theory of
nucleation can be briefly summarized as follows: The Helmholtz free energy functional is taken to be a simple form17

F @ r ~ r!# 5

E

dr f h @ r ~ r!#

1

1
2

EE

4065

dr dr8 r ~ r! r ~ r8 ! f att~ u r2r8 u ! ,~18!

where f h ( r ) is the Helmholtz free energy density ~the freeenergy per unit volume! of a uniform hard-sphere fluid,which is given by
f h~ r ! 5 r m h~ r ! 2 p h~ r ! .

~19!

The chemical potential m h ( r ) and the pressure p h ( r ) of thereference hard sphere fluid are evaluated using the highly
accurate Carnahan–Starling equation.18 f att in Eq. ~18! is theattractive part of a pair potential.
We consider the Lennard-Jones ~LJ! fluid and Yukawa
fluid in this study. In the case of the LJ fluid, f att is taken to
be the Weeks–Chandler–Anderson ~WCA! perturbation partof the LJ potential function19

FS D S D G

f att~ r ! 54 e LJ

s LJ 12 s LJ 62rr

52 e LJr,2 1/6s LJ .

r>2 1/6s LJ~20!

The two parameters e LJ and s LJ are, respectively, the depthof the potential well and the collision diameter. Reduced
variables @denoted by an asterisk ~*!# are defined in terms ofparameters: distance is given in units of s LJ , energy in unitsof e LJ , and temperature in units of e LJ /k B , where k B is theBoltzmann constant. We chose the LJ droplet–vapor systemsince it has been extensively studied both by computer
simulation9–11 and by the density–functional theory.15 In order to fairly reproduce the liquid density of the LJ fluid attemperatures lower than the critical point, we use thetemperature-dependent hard sphere diameter d(T). 20 In thecase of the Yukawa fluid, the attractive part is taken to be

f att~ r ! 52 a l 3 exp~ 2lr ! /4p lr.

~21!

Furthermore we set l to be d 21 , the inverse of the hard
sphere diameter of the repulsive part. In this system, reduced
variables are defined in terms of d and k B T c : distance isgiven in units of d, densities in units of d 23 , energy in unitsof k B T c , temperature in units of T c .
We note that Eq. ~18! is based on two key approximations: First, the pair distribution function r (2) (r,r8 ) of an
inhomogeneous fluid is replaced by r (r) r (r8 ). This is theso-called random phase approximation.21 Second, the local
density approximation is used to determine the hard spherepart of the free energy functional. These two approximations
have been well tested for weakly inhomogeneous systems
such as a vapor–liquid interface and are fairly good approximations when the system is not as close to the critical point.
The grand potential functional is given as a Legendre
transform of F
V @ r # 5F @ r # 2 m

E r

dr ~ r! .

~22!

4066

J. Chem. Phys., Vol. 109, No. 10, 8 September 1998

Koga, Zeng, and Shchekin

The density profile of a critical nucleus at m and T ~as well asthat of a planar surface at m eq and T! can be determined fromthe equilibrium conditiond V @ r # / d r ~ r! 50,

~23!

or

m 5 m h @ r ~ r!# 1

E 8f
dr

att~ u r2r8 u ! r ~ r8 ! .

~24!

We also note that the density profile of the planar surface
can be obtained from Eq. ~24! using the conventional iteration technique. However, the density profile of a critical
nucleus is not a stable solution of Eq. ~24! since it corresponds to a saddle point in the functional space. Nonetheless,
if an appropriate initial guess is made, one can find the unstable equilibrium density profile which satisfies the condition ~23! over many iterations.15
The density profile and the grand potential of the critical
droplet–vapor system were calculated over a wide range ofsupersaturation D m [ m 2 m eq at a fixed temperature: for the
LJ system T * 50.7, which is close to the triple-point temperature T *tp '0.69 of the system; for the Yukawa system
T/T c 50.6.
IV. RESULTS AND DISCUSSION
A. Behavior of the length d

The usual way of determining the length d is to calculate
R e and R s and then use the definition d [R e 2R s , which wehereafter refer to as ‘‘a direct route to d’’. The radius R e ofthe equimolar surface can be determined from the density
profile of a critical droplet by the exact relation
R 3e 5

3
Dr

E r
R

0

@ ~ r ! 2 r b # r 2 dr.

~25!

Location of the surface of tension, R s , is given by Eq. ~13!.
Another way of determining d, which we propose here,
is to use the relation
2SD

d 1 d2~ ]s s / ] m ! T
5 d 11 1.
Dr
R s 3 R 2s

~26!

This is derived from Eqs. ~14! and ~16!. s s and R s are givenby the exact expressions ~12! and ~13!, respectively. The
derivative ( ]s s / ] m ) T can be calculated via numerical differentiation. We refer to this method as ‘‘an adsorption route tod.’’
We now discuss both the result of the direct route andthat of the adsorption route. Figure 1 shows d of the LJ
droplet as a function of supersaturation Dm. When supersaturation is large (D m * .0.2), the behavior of d obtained fromthe two routes is almost indistinguishable. In this range ddecreases linearly with the decrease in Dm and changes its
sign at about D m * 50.3. The two routes, however, lead toqualitatively different results when D m * ,0.2, that is, in arange of very large droplets; d via direct route turns to bepositive and tends to diverge as Dm goes to zero, whereas dvia the adsorption route still decreases linearly and eventually converges to a small negative value d ` 520.15s LJ in
the planar surface limit. Qualitatively the same result is

FIG. 1. The length d * vs the supersaturation D m * at T * 50.7: the directroute ~dotted line! and the adsorption route ~solid line!. The reduced dimensionless variables are defined in terms of the LJ size and energy parameters.

found for the Yukawa potential system; the direct route gives
rise to the diverging behavior of d, whereas the adsorption
route leads to d ` 520.26d. The diverging behavior of d via
direct route is unphysical because both the equimolar surface(R e ) and the surface of tension (R s ) should lie within thevapor–liquid interface whose width remains finite over awhole range of supersaturation at thermodynamic states near
the triple point.
In principle, both direct and adsorption routes should
lead to the same behavior of d over a whole range of Dm
because no approximations are made in either route ~except
common approximations in the density–functional theory!.
Therefore, the problem with the direct route ~unphysical behavior of d in a range of small Dm! is a numerical one rather
than a theoretical one. The reason is the following. As Dm
goes to zero, R s and R e increase simultaneously and diverge
ultimately. In this process, absolute numerical errors of R s
and R e also increase systematically although relative errors
may remain small @see, for example, expression ~13! for R s ;the denominator Dp goes to zero.# Consequently, relative
errors of d given by subtracting R s from R e are expected tobe large and increase systematically. This problem would beeven more serious for computer simulation studies due to thestatistical errors. In fact, computer simulation results showthat the sign of d is uncertain for large droplets.11 Thus, wemust conclude that exact asymptotic behavior of d in therange of very large droplets cannot be obtained through thedirect route.
On the other hand, the adsorption route is free from thenumerical problem. The reason is the following. This approach requires ( ]s s / ] m ) T , Dr, and R s ; however, in thelarge droplet limit the ratio of the first two quantities determines d @see Eq. ~26!#. Unlike R e and R s , both ( ]s s / ] m ) T
and Dr remain finite values as Dm goes to zero. Thus, theratio ( ]s s / ] m ) T /D r can be determined accurately over awhole range of droplet sizes. As a consequence, asymptotic
behavior of d and the value of d ` can be determined without
suffering from any numerical problem.
Furthermore, it can be seen from the following that the

J. Chem. Phys., Vol. 109, No. 10, 8 September 1998

Koga, Zeng, and Shchekin

4067

FIG. 2. The surface tension ratio s s / s ` vs the supersaturation D m * atT * 50.7. The inset is the magnification at small supersaturation. Dotted linein the inset is drawn by extrapolation.

asymptotic behavior of d obtained by the adsorption-routecalculation is consistent with the properties of the planar
liquid–vapor interface. Figure 2 shows that the surface tension of the spherical surface of the LJ droplets is fully consistent with that of the planar surface; extrapolation showss s (D m 50)/ s ` 51.000 07, where s ` was obtained from thedensity–functional calculation for the planar liquid–vapor
interface. The self-consistency for the values of s s guarantees that for ( ]s s / ] m ) T . In addition, the behavior of Dr is
also self-consistent. Therefore, the asymptotic behavior of dis also consistent with the properties of the liquid–vapor interface within the framework of the density–functional
theory.
B. Range of the validity of Tolman’s equation

Once Tolman’s length d ` is determined accurately, direct comparison can be made between the behavior ofs s (R s ) given by exact thermodynamic relation ~12! and thatby Tolman’s equation ~1!. Figure 3~a! shows a comparison
for the LJ droplet–vapor system. When the reduced radius
R s* is greater than about 50, prediction of Tolman’s equationis almost indistinguishable in the scale of the plots from theresult of the exact thermodynamic formula; the ratios s (R s )/ s ` is slightly larger than unity and increases very
slowly as R s decreases. The same result was obtained in thecase of the Yukawa potential system @Fig. 3~b!#. This confirms that d ` was determined accurately in such a way thatd ` is consistent with the asymptotic behavior of s s .
However, when R s* is smaller than about 20, or N a,106 , prediction of Tolman’s equation becomes significantly different from the result of the rigorous calculation.
The surface tension s s of the exact formula shows nonmonotonic behavior; as R s decreases it reaches a maximum atabout R s* 510 and begins to decrease rapidly. On the other
hand, s s of Tolman’s equation is necessarily a monotonic
function of R s ; Fig. 3~a! shows monotonic increase in s swith decrease in R s due to the negative value of d ` . The
same significant difference was found in this range of droplet
size in the case of the Yukawa potential fluid @see Fig. 3~b!#.

FIG. 3. The surface tension ratio s s / s ` obtained via the rigorous thermodynamic route ~solid line with circles! and the prediction of Tolman’s equation ~solid line!: ~a! Lennard-Jones fluid and ~b! Yukawa fluid. Values of d `~in Tolman’s equation! were obtained via the adsorption route.

These results suggest that Tolman’s equation is valid
only for droplets with size N a >106 , at least near the triple
points of the systems. This is the main conclusion of thisarticle. Since basically the same results are obtained from thetwo distinctive potential functions ~LJ and Yukawa!, it isvery unlikely that replacement of the long-range interaction
by r 27 ~i.e., consideration of the retardation effect in thedispersion force at a very large distance22,23! would give a
qualitatively different result.
We also note that assumption B ( d 5 d ` ) becomes incorrect in a range of droplet sizes where assumption A ( d /R s
!1) is still valid; d is no longer constant when R s* ,20,
whereas d /R s 50.02 at R s* 55 ~see the solid line with points
in Fig. 5!. The breakdown of assumption B is closely related
to the nonmonotonic size dependence of surface tension
since if d changes its sign so does the term ( ]s s / ] m ) T andvice versa @see Eq. ~26!#. Also if Tolman’s length d ` is negative, assumption B should become incorrect in some range ofdroplet sizes since d .0 and ( ]s s / ] m ) T ,0 for very small
droplets. Thus, d ` ,0 is also closely related to the breakdown of assumption B and to the nonmonotonic behavior of
ss .

4068

J. Chem. Phys., Vol. 109, No. 10, 8 September 1998

FIG. 4. The surface tension ratio s s / s ` vs R s* : the rigorous thermodynamic route ~solid line with circles!; the route of Thompson et al. ~dashedline!; and the modified route ~dotted line!.

C. Examination of other thermodynamic routes

As mentioned in Sec. III, rigorous calculation of s s andR s requires the work of formation of a critical droplet DV
and the pressure difference D p as thermodynamic inputs.
Although the latter two quantities can be obtained straightforwardly by the density–functional theory, they are not easily obtained from computer simulation. For example, theequation of state of the bulk fluid is required to determine thepressures ~p a and p b ! of the reference homogeneous system
unless the droplet size is sufficiently large so that p a can be
replaced by the value at the droplet center p(0).
To circumvent the technical difficulties in the simulation, several alternatives to the rigorous thermodynamic
routes have been proposed. Since these thermodynamic
routes involve certain approximations, the extent to which
they are applicable for determining s s as a function of droplet size needs to be examined. One thermodynamic route we
examine here was originally proposed by Thompson et al. in
their molecular dynamics study of microliquid drops.10 In
this route, s s and R s are obtained by solving the two coupledequations: Tolman’s equation ~1! and Laplace’s equation~11!. Approximations in this route are the following: ~i! theuse of Tolman’s equation;24 ~ii! replacement of the constantd ` in Eq. ~1! by a variable d 5R e 2R s ; ~iii! replacement of
Dp in the rigorous thermodynamic relation ~11! by Dp(0);and ~iv! replacement of Dr in the formula ~25! by D r (0)[ r (0)2 r b . Note that all the approximations listed above
turn to be exact in the large droplet limit. This route is particularly convenient for the computer simulation since quantities such as R e , Dp(0), and Dr~0! can be calculated directly in the simulation of a droplet–vapor system. However,in order to apply this method we need to clarify the range ofvalidity for these approximations.
Figure 4 compares behavior of s s obtained via the routeof Thompson et al. with that obtained via the rigorous route
over a wide range of droplet sizes. The difference is significant for small droplets, that is, in the range R s* ,20. Similar
results at different temperatures have been obtained by Ta-

Koga, Zeng, and Shchekin

FIG. 5. The length d * vs R s* : the strict thermodynamic route ~solid linewith circles!; the route of Thompson et al. ~dashed line!; and the modified
route ~dotted line!.

lanquer and Oxtoby.7 We also find from Fig. 4 that even forvery large droplets ~e.g., R s* ;100 where Tolman’s equation
is valid and all the other approximations are totally reasonable! there exists a small but systematic difference in s s .
Again, this is due to the fact that the quantity d is determinedby the subtraction of two large quantities, R e 2R s , whichcauses larger and larger numerical errors as we go to larger
and larger droplets ~see Fig. 5!. The computer simulation ismore prone to this problem because of the statistical errors
involved in averaging.
Let us consider the possibility of improving the thermodynamic route of Thompson et al. Since both the use of Tolman’s equation and the replacement of d ` by d are key approximations in this route, let us focus on approximations~iii! and ~iv!, which can be removed provided the equation ofstate for the particular system is known. It is expected thatthe effect of approximation ~iii! on the results becomes larger
and larger with the decrease in droplet size. This is because
the size dependence of the center pressure Dp(0) is oppositethat of Dp in a range of small droplets; Dp(0) decreases
with the decrease in droplet size, whereas Dp increases
monotonically.
To examine the effects of approximations ~iii! and ~iv!
we calculated s s and the length d after replacing D p(0) andDr~0! by Dp and Dr of the rigorous route. ~We refer to thisapproach as the modified route.! Figures 4 and 5 show theresults of s s and d determined by the modified route, together with the results from the rigorous route and from theroute of Thompson et al. Contrary to our expectation, we
find that results of the modified route are rather worse thanthose of the original route of Thompson et al. except in thelarge droplet limit where both results converge to the samevalues. That is, with the decrease in R s the surface tension
increases monotonically and the values of d remain negative
over a whole range of R s . The fact that removing approximations ~iii! and ~iv! worsens the results of the route ofThompson et al. indicates that although the results obtained
via the route of Thompson et al. bear some similarities tothose obtained via the rigorous route ~e.g., as R s decreases

J. Chem. Phys., Vol. 109, No. 10, 8 September 1998

Koga, Zeng, and Shchekins s decreases for a range of small droplets, and d changes its
sign from negative to positive!, these similarities appear fortuitous.
From these analyses, we conclude that as far as the size
dependence of the surface tension and Tolman’s length d `are concerned, the thermodynamic route proposed byThompson et al. seems unsuitable.
Nijmeijer et al.11 also attempted to determine d ` fromthe molecular dynamics simulation of LJ droplets with reduced radii R s* <13. The route of Nijmeijer et al. to d ` requires the use of the relation
R e Dp ~ 0 !
2d
2252
,
s`
Re

~27!

which entails essentially the same approximations as used in
the route of Thompson et al. Thus the route of Nijmeijer
et al. should encounter the same problem emerged in theroute of Thompson et al., particularly when the droplet issmall. Although Eq. ~27! becomes exact in the large droplet
limit, again this route invokes the subtraction R e 2R s whichcauses large numerical errors for d. The limitation of thisthermodynamic route to d ` has already been pointed out by
Nijmeijer et al.11
V. CONCLUDING REMARKS

Based on the Gibbs thermodynamic theory of capillarity
combined with the density–functional theory of gas–liquid
nucleation, size dependence of the surface tension of LJ and
Yukawa droplets and that of d were investigated over a wide
range of droplet sizes. Tolman’s length d ` , namely d in theplanar surface limit, was determined unequivocally via what
we call the adsorption route to d which relies on an exact
thermodynamic relation ~26!. Direct comparison between result of rigorous calculation based on the thermodynamic relations ~12!–~13! and that from Tolman’s equation showsthat Tolman’s equation is valid only when the droplet holds
more than 106 molecules or has a radius larger than about 20
molecular diameters, which is much larger than droplets
mostly studied so far by computer simulations. We pointed
out that Tolman’s equation becomes invalid mainly due tothe breakdown of assumption B ( d 5 d ` ), rather than A ( d!R s ). As discussed in Sec. IV, the breakdown of assumption B is closely related to the nonmonotonic behavior of s s~as a function of R s or Dm! and the negative value of Tolman’s length. The nonmonotonic behavior has been foundpreviously in the framework of phenomenological gradienttheories4,5 and the square gradient theory with a doubleparabola model for free energy density.6 In the latter model,
Tolman’s length is expressed in terms of correlation lengthsof bulk gas and liquid phases; Iwamatsu concludes from theresult that in general Tolman’s length should be negative( d ` 520.38d at T50.8T c !. Also, recent calculation of d `~as the first-order curvature correction to the surface freeenergy! based on a mean field approximation shows thatd ` '20.20d at all temperatures.25 Furthermore it is shown26
that close results are obtained from three different formulas~two of Blokhuis and Bedeaux27 and one of Bykov andShchekin! for calculating Tolman’s length from the first-

4069

order curvature correction to the profile of the planner interface. Since nonmonotonic behavior of surface tension or
negative Tolman’s length is a common result of all recent
mean-field theories with different approximations, it is very
unlikely that our results are due to the random phase approximation to the perturbation term in the functional we
employed or to the potential functions ~LJ and Yukawa! forthe attractive part we employed. It is, however, a future subject to examine these results using more sophisticated approximations to the reference hard sphere part such as a nonlocal density functional approximation.28
We also discussed to what extent other thermodynamic
approaches10,11 can be used as alternatives to the rigorous
thermodynamic approach for studying the size dependenceof surface tension and of the length d. Although these approximate methods become formally exact in the large droplet limit, correct asymptotic behavior of the surface tension
and length d cannot be obtained via these routes due to theunderlying numerical-error problem in evaluating d.
The results summarized above suggest that in the application of computer simulation to the investigation of the size
dependence of surface tension of a liquid droplet, one should
consider the so-called mechanical routes9,16 or statistical mechanical routes23,29 rather than the thermodynamic routes.
The reasons are the following. The strict thermodynamic
route to the surface tension requires certain quantities which
are not readily calculated from computer simulations. Other
alternative thermodynamic routes proposed for computer
simulation are not valid for small droplets due to the approximations such as the use of Tolman’s equation, and canbe problematic for large droplets due to the numerical-error
problem. On the other hand, we note that through the statistical mechanical route as proposed by Blokhuis andBedeaux23 computer simulation can be still useful for the
calculation of d ` . This route demands information on thepair distribution functions of the planar liquid–vapor interface, which can be calculated straightforwardly via molecular simulation. Indeed results of molecular dynamics simulations of Haye and Bruin29 do not suffer from the large
numerical errors problem as appeared in the simulations
based on thermodynamic routes. However, their values of d `at several temperatures are all positive ~e.g., d `* 50.16 at
T * 50.75!, which contrast with our calculation based on the
combined theory that d ` is negative near the triple point( d ` 520.15s LJ at T50.70e LJ /k B for the LJ fluid; d `520.26d at T50.6T c for the Yukawa fluid!. As mentioned
above
other
mean-fieldordensity–functional
approaches5,6,25,26 also give rise to negative values of d ` .
Further investigation is needed to resolve the difference inprediction of the sign from the simulation and from the theoretical approaches.
The density–functional theory of nucleation when combined with the strict thermodynamic route to the surface tension is shown to be a very useful tool for examining the size
dependence of surface tension and of length d of a critical
droplet. Despite the fact that the absolute value of the surface
tension is dependent on the underlying approximations ~suchas the local density approximation!, this combined theory has
the virtue of self-consistency in calculation over a whole

4070

J. Chem. Phys., Vol. 109, No. 10, 8 September 1998

range of droplet sizes, including the planar surface limit.
Consequently, the density–functional theory provides aunique means for determining Tolman’s length. Nevertheless, the elusiveness of its sign still demands further studies.
ACKNOWLEDGMENTS

The authors are grateful to Dr. D. J. Diestler for helpful
discussions. XCZ thanks the NSF and ONR for support ofthis work. AKS also thanks RFBR ~Grant No. 96-02-18959!
for support of this work.
J. W. Gibbs, The Scientific Papers of J. Willard Gibbs ~Ox Bow, Woodbridge, NJ, 1993!.
2
R. C. Tolman, J. Chem. Phys. 17, 333 ~1949!.
3
A. H. Falls, L. E. Scriven, and H. T. Davis, J. Chem. Phys. 75, 3986
~1981!.
4
R. Guermeur, F. Biquard, and C. Jacolin, J. Chem. Phys. 82, 2040 ~1985!.
5
V. G. Baidakov and G. Sh. Boltachev, Russ. J. Phys. Chem. 69, 468~1995!.
6
M. Iwamatsu, J. Phys.: Condens. Matter 6, L173 ~1994!.
7
V. Talanquer and D. W. Oxtoby, J. Chem. Phys. 99, 2865 ~1995!.
8
A. Laaksonen and R. McGraw, Europhys. Lett. 35, 367 ~1996!.
9
A. I. Rusanov and E. N. Brodskaya, J. Colloid Interface Sci. 62, 542~1977!.
10
S. M. Thompson, K. E. Gubbins, J. P. R. B. Walton, R. A. R. Chantry, andJ. S. Rowlinson, J. Chem. Phys. 81, 530 ~1984!.
11
M. J. P. Nijmeijer, C. Bruin, A. B. van Woerkom, A. F. Bakker, and J. M.
J. van Leeuwen, J. Chem. Phys. 96, 565 ~1992!.
12
E. N. Brodskaya, J. C. Eriksson, A. Laaksonen, and A. I. Rusanov, J.
Colloid Interface Sci. 180, 86 ~1996!.
1

Koga, Zeng, and Shchekin
In fact, Tolman correctly noted in the discussion of the validity of Eq. ~1!
that d can be expected to vary appreciably when droplets are smallenough. Also Tolman predicted the rapid decrease in surface tension; theconclusion was, however, derived from Eq. ~1!.
14
D. W. Oxtoby and R. Evans, J. Chem. Phys. 89, 7521 ~1988!.
15
X. C. Zeng and D. W. Oxtoby, J. Chem. Phys. 94, 4472 ~1991!.
16
J. S. Rowlinson and B. Widom, Molecular Theory of Capillarity ~Oxford,
New York, 1989!.
17
P. Tarazona and R. Evans, Mol. Phys. 48, 799 ~1983!.
18
N. F. Carnahan and K. E. Starling, J. Chem. Phys. 51, 635 ~1969!.
19
J. D. Weeks, D. Chandler, and H. C. Anderson, J. Chem. Phys. 54, 5237
~1971!.
20
B. Q. Lu, R. Evans, and M. M. Telo da Gama, Mol. Phys. 55, 1913
~1985!.
21
R. Evans, Adv. Phys. 28, 143 ~1979!.
22
J. S. Rowlinson, J. Phys.: Condens. Matter 6, A1 ~1994!.
23
E. M. Blokhuis and D. Bedeaux, J. Chem. Phys. 97, 3576 ~1992!.
24
Strictly speaking, Thompson et al. used the truncated form of Tolman’sequation in their calculation; we used the original form of Tolman’s equation. However, we confirmed that our conclusion does not depend on the
choice.
25
A. E. van Giessen, E. M. Blokhuis, and D. J. Bukman, J. Chem. Phys.
108, 1148 ~1998!.
26
V. G. Bykov and A. K. Shchekin, in Proceedings of International Workshop, ‘‘Nucleation and Non-linear Problems in the First-Order Phase
Transition,’’ St. Petersburg, 1998.
27
E. M. Blokhuis and D. Bedeaux, Physica A 184, 42 ~1992!.
28
W. E. Brown, R. J. F. Decarvalho, and R. Evans, Mol. Phys. 88, 579~1996!.
29
M. J. Haye and C. Bruin, J. Chem. Phys. 100, 556 ~1994!.
13