arXiv:0704.0157v1  [cond-mat.stat-mech]  2 Apr 2007
7 Alternative Approaches to the Equilibrium
Properties of Hard-Sphere Liquids
M. López de Haro1, S. B. Yuste2 and A. Santos3
1 Centro de Investigación en Enerǵıa, Universidad Nacional Autónoma de México
(U.N.A.M.), Temixco, Morelos 62580, Mexico
malopez@servidor.unam.mx
2 Departamento de F́ısica, Universidad de Extremadura, E-06071 Badajoz, Spain
santos@unex.es
3 Departamento de F́ısica, Universidad de Extremadura, E-06071 Badajoz, Spain
andres@unex.es
An overview of some analytical approaches to the computation of the struc-
tural and thermodynamic properties of single component and multicompo-
nent hard-sphere fluids is provided. For the structural properties, they yield
a thermodynamically consistent formulation, thus improving and extending
the known analytical results of the Percus–Yevick theory. Approximate ex-
pressions for the contact values of the radial distribution functions and the
corresponding analytical equations of state are also discussed. Extensions of
this methodology to related systems, such as sticky hard spheres and square-
well fluids, as well as its use in connection with the perturbation theory of
fluids are briefly addressed.
1 Introduction
In the statistical thermodynamic approach to the theory of simple liquids,
there is a close connection between the thermodynamic and structural prop-
erties [1–4]. These properties depend on the intermolecular potential of the
system, which is generally assumed to be well represented by pair interactions.
The simplest model pair potential is that of a hard-core fluid (rods, disks,
spheres, hyperspheres) in which attractive forces are completely neglected. In
fact, it is a model that has been most studied and has rendered some analytical
results, although up to this day no general (exact) explicit expression for the
equation of state is available, except for the one-dimensional case. Something
similar applies to the structural properties. An interesting feature concerning
the thermodynamic properties is that in hard-core systems the equation of
state depends only on the contact values of the radial distribution functions.
In the absence of a completely analytical approach, the most popular methods
http://arxiv.org/abs/0704.0157v1
2 M. López de Haro, S. B. Yuste and A. Santos
to deal with both kinds of properties of these systems are integral equation
theories and computer simulations.
It is well known that in real gases and liquids at high temperatures the
state and thermodynamic properties are determined almost entirely by the
repulsive forces among molecules. At lower temperatures, attractive forces
become significant, but even in this case they affect very little the configu-
ration of the system at moderate and high densities. These facts are taken
into account in the application of the perturbation theory of fluids, where
hard-core fluids are used as the reference systems in the computation of the
thermodynamic and structural properties of real fluids. However, successful
results using perturbation theory are rather limited due to the fact that, as
mentioned above, there are in general no exact (analytical) expressions for
the thermodynamic and structural properties of the reference systems which
are in principle required in the calculations. On the other hand, in the realm
of soft condensed matter the use of the hard-sphere model in connection, for
instance, with sterically stabilized colloidal systems is quite common. This is
due to the fact that nowadays it is possible to prepare (almost) monodisperse
spherical colloidal particles with short-ranged harshly repulsive interparticle
forces that may be well described theoretically with the hard-sphere potential.
This chapter presents an overview of the efforts we have made over the
last few years to compute the thermodynamic and structural properties of
hard-core systems using relatively simple (approximate) analytical methods.
It is structured as follows. In Section 2 we describe our proposals to derive
the contact values of the radial distribution functions of a multicomponent
mixture (with an arbitrary size distribution, either discrete or continuous) of
d-dimensional hard spheres from the use of some consistency conditions and
the knowledge of the contact value of the radial distribution function of the
corresponding single component system. In turn, these contact values lead to
equations of state both for additive and non-additive hard spheres. Some con-
sequences of such equations of state, in particular the demixing transition, are
briefly analyzed. This is followed in Section 3 by the description of the Ratio-
nal Function Approximation method to obtain analytical expressions for the
structural quantities of three-dimensional single component and multicompo-
nent fluids. The only required inputs in this approach are the contact values of
the radial distribution functions and so the connection with the work of the
previous section follows naturally. Structural properties of related systems,
like sticky hard spheres or square-well fluids, that may also be tackled with
the same philosophy are also discussed in Section 4. Section 5 provides an
account of the reformulation of the perturbation theory of liquids using the
results of the Rational Function Approximation method for a single compo-
nent hard-sphere fluid and its illustration in the case of the Lennard–Jones
fluid. In the final section, we provide some perspectives of the achievements
obtained so far and of the challenges that remain ahead.
Alternative Approaches to Hard-Sphere Liquids 3
2 Contact Values and Equations of State for Mixtures
As stated in the Introduction, a nice feature of hard-core fluids is that the
expressions of all their thermodynamic properties in terms of the radial distri-
bution functions (RDF) are particularly simple. In fact, for these systems the
internal energy reduces to that of the ideal gas and in the pressure equation
it is only the contact values rather than the full RDF which appear explicitly.
In this section we present our approach to the derivation of the contact values
of hard-core fluid mixtures in d dimensions.
2.1 Additive Systems in d Dimensions
If σij denotes the distance of separation at contact between the centers of
two interacting fluid particles, one of species i and the other of species j,
the mixture is said to be additive if σij is just the arithmetic mean of the
hard-core diameters of each species. Otherwise, the system is non-additive.
We deal in this subsection and in Subsection 2.2 with additive systems, while
non-additive hard-core mixtures will be treated in Subsection 2.3.
Definitions
Let us consider an additive mixture of hard spheres (HS) in d dimensions with
an arbitrary number N of components. In fact, our discussion will remain valid
for N → ∞, i.e., for polydisperse mixtures with a continuous distribution of
sizes.
The additive hard core of the interaction between a sphere of species i and
a sphere of species j is σij =
(σi + σj), where the diameter of a sphere of
species i is σii = σi. Let the number density of the mixture be ρ and the mole
fraction of species i be xi = ρi/ρ, where ρi is the number density of species i.
From these quantities one can define the packing fraction η = vdρMd, where
vd = (π/4)
d/2/Γ (1 + d/2) is the volume of a d-dimensional sphere of unit
diameter and
Mn ≡ 〈σn〉 =
i (1)
denotes the nth moment of the diameter distribution.
In a HS mixture, the knowledge of the contact values gij(σij) of the RDF
gij(r), where r is the distance, is important for a number of reasons. For
example, the availability of gij(σij) is sufficient to get the equation of state
(EOS) of the mixture via the virial expression
Z(η) = 1 +
i,j=1
xixjσ
ijgij(σij), (2)
4 M. López de Haro, S. B. Yuste and A. Santos
where Z = p/ρkBT is the compressibility factor of the mixture, p being the
pressure, kB the Boltzmann constant, and T the absolute temperature.
The exact form of gij(σij) as functions of the packing fraction η, the set
of diameters {σk}, and the set of mole fractions {xk} is only known in the
one-dimensional case, where one simply has [5]
gij(σij) =
, (d = 1). (3)
Consequently, for d ≥ 2 one has to resort to approximate theories or empirical
expressions. For hard-disk mixtures, an accurate expression is provided by
Jenkins and Mancini’s (JM) approximation [6, 7],
gJMij (σij) =
(1− η)2
σiσjM1
σijM2
, (d = 2). (4)
The associated compressibility factor is
ZJM(η) =
1 + η/8
(1 − η)2
, (d = 2). (5)
In the case of three-dimensional systems, some important analytical expres-
sions for the contact values and the corresponding compressibility factor also
exist. For instance, the expressions which follow from the solution of the
Percus–Yevick (PY) equation of additive HS mixtures by Lebowitz [8] are
gPYij (σij) =
(1− η)2
σiσjM2
σijM3
, (d = 3), (6)
ZPY(η) =
(1− η)2
(1 − η)2
, (d = 3). (7)
Also analytical are the results obtained from the Scaled Particle Theory (SPT)
[9–12],
gSPTij (σij) =
(1 − η)2
σiσjM2
σijM3
(1 − η)3
σiσjM2
σijM3
, (d = 3),
ZSPT(η) =
(1− η)2
(1− η)3
, (d = 3). (9)
Neither the PY nor the SPT lead to particularly accurate values and so
Boubĺık [13] and, independently, Grundke and Henderson [14] and Lee and
Levesque [15] proposed an interpolation between the PY and the SPT contact
values, that we will refer to as the BGHLL values:
gBGHLLij (σij) =
(1− η)2
σiσjM2
σijM3
(1− η)3
σiσjM2
σijM3
, (d = 3).
Alternative Approaches to Hard-Sphere Liquids 5
This leads through Eq. (2) to the widely used and rather accurate Boubĺık–
Mansoori–Carnahan–Starling–Leland (BMCSL) EOS [13,16] for HS mixtures:
ZBMCSL(η) =
(1− η)2
η2(3− η)
(1 − η)3
, (d = 3). (11)
Refinements of the BGHLL values have been subsequently introduced, among
others, by Henderson et al. [17–22], Matyushov and Ladanyi [23], and Barrio
and Solana [24] to eliminate some drawbacks of the BMCSL EOS in the so-
called colloidal limit of binary HS mixtures. On a different path, but also
having to do with the colloidal limit, Viduna and Smith [25] have proposed
a method to obtain contact values of the RDF of HS mixtures from a given
EOS. However, none of these proposals may be easily generalized so as to be
valid for any dimensionality and any number of components. Therefore, if one
wants to have a more general framework able to deal with arbitrary d and N
an alternative strategy is called for.
Universality Ansatz
In order to follow our alternative strategy, it is useful to make use of exact
limit results that can help one in the construction of approximate expressions
for gij(σij). Let us consider first the limit in which one of the species, say i,
is made of point particles, i.e., σi → 0. In that case, gii(σi) takes the ideal
gas value, except that one has to take into account that the available volume
fraction is 1− η. Thus,
gii(σi) =
. (12)
An even simpler situation occurs when all the species have the same size,
{σk} → σ, so that the system becomes equivalent to a single component
system. Therefore,
{σk}→σ
gij(σij) = gs, (13)
where gs is the contact value of the RDF of the single component fluid at
the same packing fraction η as that of the mixture. Table 1 lists some of
the most widely used proposals for the contact value gs and the associated
compressibility factor
Zs = 1 + 2
d−1ηgs (14)
in the case of the single component HS fluid.
Equations (12) and (13) represent the simplest and most basic conditions
that gij(σij) must satisfy. There is a number of other less trivial consistency
conditions [11, 17, 19, 20, 23, 24,32–34], some of which will be used later on.
In order to proceed, in line with a property shared by earlier proposals
[see, in particular, Eqs. (4), (6), (8), and (10)], we assume that, at a given
packing fraction η, the dependence of gij(σij) on the parameters {σk} and
{xk} takes place only through the scaled quantity
6 M. López de Haro, S. B. Yuste and A. Santos
Table 1. Some expressions of gs and Zs for the single component HS fluid. In the
SHY proposal, ηcp = (
3/6)π is the crystalline close-packing fraction for hard disks.
In the LM proposal, b3 and b4 are the (reduced) third and fourth virial coefficients,
ζ(η) = 1.2973(59)−0.062(13)η/ηcp for d = 4, and ζ(η) = 1.074(16)+0.163(45)η/ηcp
for d = 5, where the values of the close-packing fractions are ηcp = π
2/16 ≃ 0.617
and ηcp = π
2/30 ≃ 0.465 for d = 4 and d = 5, respectively.
d gs Zs Label Ref.
1− 7η/16
(1− η)2
1 + η2/8
(1− η)2
H [26]
1− η(2ηcp − 1)/2η2cp
1− 2η + η2(ηcp − 1)/2η2cp
1− 2η + η2(ηcp − 1)/2η2cp
SHY [27]
2 gHs −
27(1− η)4
ZHs −
26(1− η)4
L [28]
1 + η/2
(1− η)2
1 + 2η + 3η2
(1− η)2
PY [29]
1− η/2 + η2/4
(1− η)3
1 + η + η2
(1− η)3
SPT [9]
1− η/2
(1− η)3
1 + η + η2 − η3
(1− η)3
CS [30]
1 + [21−db3 − ζ(η)b4/b3]η
1− ζ(η)(b4/b3)η + [ζ(η)− 1] 21−db4η2
1 + 2d−1ηgLMs LM [31]
zij ≡
. (15)
More specifically, we assume
gij(σij) = G(η, zij), (16)
where the function G(η, z) is universal in the sense that it is a common func-
tion for all the pairs (i, j), regardless of the composition and number of compo-
nents of the mixture. Of course, the function G(η, z) is in principle different for
each dimensionality d. To clarify the implications of this universality ansatz,
let us imagine two mixtures M and M′ having the same packing fraction
η but strongly differing in the set of mole fractions, the sizes of the parti-
cles, and even the number of components. Suppose now that there exists a
pair (i, j) in mixture M and another pair (i′, j′) in mixture M′ such that
zij = zi′j′ . Then, according to Eq. (16), the contact value of the RDF for the
pair (i, j) in mixture M is the same as that for the pair (i′, j′) in mixture
M′, i.e., gij(σij) = gi′j′ (σi′j′ ). In order to ascribe a physical meaning to the
parameter zij , note that the ratio Md−1/Md can be understood as a “typ-
ical” inverse diameter (or curvature) of the particles of the mixture. Thus,
z−1ij =
(σ−1i + σ
j )/(Md−1/Md) represents the arithmetic mean curvature,
in units of Md−1/Md, of a particle of species i and a particle of species j.
Alternative Approaches to Hard-Sphere Liquids 7
Once the ansatz (16) is adopted, one may use the limits in (12) and (13)
to get G(η, z) at z = 0 and z = 1, respectively. Since zii → 0 in the limit
σi → 0, insertion of Eq. (12) into (16) yields
G(η, 0) = 1
≡ G0(η). (17)
Next, if all the diameters are equal, zij → 1, so that Eq. (13) implies that
G(η, 1) = gs. (18)
Linear Approximation
As the simplest approximation [35], one may assume a linear dependence of
G on z that satisfies the basic requirements (17) and (18), namely
G(η, z) = 1
z. (19)
Inserting this into Eq. (16), one has
ge1ij (σij) =
. (20)
Here, the label “e1” is meant to indicate that (i) the contact values used are
an extension of the single component contact value gs and that (ii) G(η, z) is
a linear polynomial in z. This notation will become handy below. Although
the proposal (20) is rather crude and does not produce especially accurate
results for gij(σij) when d ≥ 3, it nevertheless leads to an EOS that exhibits
an excellent agreement with simulations in 2, 3, 4, and 5 dimensions, provided
that an accurate gs is used as input [35–39]. This EOS may be written as
Ze1(η) = 1 +
2d−1(Ω0 −Ω1) + [Zs(η)− 1]Ω1, (21)
where the coefficients Ωm depend only on the composition of the mixture and
are defined by
Ωm = 2
−(d−m)
Mmd−1
Mm+1d
Mn+mMd−n. (22)
In particular, for d = 2 and d = 3,
Ze1(η) =
Zs(η) −
, (d = 2), (23)
8 M. López de Haro, S. B. Yuste and A. Santos
Ze1(η) =
Zs(η)−
, (d = 3). (24)
As an extra asset, from Eq. (21) one may write the virial coefficients of
the mixture Bn, defined by
Z = 1 +
Bn+1ρ
n, (25)
in terms of the (reduced) virial coefficients of the single component fluid bn
defined by
Zs = 1 +
bn+1η
n. (26)
The result is
Bn = v
Ω1bn + 2
d−1(Ω0 −Ω1)
. (27)
In the case of binary mixtures, these coefficients are in very good agreement
with the available exact and simulation results [35,37], except when the mix-
ture involves components of very disparate sizes, especially for high dimen-
sionalities. One may perform a slight modification such that this deficiency is
avoided and thus get a modified EOS [37, 40]. For d = 2 and d = 3 it reads
Z(η) = Zs(η) + x1
1− η2
1− η2
− Zs(η)
σ2 − σ1
1− η1
1− η1
− Zs(η)
σ1 − σ2
, (d = 2, 3),
where ηi = vdρiσ
i is the partial volume packing fraction due to species i. In
contrast to most of the approaches (PY, SPT, BMCSL, e1, . . . ), the proposal
(28) expresses Z(η) in terms not only of Zs(η) but also involves Zs
and Zs
. Equation (28) should in principle be useful in particular for
binary mixtures involving components of very disparate sizes. However, it is
slightly less accurate than the one given in Eq. (21) for ordinary mixtures [37].
Quadratic Approximation
In order to improve the proposal contained in Eq. (20), in addition to the
consistency requirements (12) and (13), one may consider the condition stem-
ming from a binary mixture in which one of the species (say i = 1) is much
Alternative Approaches to Hard-Sphere Liquids 9
larger than the other one (i.e., σ1/σ2 → ∞), but occupies a negligible volume
(i.e., x1(σ1/σ2)
d → 0). In that case, a sphere of species 1 is felt as a wall by
particles of species 2, so that [17, 20, 41]
σ1/σ2→∞
x1(σ1/σ2)d→0
g12(σ12)− 2d−1ηg22(σ2)
= 1. (29)
Hence, in the limit considered in Eq. (29), we have z22 → 1, z12 → 2. Conse-
quently, under the universality ansatz (16), one may rewrite Eq. (29) as
G(η, 2) = 1 + 2d−1ηG(η, 1). (30)
Thus, Eqs. (17), (18), and (30) provide complete information on the function
G at z = 0, z = 1, and z = 2, respectively, in terms of the contact value gs of
the single component RDF.
The simplest functional form of G that complies with the above consistency
conditions is a quadratic function of z [42]:
G(η, z) = G0(η) + G1(η)z + G2(η)z2, (31)
where the coefficients G1(η) and G2(η) are explicitly given by
G1(η) = (2− 2d−2η)gs −
2− η/2
1− η , (32)
G2(η) =
1− η/2
− (1− 2d−2η)gs. (33)
Therefore, the explicit expression for the contact values is
ge2ij (σij) =
(2− 2d−2η)gs −
2− η/2
1− η/2
− (1− 2d−2η)gs
. (34)
Following the same criterion as the one used in connection with Eq. (20), the
label “e2” is meant to indicate that (i) the resulting contact values represent
an extension of the single component contact value gs and that (ii) G(η, z)
is a quadratic polynomial in z. Of course, the quadratic form (31) is not the
only choice compatible with conditions (17), (18), and (30). For instance, a
rational function was also considered in Ref. [42]. However, although it is
rather accurate, it does not lead to a closed form for the EOS. In contrast,
when Eq. (34) is inserted into Eq. (2), one gets a closed expression for the
compressibility factor in terms of the packing fraction η and the first few
moments Mn, n ≤ d. The result is
Ze2(η) = 1 + 2
d−2 η
1− η [2(Ω0 − 2Ω1 +Ω2) + (Ω1 −Ω2)η]
+ [Zs(η)− 1]
2Ω1 −Ω2 + 2d−2(Ω2 −Ω1)η
, (35)
10 M. López de Haro, S. B. Yuste and A. Santos
where the quantities Ωm are defined in Eq. (22). Quite interestingly, in the
two-dimensional case Eq. (35) reduces to Eq. (23), i.e.,
Ze1(η) = Ze2(η), (d = 2). (36)
This illustrates the fact that two different proposals for the contact values
gij(σij) can yield the same EOS when inserted into Eq. (2). On the other
hand, for three-dimensional mixtures Eq. (35) becomes
Ze2(η) =
1− η + M
Zs(η)−
, (d = 3),
which differs from Eq. (24). In fact,
Ze1(η)− Ze2(η) =
1 + η
− (1 − 2η)Zs(η)
, (d = 3).
Specific Examples
In this subsection, rather than carrying out an exhaustive comparison with
the wealth of results available in the literature, we will consider only a few
representative examples. In particular, for d = 3, we will restrict ourselves
to a comparison with classical proposals (say BGHLL, PY, and SPT for the
contact values). The comparison with more recent ones may be found in Refs.
[35, 42, 43].
Thus far the development has been rather general since gs remains free in
Eqs. (20) and (34). In order to get specific results, it is necessary to fix gs [cf.
Table 1]. In the one-dimensional case, one has gs = 1/(1− η) and so one gets
the exact result (3) after substitution into Eq. (20). Similarly Eqs. (32) and
(33) lead to G1 = G2 = 0 and so we recover again the exact result.
If in the two-dimensional case we take Henderson’s value [26] gs = g
then the linear approximation (20) reduces to the JM approximation, Eq. (4).
This equivalence can be symbolically represented as geH1ij = g
ij , where the
label “eH1” refers to the extension of Henderson’s single component value in
the linear approximation. While gJMij is very accurate, even better results are
provided by the quadratic form (34), especially if Luding’s value [28] gs = g
is used [44].
In the three-dimensional case, Eq. (20) is of the form of the solution of the
PY equation [8]. In fact, insertion of gs = g
s leads to Eq. (6), i.e., g
gPYij . Similarly, if the SPT expression [9] gs = g
s is used for the single
component contact value in the quadratic approximation (34), we reobtain
the SPT expression for the mixture, Eq. (8). In other words, geSPT2ij = g
On the other hand, if the much more accurate CS [30] expression gs = g
used as input, we arrive at the following expression:
Alternative Approaches to Hard-Sphere Liquids 11
geCS2ij =
η(1 − η/3)
(1 − η)2
σiσjM2
σijM3
η2(1− η/2)
(1− η)3
σiσjM2
σijM3
, (d = 3),
which is different from the BGHLL one, Eq. (10), improves the latter for
zij > 1, and leads to similar results for zij < 1, as comparison with computer
simulations shows [42]. The four approximations (6), (8), (10), and (39) are
consistent with conditions (12) and (13), but only the SPT and eCS2 are also
consistent with condition (29). It should also be noted that if one considers
a binary mixture in the infinite solute dilution limit, namely x1 → 0, so that
z12 → 2/(1 + σ2/σ1), Eq. (39) yields the same result for g12(σ12) as the one
proposed by Matyushov and Ladanyi [23] for this quantity on the basis of
exact geometrical relations. However, the extension that the same authors
propose when there is a non-vanishing solute concentration, i.e., for x1 6= 0,
is different from Eq. (39).
Equation (34) can also be used in the case of hyperspheres (d ≥ 4) [42]. In
particular, a very good agreement with available computer simulations [38] is
obtained for d = 4 and d = 5 by using Luban and Michels [31] value gs = g
0.0 0.1 0.2 0.3 0.4 0.5
-0.02
 eCS1
 eCS2
Fig. 1. Deviation of the compressibility factor from the BMCSL value, as a function
of the packing fraction η for an equimolar three-dimensional binary mixture with
σ2/σ1 = 0.6. The open (Ref. [18]) and closed (Ref. [45]) circles are simulation data.
The lines are the PY EOS (– · · –), the SPT EOS (– · – ·), the eCS1 EOS (· · · ), and
the eCS2 EOS (– – –).
12 M. López de Haro, S. B. Yuste and A. Santos
0.00 0.05 0.10 0.15 0.20 0.25
 4D, / =1/2
 4D, / =1/3
 5D, / =2/5
Fig. 2. Compressibility factor for three equimolar mixtures in 4D and 5D systems.
Lines are the eLM1 predictions, while symbols are simulation data [38].
Now we turn to the compressibility factors (21) and (35), which are ob-
tained from the contact values (20) and (34), respectively. Since they depend
on the details of the composition through the d first moments, they are mean-
ingful even for continuous polydisperse mixtures.
As said above, in the two-dimensional case both Eqs. (21) and (35) reduce
to Eq. (23), which yield very accurate results when a good Zs is used as
input [39, 42, 44]. For three-dimensional mixtures, insertion of Zs = Z
Eqs. (24) and (37) yields
ZeCS1(η) = ZBMCSL(η) +
(1 − η)3M23
M1M3 −M22
, (d = 3), (40)
ZeCS2(η) = ZBMCSL(η)−
(1 − η)2M23
M1M3 −M22
, (d = 3), (41)
where ZBMCSL(η) is given by Eq. (11). Note that ZeCS1(η) > ZBMCSL(η) >
ZeCS2(η). Since simulation data indicate that the BMCSL EOS tends to un-
derestimate the compressibility factor, it turns out that, as illustrated in Fig.
1 for an equimolar binary mixture with σ2/σ1 = 0.6, the performance of ZeCS1
is, paradoxically, better than that of ZeCS2 [42], despite the fact that the un-
derlying linear approximation for the contact values is much less accurate than
the quadratic approximation. This shows that a rather crude approximation
such as Eq. (20) may lead to an extremely good EOS [35, 37–39], which, as
Alternative Approaches to Hard-Sphere Liquids 13
clearly seen in Fig. 1, represents a substantial improvement over the classical
proposals. Interestingly, the EOS corresponding to ZeCS1 has recently been
independently derived as the second order approximation of the Fundamental
Measure Theory for the HS fluid by Hansen-Goos and Roth [46].
In the case of d = 4 and d = 5, use of Zs(η) = Z
s (η) in Eq. (21)
produces a simple extended EOS of a mixture of hard additive hyperspheres
in these dimensionalities. The accuracy of these two EOS for hard hypersphere
mixtures in the fluid region has been confirmed by simulation data [38] for a
wide range of compositions and size ratios. In Fig. 2, this accuracy is explicitly
exhibited in the case of three equimolar mixtures, two in 4D and one in 5D.
2.2 A More Consistent Approximation for Three-Dimensional
Additive Mixtures
Up to this point, we have considered an arbitrary dimensionality d and have
constructed, under the universality assumption (16), the acurate quadratic
approximation (34), which fulfills the consistency conditions (12), (13), and
(29). However, there exist extra consistency conditions that are not necessarily
satisfied by (34). In particular, when the mixture is in contact with a hard wall,
the state of equilibrium imposes that the pressure evaluated near the wall by
considering the impacts with the wall must be the same as the pressure in the
bulk evaluated from the particle-particle collisions. This consistency condition
is especially important if one is interested in deriving accurate expressions for
the contact values of the particle-wall correlation functions.
Since a hard wall can be seen as a sphere of infinite diameter, the contact
value gwj of the correlation function of a sphere of diameter σj with the wall
can be obtained from gij(σij) as
gwj = lim
gij(σij). (42)
Note that gwj provides the ratio between the density of particles of species
j adjacent to the wall and the density of those particles far away from the
wall. The sum rule connecting the pressure of the fluid and the above contact
values is [47]
Zw(η) =
xjgwj, (43)
where the subscript w in Zw has been used to emphasize that Eq. (43) repre-
sents a route alternative to the virial one, Eq. (2), to get the EOS of the HS
mixture. The condition Z = Zw is equivalent to (29) in the special case where
one has a single fluid in the presence of the wall. However, in the general case
of a mixture plus a wall, the condition Z = Zw is stronger than Eq. (29).
In the two-dimensional case, it turns out that the quadratic approximation
(34) already satisfies the requirement Z = Zw, regardless of the density and
composition of the mixture [44]. However, this is not the case for d ≥ 3.
14 M. López de Haro, S. B. Yuste and A. Santos
Our problem now consists of computing gij(σij) and the associated gwj for
the HS mixture in the presence of a hard wall, so that the condition Z = Zw
is satisfied for an arbitrary mixture [43]. Due to the mathematical complexity
of the problem, here we will restrict ourselves to three-dimensional systems
(d = 3). Similarly to what we did in the preceding subsection, we consider
a class of approximations of the universal type (16), so that conditions (12)
and (13) lead again to Eqs. (17) and (18), respectively. Notice that Eq. (16)
implies in particular that
gwj = G(η, zwj), zwj = 2σj
. (44)
Assuming that z = 0 is a regular point and taking into account condition
(17), G(η, z) can be expanded in a power series in z:
G(η, z) = G0(η) +
Gn(η)zn. (45)
After simple algebra, using the ansatz (16) and Eq. (45) in Eqs. (2) (with
d = 3) and (43) one gets
Z = G0 + 3η
G0 + 4η
Mn+13
i,j=1
xixjσ
ij , (46)
Zw = G0 +
Mn. (47)
Notice that if the series (45) is truncated after a given order n ≥ 3, Zw is given
by the first n moments of the size distribution only. On the other hand, Z still
involves an infinite number of moments if the truncation is made after n ≥ 4
due to the presence of terms like
i,j xixjσ
4/σij ,
i,j xixjσ
5/σ2ij , . . . .
Therefore, if we want the consistency condition Z = Zw to be satisfied for any
discrete or continuous polydisperse mixture, either the whole infinite series
(45) needs to be considered or it must be truncated after n = 3. The latter is
of course the simplest possibility and thus we make the approximation
G(η, z) = G0(η) + G1(η)z + G2(η)z2 + G3(η)z3. (48)
As a consequence, Z and Zw depend functionally on the size distribution of
the mixture only through the first three moments (which is in the spirit of
Rosenfeld’s Fundamental Measure Theory [48]).
Using the approximation (48) in Eqs. (46) and (47) we are led to
Z = G0 + η
(3G0 + 2G1) + 2
(G1 + 2G2 + 2G3)
, (49)
Alternative Approaches to Hard-Sphere Liquids 15
Zw = G0 + 2
G1 + 4
(G2 + 2G3) . (50)
Thus far, the dependence of both Z and Zw on the momentsM1,M2, andM3
is explicit and we only lack the packing-fraction dependence of G1, G2, and
G3. From Eqs. (49) and (50) it follows that the difference between Z and Zw
is given by
Z−Zw =
[3ηG0 − 2(1− η)G1]+2
[ηG1 − 2(1− η)G2 − 2(2− η)G3] .
Therefore, Z = Zw for any dispersity provided that
G1(η) =
2 (1− η)2
, (52)
G2(η) =
4 (1− η)3
− 2− η
G3(η), (53)
where use has been made of the definition of G0, Eq. (17). To close the problem,
we use the equal size limit given in Eq. (18), which yields G0+G1+G2+G3 = gs.
After a little algebra we are led to
G2(η) = (2− η)gs −
2 + η2/4
(1− η)2
, (54)
G3(η) = (1− η)
gSPTs − gs
. (55)
This completes the derivation of our improved approximation, which we will
call “e3”, following the same criterion as the one used to call “e1” and “e2” to
the approximations (20) and (34), respectively. In Eq. (55), gSPTs is the SPT
contact value for a single fluid, whose expression appears in Table 1. From
Eq. (55) it is obvious that the choice gs = g
s makes our e3 approximation
to become the e2 approximation, both reducing to the SPT for mixtures, Eq.
(8). This means that the SPT is fully internally consistent with the require-
ment Z = Zw, although it has the shortcoming of not being too accurate in
the single component case. The e3 proposal, on the other hand, satisfies the
condition Z = Zw and has the flexibility of accommodating any desired gs.
For the sake of concreteness, let us write explicitly the contact values in
the e3 aproximation:
ge3ij (σij) =
2 (1− η)2
(2 − η)gs −
2 + η2/4
(1− η)2
+ (1− η)
gSPTs − gs
, (56)
16 M. López de Haro, S. B. Yuste and A. Santos
ge3wj =
1− η +
(1− η)2
σj + 4
(2 − η)gs −
2 + η2/4
(1− η)2
+8(1− η)
gSPTs − gs
. (57)
With the above results the compressibility factor may be finally written in
terms of Zs as
Ze3(η) =
(1− η)
(1− η)2
Zs(η)−
. (58)
A few comments are in order at this stage. First, from Eq. (49) we can
observe that, for the class of approximations (48), the compressibility factor
Z does not depend on the individual values of the coefficients G2 and G3,
but only on their sum. As a consequence, two different approximations of
the form (48) sharing the same density dependence of G1 and G2 + G3 also
share the same virial EOS. For instance, if one makes the choice gs = g
then ZePY3 = ZPY, even though g
ij (σij) 6= gPYij (σij). Furthermore, if one
makes the more accurate choice gs = g
s , then ZeCS3 = ZBMCSL, but again
geCS3ij (σij) 6= gBGHLLij (σij). The eCS3 contact values are
geCS3ij (σij) =
2 (1− η)2
η2(1 + η)
4(1− η)3
4(1− η)2
, (59)
geCS3wj =
(1− η)2
η2(1 + η)
(1− η)3
(1 − η)2
. (60)
In Figs. 3 and 4 we display the performance of the contact values as given
by Eqs. (59) and (60), respectively, by comparison with results of computer
simulations for both discrete and polydisperse mixtures. In both figures we
have also included the results that follow from the classical proposals as well
as those of the eCS1 and eCS2 approximations. It is clear that for the wall-
particle contact values the eCS3 approximation yields the best performance,
while for the particle-particle contact values both the eCS2 and eCS3 are of
comparable accuracy. A further feature to be pointed out is that the practical
collapse on a common curve of the simulation data in Figs. 3 and 4 provide a
posteriori support for the universality ansatz made in Eq. (16).
As mentioned earlier, there exist extra consistency conditions (see for in-
stance Ref. [12]) that one might use as well within our approach. Assuming
Alternative Approaches to Hard-Sphere Liquids 17
0.0 0.5 1.0 1.5 2.0
 eCS1
 eCS2 
 eCS3 
Fig. 3. Plot of the difference gij(σij)− gBGHLLij (σij) as a function of the parameter
zij = (σiσj/σij)M2/M3 for hard spheres (d = 3) at a packing fraction η = 0.49.
The symbols are simulation data for the single fluid (circle, Ref. [36]), three binary
mixtures (squares, Ref. [49]) with σ2/σ1 = 0.3 and x1 = 0.0625, 0.125, and 0.25, and
a ternary mixture (triangles, Ref. [50]) with σ2/σ1 =
, σ3/σ1 =
, and x1 = 0.1,
x2 = 0.2. The lines are the PY approximation (– · · –), the SPT approximation
(– · – ·), the eCS1 approximation (· · · ), the eCS2 approximation (– – –), and the
eCS3 approximation (—).
that the ansatz (16) still holds, some of these conditions are related to the
derivatives of G with respect to z, namely
∂G(η, z)
2(1− η)2
, (61)
∂2G(η, z)
gPYs −
, (62)
∂3G(η, z)
= 0. (63)
Interestingly enough, as shown by Eq. (52), condition (61) is already satisfied
by our e3 approximation without having to be imposed. On the other hand,
condition (63) implies G3 = 0 in the e3 scheme and thus it is only satisfied if
gs = g
s , in which case we recover the SPT. Condition (62) is not fulfilled
either by the SPT or by the e3 approximation (except for a particular ex-
pression of gs which is otherwise not very accurate). Thus, fulfilling the extra
18 M. López de Haro, S. B. Yuste and A. Santos
0.0 0.5 1.0 1.5 2.0 2.5
 eCS1
 eCS2
 eCS3
Fig. 4. Plot of the difference gwj − gBGHLLwj as a function of the parameter zwj/2 =
σjM2/M3 for hard spheres (d = 3) at a packing fraction η = 0.4. The symbols are
simulation data for a polydisperse mixture with a narrow top-hat distribution (open
squares, Ref. [51]), a polydisperse mixture with a wide top-hat distribution (open
circles, Ref. [51]), a polydisperse mixture with a Schulz distribution (open triangles,
Ref. [51]), and a binary mixture (closed circles, Ref. [52]). The lines are the PY
approximation (– · · –), the SPT approximation (– · – ·), the eCS1 approximation
(· · · ), the eCS2 approximation (– – –), and the eCS3 approximation (—).
conditions (62) and (63) with a free gs requires either considering a higher
order polynomial in z (in which case the consistency condition Z = Zw can-
not be satisfied for arbitrary mixtures, as discussed before) or not using the
universality ansatz at all. In the first case, we have checked that a quartic or
even a quintic polynomial does not improve matters, whereas giving up the
universality assumption increases significantly the number of parameters to
be determined and seems not to be adequate in view of the behavior observed
in the simulation data.
An additional comment has to do with the restriction to d = 3 in this
subsection. As noted before, the approximation e1 reduces to the exact result
(3) for d = 1. For d = 2, the approximation e2 already fulfills the condition
Z = Zw and so there is no real need to go further in that case. Since we
have needed the approximation e3 to satisfy Z = Zw for d = 3, it is tempting
to speculate that a polynomial form for G(z) of degree d could be found to
be consistent with the condition Z = Zw for d ≥ 4. However, a detailed
analysis shows that this is not the case for an arbitrary mixture, since the
Alternative Approaches to Hard-Sphere Liquids 19
number of conditions exceeds the number of unknowns, unless the universality
assumption is partially relaxed.
As a final comment, let us stress that, although the discussion in this
section has referred, for the sake of simplicity, to discrete mixtures, all the
dependence on the details of the composition occurs through a finite number
of moments, so that the results remain meaningful even for continuous poly-
disperse mixtures [53]. In that case, instead of a set of mole fractions {xi} and
a set of diameters {σi}, one has to deal with a distribution function w(σ) such
that w(σ)dσ is the fraction of particles with a diameter comprised between σ
and σ + dσ. Therefore, the moments (1) are now defined as
dσ σnw(σ), (64)
and with such a change the results we have derived for discrete mixtures also
hold for polydisperse systems.
2.3 Non-Additive Systems
Non-additive hard-core mixtures, where the distance of closest approach be-
tween particles of different species is no longer the arithmetic mean of the
diameters of both particles, have received much less attention than additive
mixtures, in spite of their in principle more versatility to deal with interesting
aspects occurring in real systems (such as fluid-fluid phase separation) and
of their potential use as reference systems in perturbation calculations on the
thermodynamic and structural properties of, say, Lennard–Jones mixtures.
Nevertheless, the study of non-additive systems goes back fifty years [54–56]
and is still a rapidly developing and challenging problem.
As mentioned in the paper by Ballone et al. [57], where the relevant
references may be found, experimental work on alloys, aqueous electrolyte
solutions, and molten salts suggests that hetero-coordination and homo-
coordination may be interpreted in terms of excluded volume effects due to
non-additivity of the repulsive part of the intermolecular potential. In particu-
lar, positive non-additivity leads naturally to demixing in HS mixtures, so that
some of the experimental findings of phase separation in the above mentioned
(real) systems may be accounted for by using a model of a binary mixture of
(positive) non-additive HS. On the other hand, negative non-additivity seems
to account well for chemical short-range order in amorphous and liquid binary
mixtures with preferred hetero-coordination [58].
Some Preliminary Definitions
Let us consider an N -component mixture of non-additive HS in d dimensions.
In this case, σij =
(σi + σj)(1 + ∆ij), where ∆ij ≥ −1 is a symmetric
matrix with zero diagonal elements (∆ii = 0) that characterizes the degree
20 M. López de Haro, S. B. Yuste and A. Santos
of non-additivity of the interactions. If ∆ij > 0 the non-additivity character
of the ij interaction is said to be positive, while it is negative if ∆ij < 0. In
the case of a binary mixture (N = 2), the only non-additivity parameter is
∆ ≡ ∆12 = ∆21. The virial EOS (2) remains being valid in the non-additive
case.
The contact values gij(σij) can be expanded in a power series in density
gij(σij) = 1 + vdρ
xkck;ij + (vdρ)
k,ℓ=1
xkxℓckℓ;ij +O(ρ3). (65)
The coefficients ck;ij , ckℓ;ij , . . . are independent of the composition of the mix-
ture, but they are in general complicated nonlinear functions of the diameters
σij , σik, σjk, σkℓ, . . . . Insertion of the expansion (65) into Eq. (2) yields the
virial expansion of Z, namely
Z(ρ) = 1 +
Bn(vdρ)
= 1 + vdρ
i,j=1
Bijxixj + (vdρ)
i,j,k=1
Bijkxixjxk
+(vdρ)
i,j,k,ℓ=1
Bijkℓxixjxkxℓ +O(ρ4). (66)
Note that, for further convenience, we have introduced the coefficients Bn ≡
−(n−1)
d Bn, where Bn are the usual virial coefficients [cf. Eq. (25)]. The
composition-independent second, third, and fourth (barred) virial coefficients
are given by
Bij = 2
d−1σdij , (67)
Bijk =
ck;ijσ
ij + cj;ikσ
ik + ci;jkσ
, (68)
Bijkℓ =
ckℓ;ijσ
ij + cjℓ;ikσ
ik + ciℓ;jkσ
jk + cjk,iℓσ
iℓ + cik,jℓσ
+cij;kℓσ
. (69)
A Simple Proposal for the Equation of State of d-Dimensional
Non-Additive Mixtures
Our goal now is to generalize the e1 proposal given by Eq. (20) to the non-
additive case [59]. We will not try to extend the e2 and e3 proposals, Eqs. (34)
and (56), because of two reasons. First, given the inherent complexity of non-
additive systems, we want to keep the approach as simple as possible. Second,
Alternative Approaches to Hard-Sphere Liquids 21
we are more interested in the EOS than in the contact values themselves and,
as mentioned earlier, the e1 proposal provides excellent EOS, at least in the
additive case, despite the simplicity of the corresponding contact values.
As the simplest possible extension, we impose again the point particle and
equal size consistency conditions, Eqs. (12) and (13), and thus keep in this
case also the ansatz (16) and the linear structure of Eq. (19). However, instead
of using Eq. (15), we determine the parameters zij as to reproduce Eq. (65)
to first order in the density. The result is readily found to be [59]
zij =
)−1(∑
k xkck;ij
. (70)
Here b2 = 2
d−1 and b3 are the second and third virial coefficients for the
single component fluid, as defined by Eq. (26). The proposal of Eq. (19) sup-
plemented by Eq. (70) is, by construction, accurate for densities low enough
as to justify the truncated approximation gij(σij) ≈ 1 + vdρ
k xkck;ij . On
the other hand, the limitations of this truncated expansion for moderate and
large densities may be compensated by the use of gs. When Eqs. (16), (19),
and (70) are inserted into Eq. (2) one gets
Z(η) = 1 +
b3MdB2 − b2B3
(b3 − b2)M2d
+ [Zs(η) − 1]
B3 −MdB2
(b3 − b2)M2d
. (71)
Equation (71) is the sought generalization of Eq. (21) to non-additive hard-
core systems. As in the additive case, the the density dependence in the EOS
of the mixture is rather simple: Z(η)− 1 is expressed as a linear combination
of η/(1 − η) and Zs(η) − 1, with coefficients such that the second and third
virial coefficients are reproduced. Again, Eq. (71) is bound to be accurate for
sufficiently low densities, while the limitations of the truncated expansion for
moderate and large densities are compensated by the use of the EOS of the
pure fluid.
The exact second virial coefficient B2 is known from Eq. (67). In principle,
one should use the exact coefficients ck;ij to compute B3. However, to the
best of our knowledge they are only known for d ≤ 3. Since our objective is to
have a proposal which is explicit for any d, we can make use of a reasonable
approximation for them [59], as described below.
An Approximate Proposal for ck;ij
The values of the coefficients ck;ij are exactly known for d = 1 and d = 3 and
from these results one may approximate them in d dimensions as [59]
ck;ij = σ
k;ij +
σd−1k;ij
σi;jkσj;ik, (72)
where we have called
22 M. López de Haro, S. B. Yuste and A. Santos
σk;ij ≡ σik + σjk − σij (73)
and it is understood that σk;ij ≥ 0 for all sets ijk. Clearly, σi;ij = σi. For a
binary mixture Eq. (72) yields
c1;11 = (b3/b2)σ
c2;11 = (2σ12 − σ1)d + (b3/b2 − 1)σ1(2σ12 − σ1)d−1,
c1;12 = σ
1 + (b3/b2 − 1) (2σ12 − σ1)σd1/σ12.
Of course, Eqs. (72) and (74) reduce to the exact results for d = 1 (b2 = b3 = 1)
and for d = 3 (b2 = 4, b3 = 10).
The quantities σk;ij may be given a simple geometrical interpretation.
Assume that we have three spheres of species i, j, and k aligned in the sequence
ikj. In such a case, the distance of closest approach between the centers of
spheres i and j is σik + σjk. If the sphere of species k were not there, that
distance would of course be σij . Therefore σk;ij as given by Eq. (73) represents
a kind of effective diameter of sphere k, as seen from the point of view of the
interaction between spheres i and j.
Inserting Eq. (72) into Eq. (70), one gets
zij =
)−1(∑
k xkσ
k xkσ
k;ij σi;jkσj;ik
Mdσij
. (75)
It can be easily checked that in the additive case (σk;ij → σk), Eq. (75) reduces
to Eq. (15).
Equations (72) and (74) are restricted to the situation σk;ij ≥ 0 for any
choice of i, j, and k, i.e., 2σ12 ≥ max(σ1, σ2) in the binary case. This excludes
the possibility of dealing with mixtures with extremely high negative non-
additivity in which one sphere of species k might “fit in” between two spheres
of species i and j in contact. Since for d = 3 and N = 2 the coefficients ck;ij
are also known for such mixtures [60], we may extend our proposal to deal
with these cases:
c1;11 = (b3/b2)σ
c2;11 = σ̂
2 + (b3/b2 − 1)σ1σ̂d−12 ,
c1;12 = (2σ12 − σ̂2)d + (b3/b2 − 1) σ̂2σd1/σ12,
where we have defined
σ̂2 = max (2σ12 − σ1, 0) . (77)
With such an extension, we recover the exact values of ck;ij for a binary
mixture of hard spheres (d = 3), even if σ1 > 2σ12 or σ2 > 2σ12.
The EOS (71) becomes explicit when B3 is obtained from Eq. (68) by
using the approximation (72). The resulting virial coefficient is the exact one
for d = 1 and d = 3. For hard disks (d = 2), it turns out that the approximate
third virial coefficient is practically indistinguishable from the exact one [59].
Alternative Approaches to Hard-Sphere Liquids 23
-0.6 -0.4 -0.2 0.0 0.2 0.4 0.6
Fig. 5. Plot of the compressibility factor versus the non-additivity parameter ∆ for
a symmetric binary mixture of non-additive hard spheres (d = 3) at η = π/30 and
two different compositions. The solid lines are our proposal, Eq. (71), with Zs = Z
while the dashed lines are Hamad’s proposal (Refs. [61–63]). The symbols are results
from Monte Carlo simulations (Refs. [64,65]).
When the approximate B3 is used, Eq. (71) reduces to Eq. (21) in the additive
case.
From the comparison with simulation results, both for the compressibility
factor and higher order virial coefficients, we find that the EOS (71) does a
good job for non-additive mixtures, thus representing a reasonable compro-
mise between simplicity and accuracy, provided that Zs is accurate enough.
This is illustrated in Fig. 5, where the proposal (71) with Zs = Z
s and a
similar proposal by Hamad [61–63] are compared with simulation data [64,65]
for some three-dimensional symmetric mixtures. A more extensive compari-
son [59] shows that Eq. (71) seems to work better (especially as the density
is increased) in the case of positive non-additivities, at least for d = 1, d = 2,
and d = 3, but its performance is also reasonably good in highly asymmetric
mixtures, even for negative ∆. Of course the full assessment of this proposal
is still pending since it involves many facets (non-additivity parameters, size
ratios, density, and composition). Without this full assessment and given its
rather satisfactory performance so far, going beyond the approximation given
by Eq. (19) (taking similar steps to the ones described in Subsections 2.1 and
2.2 for additive systems) does not seem to be necessary at this stage, although
it is in principle feasible.
24 M. López de Haro, S. B. Yuste and A. Santos
2.4 Demixing
Demixing is a common phase transition in fluid mixtures usually originated
on the asymmetry of the interactions (e.g., their strength and/or range) be-
tween the different components in the mixture. In the case of athermal systems
such as HS mixtures in d dimensions, if fluid-fluid separation occurs, it would
represent a neat example of an entropy-driven phase transition, i.e., a phase
separation based only on the size asymmetry of the components. The exis-
tence of demixing in binary additive three dimensional HS mixtures has been
studied theoretically since decades, and the issue is still controversial. In this
subsection we will present our results following different but related routes
that attempt to clarify some aspects of this problem.
Binary Mixtures of Additive d-Dimensional Spheres (d = 3, d = 4
and d = 5)
Now we look at the possible instability of a binary fluid mixture of HS of
diameters σ1 and σ2 (σ1 > σ2) in d dimensions by looking at the Helmholtz
free energy per unit volume, f , which is given by
= −1 +
xi ln
Z(η′)− 1
, (78)
where λi is the thermal de Broglie wavelength of species i. We locate the
spinodals through the condition f11f22−f212 = 0, with fij ≡ ∂2f/∂ρi∂ρj . Due
to the spinodal instability, the mixture separates into two phases of different
composition. The coexistence conditions are determined through the equality
of the pressure p and the two chemical potentials µ1 and µ2 in both phases
(µi = ∂f/∂ρi), leading to binodal (or coexistence) curves.
We begin with the case d = 3. It is well known that the BMCSL EOS, Eq.
(11), does not lead to demixing. However, other EOS for HS mixtures have
been shown to predict demixing [41, 66], including the EOS that is obtained
by truncating the virial series after a certain number of terms [67, 68]. In
particular, it turns out that both Z = ZeCS1, Eq. (40), and Z = ZeCS2,
Eq. (41), lead to demixing for certain values of the parameter γ ≡ σ2/σ1
that measures the size asymmetry. The critical values of the pressure, the
composition, and the packing fraction are presented in Table 2 for a few
values of γ.
As discussed earlier, the eCS1 EOS and, to a lesser extent, the eCS2 EOS
are both in reasonably good agreement with the available simulation results
for the compressibility factor [18, 36, 45] and lead to the exact second and
third virial coefficients but differ in the predictions for Bn with n ≥ 4. The
scatter in the values for the critical constants shown in Table 2 is evident and
so there is no indication as to whether one should prefer one equation over
the other in connection with this problem. Notice, for instance, that the eCS2
Alternative Approaches to Hard-Sphere Liquids 25
Table 2. Critical constants pcσ
1/kBT , x1c, and ηc for different γ-values as obtained
from the two extended CS equations (40) and (41).
eCS1 eCS2
γ pcσ
1/kBT x1c ηc pcσ
1/kBT x1c ηc
0.05 3599 0.0093 0.822 1096 0.0004 0.204
0.1 1307 0.0203 0.757 832.0 0.0008 0.290
0.2 653.4 0.0537 0.725 — — —
0.3 581.9 0.0998 0.738 — — —
0.4 663.4 0.1532 0.766 — — —
does not predict demixing for γ ≥ 0.2, while both the values of the critical
pressures and packing fractions for which it occurs according to the eCS1 EOS
suggest that the transition might be metastable with respect to a fluid-solid
transition.
Now we turn to the cases d = 4 and d = 5. Here we use the extended
Luban–Michels equation (eLM1) described in Subsection 2.1 [see Eq. (21)
and Table 1]. As seen in Fig. 6, the location of the critical point tends to go
down and to the right in the η2 vs η1 plane as γ decreases for d = 4 [69].
On the other hand, while it also tends to go down as γ decreases if d = 5,
its behavior in the η2 vs η1 plane is rather more erratic in this case. Also,
the value of the critical pressure pc (in units of kBT/σ
1) is not a monotonic
function of γ; its minimum value lies between γ = 1/3 and γ = 1/2 when
d = 4, and it is around γ = 3/5 for d = 5. This non-monotonic behavior is
also observed for three-dimensional HS [66, 68].
It is conceivable that the demixing transition in binary mixtures of hard
hyperspheres in four and five dimensions described above may be metastable
with respect to a fluid-solid transition, as it may also be the case of 3D HS.
In fact, the value of the pressure at the freezing transition for the single
component fluid is [31] pfσ
d/kBT ≃ 12.7 (d = 3), 11.5 (d = 4), and 12.2
(d = 5), i.e., pfσ
d/kBT does not change appreciably with the dimensionality
but is clearly very small in comparison with the critical pressures pcσ
1/kBT
we obtain for the mixture; for instance, pcσ
1/kBT ≃ 600 (d = 3, γ = 3/10),
300 (d = 4, γ = 1/3) and 123 (d = 5, γ = 3/5). However, one should also
bear in mind that, if the concentration x1 of the bigger spheres decreases, the
value of the pressure at which the solid-fluid transition in the mixture occurs
in 3D is also considerably increased with respect to pf [cf. Fig. 6 of Ref. [66]].
Thus, for concentrations x1 ≃ 0.01 corresponding to the critical point of the
fluid-fluid transition, the maximum pressure of the fluid phase greatly exceeds
pf. If a similar trend with composition also holds in 4D and 5D, and given that
the critical pressures become smaller as the dimensionality d is increased, it is
not clear whether the competition between the fluid-solid and the fluid-fluid
transitions in these dimensionalities will always be won by the former. The
point clearly deserves further investigation.
26 M. López de Haro, S. B. Yuste and A. Santos
0.0 0.1 0.2 0.3 0.4 0.5 0.6
1400  2/ 1=3/4
4 1/k
0.0 0.1 0.2 0.3 0.4
5 1/k
Fig. 6. Spinodal curves (upper panels: lines) and binodal curves (upper panels: open
symbols; lower panels: lines) in a 4D system (left panels) and in a 5D system (right
panels). The closed symbols are the critical consolute points.
An interesting feature must be mentioned. There is a remarkable similarity
between the binodal curves represented in the pσdi –η1 and in the µi–η1 planes
[69]. By eliminating η1 as if it were a parameter, one can represent the binodal
curves in a µi vs pσ
i plane. Provided the origin of the chemical potentials
is such as to make λi = σi, the binodals in the µi–pσ
i plane practically
collapse into a single curve (which is in fact almost a straight line) for each
dimensionality (d = 3, d = 4, and d = 5) [69]. A closer analysis of this
Alternative Approaches to Hard-Sphere Liquids 27
phenomenon shows, however, that it is mainly due to the influence on µi
of terms which are quantitatively dominant but otherwise irrelevant to the
coexistence conditions.
Binary Mixtures of Non-Additive Hard Hyperspheres in the Limit
of High Dimensionality
Let us now consider a binary mixture of non-additive HS of diameters σ1
and σ2 in d dimensions. Thus in this case σ12 ≡ 12 (σ1 + σ2)(1 +∆) where as
before ∆ may be either positive or negative. Further assume (something that
will become exact in the limit d → ∞ [70]) that the EOS of the mixture is
described by the second virial coefficient only, namely
p = ρkBT [1 +B2(x1)ρ] , (79)
where, according to Eq. (67),
B2(x1) = vd2
1 + x
2 + 2x1x2σ
. (80)
The Helmholtz free energy per unit volume is given by f/ρkBT = −1 +∑2
i=1 xi ln
+B2ρ, where Eq. (78) has been used. The Gibbs free energy
per particle is
g = (f + p)/ρ =
xi ln
+ 2B2(x1)ρ, (81)
where without loss of generality we have set kBT = 1. Given a size ratio γ, a
value of ∆, and a dimensionality d, the consolute critical point (x1c, pc) is the
solution to
∂2g/∂x21
∂3g/∂x31
= 0, provided of course it exists. Then,
one can get the critical density ρc from Eq. (79).
We now introduce the scaled quantities [71]
p̃ ≡ 2d−1vdd−2pσd1/kBT, u ≡ d−1B2ρ. (82)
Consequently, Eqs. (79) and (81) can be rewritten as
p̃ = u
u+ d−1
/B̃2, (83)
xi ln (xiΛi) + ln
Adu/B̃2
+ 2du, (84)
where B̃2 ≡ B2/2d−1vdσd1 , Λi ≡ (λi/σ1)d, and Ad ≡ d/2d−1vd. Next we take
the limit d → ∞ and assume that the volume ratio γ̃ ≡ γd is kept fixed and
that there is a (slight) non-additivity ∆ = d−2∆̃ such that the scaled non-
additivity parameter ∆̃ is also kept fixed in this limit. Thus, the second virial
coefficient can be approximated by
28 M. López de Haro, S. B. Yuste and A. Santos
B̃2 = B̃
2 +B̃
−1+O(d−2), B̃(0)2 =
x1 + x2γ̃
2 = x1x2γ̃
1/2J,
J ≡ 1
(ln γ̃)
+ 2∆̃. (86)
Let us remark that, in order to find a consolute critical point, it is essential
to keep the term of order d−1 if ∆̃ ≤ 0. The EOS (83) can then be inverted
to yield
u = u(0)+u(1)d−1+O(d−2), u(0) =
2 , u
(1) = −1
1− u(0) B̃
In turn, the Gibbs free energy (84) becomes
g = g(0)d+ g(1) +O(d−1),
g(0) = 2u(0), g(1) =
i=1 xi ln (xiΛi) + ln
(0)/B̃
+ 2u(1),
while the chemical potentials µ1 = g+x2 (∂g/∂x1)p and µ2 = g−x1 (∂g/∂x1)p
are given by
µi = µ
i d+ µ
i +O(d−1), µ
1 = 2p̃
1 = ln
Adx1Λ1
p̃/B̃
2 + (x2/x1)(γ̃p̃)
1/2B̃
2 /B̃
where µ2 is obtained from µ1 by the changes x1 ↔ x2, Λ1 → Λ2/γ̃, γ̃ → 1/γ̃,
p̃→ p̃γ̃, B̃2 → B̃2/γ̃.
The coordinates of the critical point are readily found to be
x1c =
γ̃3/4
1 + γ̃3/4
, p̃c =
1 + γ̃1/4
4γ̃J2
. (90)
Note that x1c is independent of ∆̃. The coexistence curve, which has to be
obtained numerically, follows from the conditions µ
i (xA, p̃) = µ
i (xB , p̃)
(i = 1, 2) where x1 = xA and x1 = xB are the mole fractions of the co-
existing phases. Once the critical consolute point has been identified in the
pressure/concentration plane, we can obtain the critical density. The domi-
nant behaviors of B̃2 and u at the critical point are
2 (x1c) =
1− γ̃1/4 + γ̃1/2
)2 , u
1 + γ̃1/4
1− γ̃1/4 + γ̃1/2
. (91)
Hence, the critical density readily follows after substitution in the scaling
relation given in Eq. (82). It is also convenient to consider the scaled version
η̃ ≡ d−12dη of the packing fraction η = vdρσd1 (x1 + x2γ̃). At the critical point,
it takes the nice expression
Alternative Approaches to Hard-Sphere Liquids 29
0.00 0.04 0.08 0.12 0.16
p  =-0.1
 =0.1
=0.01
Fig. 7. Binodal curves in the planes ep vs x1 and eη vs x1 corresponding to eγ = 0.01
and e∆ = −0.1, e∆ = 0, and e∆ = 0.1.
η̃c =
γ̃1/8 + γ̃−1/8
. (92)
The previous results clearly indicate that a demixing transition is possible
not only for additive or positively non-additive mixtures but even for negative
non-additivities. The only requirement is J > 0, i.e., ∆̃ > − 1
(ln γ̃)
equivalently, ∆ > − 1
(ln γ)
. Figure 7 shows the binodal curves corresponding
to γ̃ = 0.01 and ∆̃ = −0.1 (negative non-additivity), ∆̃ = 0 (additivity), and
∆̃ = 0.1 (positive non-additivity).
While the high dimensionality limit has allowed us to address the prob-
lem in a mathematically simple and clear-cut way, the possibility of demixing
with negative non-additivity is not an artifact of that limit. As said before,
demixing is known to occur for positive non-additive binary mixtures of HS
in three dimensions and there is compelling evidence on the existence of this
phenomenon in the additive case, at least in the metastable fluid region. Even
though in a three-dimensional mixture the EOS is certainly more complicated
than Eq. (79) and the demixing transition that we have just discussed for neg-
ative non-additivity is possibly metastable with respect to the freezing transi-
tion, the main effects at work (namely the competition between depletion due
to size asymmetry and hetero-coordination due to negative non-additivity)
are also present. In fact, it is interesting to point out that Roth et al. [72],
using the approximation of an effective single component fluid with pair inter-
30 M. López de Haro, S. B. Yuste and A. Santos
actions to describe a binary mixture of non-additive 3D HS and employing an
empirical rule based on the effective second virial coefficient, have also sug-
gested that demixing is possible for small negative non-additivity and high
size asymmetry. Our exact results lend support to this suggestion and con-
firm that, in some cases, the limit d→ ∞ highlights features already present
in real systems.
3 The Rational Function Approximation (RFA) Method
for the Structure of Hard-Sphere Fluids
The RDF g(r) and its close relative the (static) structure factor S(q) are the
basic quantities used to discuss the structure of a single component fluid [1–4].
The latter quantity is defined as
S(q) = 1 + ρh̃(q), (93)
where
h̃(q) =
dr e−iq·rh(r) (94)
is the Fourier transform of the total correlation function h(r) ≡ g(r)−1, i being
the imaginary unit. An important related quantity is the direct correlation
function c(r), which is defined in Fourier space through the Ornstein–Zernike
(OZ) relation [1–4]
c̃(q) =
h̃(q)
1 + ρh̃(q)
, (95)
where c̃(q) is the Fourier transform of c(r)
The usual approach to obtain g(r) is through one of the integral equa-
tion theories, where the OZ equation is complemented by a closure relation
between c(r) and h(r) [1]. However, apart from requiring in general hard nu-
merical labor, a disappointing aspect is that the substitution of the (necessar-
ily) approximate values of g(r) obtained from them in the (exact) statistical
mechanical formulae may lead to the thermodynamic inconsistency problem.
The two basic routes to obtain the EOS of a single component fluid of HS
are the virial route, Eq. (14), and the compressibility route
χs ≡ kBT
= [1− ρc̃(0)]−1 = S(0)
= 1 + 2ddησ−d
dr rd−1h(r). (96)
Thermodynamic consistency implies that
χ−1s (η) =
[ηZs(η)], (97)
Alternative Approaches to Hard-Sphere Liquids 31
but, in general, this condition is not satisfied by an approximate RDF. In
the case of a HS mixture, the virial route is given by Eq. (2), while the
compressibility route is indicated below [cf. Eq. (145)].
In this section we describe the RFA method, which is an alternative to
the integral equation approach and in particular leads by construction to
thermodynamic consistency.
3.1 The Single Component HS Fluid
We begin with the case of a single component fluid of HS of diameter σ.
The following presentation is equivalent to the one given in Refs. [73, 74],
where all details can be found, but more suitable than the former for direct
generalization to the case of mixtures.
The starting point will be the Laplace transform
G(s) =
dr e−srrg(r) (98)
and the auxiliary function Ψ(s) defined through
G(s) =
[ρ+ esσΨ(s)]
. (99)
The choice of G(s) as the Laplace transform of rg(r) and the definition of
Ψ(s) from Eq. (99) are suggested by the exact form of g(r) to first order in
density [73].
Since g(r) = 0 for r < σ while g(σ+) = finite, one has
g(r) = Θ(r − σ)
g(σ+) + g′(σ+)(r − σ) + · · ·
, (100)
where g′(r) ≡ dg(r)/dr. This property imposes a constraint on the large s
behavior of G(s), namely
eσssG(s) = σg(σ+) +
g(σ+) + σg′(σ+)
s−1 +O(s−2). (101)
Therefore, lims→∞ e
sσsG(s) = σg(σ+) = finite or, equivalently,
s−2Ψ(s) =
2πσg(σ+)
= finite. (102)
On the other hand, according to Eq. (96) with d = 3,
χs = 1− 24ησ−3 lim
dr e−srr [g(r) − 1]
= 1− 24ησ−3 lim
G(s)− s−2
. (103)
Since the (reduced) isothermal compressibility χs is also finite, one has∫∞
dr r2 [g(r) − 1] = finite, so that the weaker condition
dr r [g(r)− 1] =
lims→0[G(s)− s−2] = finite must hold. This in turn implies
32 M. López de Haro, S. B. Yuste and A. Santos
Ψ(s) = −ρ+ ρσs−
ρσ2s2+
ρσ3 +
ρσ3 +
σs4+O(s5).
(104)
First-Order Approximation (PY Solution)
An interesting aspect to be remarked is that the minimal input we have just
described on the physical requirements related to the structure and thermo-
dynamics of the system is enough to determine the small and large s limits of
Ψ(s), Eqs. (102) and (104), respectively. While infinite choices for Ψ(s) would
comply with such limits, a particularly simple form is a rational function. In
particular, the rational function having the least number of coefficients to be
determined is
Ψ(s) =
E(0) + E(1)s+ E(2)s2 + E(3)s3
L(0) + L(1)s
, (105)
where one of the coefficients can be given an arbitrary non-zero value. We
choose E(3) = 1. With such a choice and in view of Eq. (104), one finds
E(0) = −ρL(0), E(1) = −ρ(L(1) − σL(0)), E(2) = ρ(σL(1) − 1
σ2L(0)), and
L(0) = 2π
1 + 2η
(1− η)2
, (106)
L(1) = 2πσ
1 + η/2
(1− η)2
. (107)
Upon substitution of these results into Eqs. (99) and (105), we get
G(s) =
L(0) + L(1)s
ϕ2(σs)σ3L(0) + ϕ1(σs)σ2L(1)
] , (108)
where
ϕn(x) ≡ x−(n+1)
(−x)m
− e−x
. (109)
In particular,
ϕ0(x) =
1− e−x
, ϕ1(x) =
1− x− e−x
, ϕ2(x) =
1− x+ x2/2− e−x
(110)
Note that limx→0 ϕn(x) = (−1)n/(n+ 1)!.
It is remarkable that Eq. (108), which has been derived here as the sim-
plest rational form for Ψ(s) complying with the requirements (102) and (104),
coincides with the solution to the PY closure, c(r) = 0 for r > σ, of the OZ
equation [29]. Application of Eq. (102) yields the PY contact value gPYs and
compressibility factor ZPYs shown in Table 1. Analogously, Eq. (103) yields
χPYs =
(1− η)4
(1 + 2η)2
. (111)
It can be easily checked that the thermodynamic relation (97) is not satisfied
by the PY theory.
Alternative Approaches to Hard-Sphere Liquids 33
Second-Order Approximation
In the spirit of the RFA, the simplest extension of the rational approximation
(105) involves two new terms, namely αs4 in the numerator and L(2)s2 in the
denominator, both of them necessary in order to satisfy Eq. (102). Such an
addition leads to
Ψ(s) =
E(0) + E(1)s+ E(2)s2 + E(3)s3 + αs4
L(0) + L(1)s+ L(2)s2
. (112)
Applying Eq. (104), it is possible to express E(0), E(1), E(2), E(3), L(0), and
L(1) in terms of α and L(2). This leads to
G(s) =
L(0) + L(1)s+ L(2)s2
1 + αs− ρ
ϕ2(σs)σ3L(0) + ϕ1(σs)σ2L(1) + ϕ0(σs)σL(2)
(113)
where
L(0) = 2π
1 + 2η
(1− η)2
, (114)
L(1) = 2πσ
1 + 1
(1 − η)2
1 + 2η
α− 3ηL
. (115)
Thus far, irrespective of the values of the coefficients L(2) and α, the condi-
tions lims→∞ e
sσsG(s) = finite and lims→0[G(s) − s−2] = finite are satisfied.
Of course, if L(2) = α = 0, one recovers the PY approximation. More gen-
erally, we may determine these coefficients by prescribing the compressibility
factor Zs (or equivalently the contact value gs) and then, in order to ensure
thermodynamic consistency, compute from it the isothermal compressibility
χs by means of Eq. (97). From Eqs. (102) and (103) one gets
L(2) = 2πασgs, (116)
1− 12η
1 + 2
αL(2)
. (117)
Clearly, upon substitution of Eqs. (114) and (116) into Eq. (117) a quadratic
algebraic equation for α is obtained. The physical root is
α = − 12η(1 + 2η)E4
(1− η)2 + 36η [1 + η − Zs(1− η)]E4
, (118)
where
Zs − 13
Zs − 13
Zs − ZPYs
)]1/2}
. (119)
The other root must be discarded because it corresponds to a negative value
of α, which, according to Eq. (116), yields a negative value of L(2). This would
34 M. López de Haro, S. B. Yuste and A. Santos
imply the existence of a positive real value of s at which G(s) = 0 [73, 74],
which is not compatible with a positive definite RDF. However, according to
the form of Eq. (119) it may well happen that, once Zs has been chosen, there
exists a certain packing fraction ηg above which α is no longer positive. This
may be interpreted as an indication that, at the packing fraction ηg where α
vanishes, the system ceases to be a fluid and a glass transition in the HS fluid
occurs [74–76].
Expanding (113) in powers of s and using Eq. (101) one can obtain the
derivatives of the RDF at r = σ+ [77]. In particular, the first derivative is
g′(σ+) =
L(1) − L(2)
, (120)
which may have some use in connection with perturbation theory [15].
It is worthwhile to point out that the structure implied by Eq. (113) coin-
cides in this single component case with the solution of the Generalized Mean
Spherical Approximation (GMSA) [78], where the OZ relation is solved under
the ansatz that the direct correlation function has a Yukawa form outside the
core.
For a given Zs, once G(s) has been determined, inverse Laplace trans-
formation yields rg(r). First, note that Eq. (99) can be formally rewritten
G(s) = −
ρn−1 [−Ψ(s)]−n e−nsσ. (121)
Thus, the RDF is then given by
g (r) =
ρn−1ψn (r − nσ)Θ (r − nσ) , (122)
with Θ (x) denoting the Heaviside step function and
ψn (r) = −L−1
s [−Ψ (s)]−n
, (123)
L−1 denoting the inverse Laplace transform. Explicitly, using the residue the-
orem,
ψn (r) = −
(n−m)!(m− 1)!
rn−m, (124)
where
a(i)mn = lim
s [−Ψ (s) /(s− si)]−n , (125)
si (i = 1, . . . , 4) being the poles of 1/Ψ(s), i.e., the roots of E
(0) + E(1)s +
E(2)s2 + E(3)s3 + αs4 = 0. Explicit expressions of g(r) up to the second
coordination shell σ ≤ r ≤ 3σ can be found in Ref. [79].
Alternative Approaches to Hard-Sphere Liquids 35
On the other hand, the static structure factor S(q) [cf. Eq. (93)] and the
Fourier transform h̃(q) may be related to G(s) by noting that
h̃(q) =
dr r sin(qr)h(r) = −2π G(s)−G(−s)
. (126)
Therefore, the basic structural quantities of the single component HS fluid,
namely the RDF and the static structure factor, may be analytically deter-
mined within the RFA method once the compressibility factor Zs, or equiva-
lently the contact value gs, is specified. In Fig. 8 we compare simulation data
of g(r) for a density ρσ3 = 0.9 [80] with the RFA prediction and a recent
approach by Trokhymchuk et al. [81], where Zs = Z
s [cf. Table 1] and the
associated compressibility
χCSs =
(1− η)4
1 + 4η + 4η2 − 4η3 + η4
(127)
are taken in both cases. Both theories are rather accurate, but the RFA cap-
tures better the maxima and minima of g(r) [82].
It is also possible to obtain within the RFA method the direct correlation
function c(r). Using Eqs. (95) and (126), and applying the residue theorem,
one gets, after some algebra,
1 2 3 4
1 2 3 4
Fig. 8. Radial distribution function of a single component HS fluid for ρσ3 = 0.9.
The solid lines represent simulation data [80]. The dashed lines represent the results
of the approach of Ref. [81], while the dotted lines refer to those of the RFA method.
The inset shows the oscillations of g(r) in more detail.
36 M. López de Haro, S. B. Yuste and A. Santos
c(r) =
+K0 +K1r +K3r
Θ(1 − r) +K
(128)
where
12αηL(2)/π + 1− 12α(1 + 2α)η/(1− η), (129)
4α2(1 − η)4κ6
2 [1 + 2(1 + 3α)η]± [2 + η + 2α(1 + 2η)]κ
+(1− η)
κ2 − η (12 + (κ± 6)κ)
L(2)/π
12η [1 + 2(1 + 3α)η]
±6η [3η − 2α(1− 4η)]κ− 6η(1 + 2α)(1− η)κ2 − (1 − η)2κ3(ακ∓ 1)
+6η(1− η)
κ2 − η (12 + (κ± 6)κ)
L(2)/π
, (130)
K−1 = −
κ +K−e
−κ +K0 +K1 +K3
, (131)
K0 = −
1 + 2 (1 + 3α) η − 6η (1− η)L(2)/π
ακ (1− η)2
, (132)
2α2κ2 (1− η)4
[2 + η + 2α(1 + 2η)]
2 − 4 (1− η) [1 + η
×(7 + η + 6α (2 + η))]L(2)/π + 12η (2 + η) (1− η)2L(2)
,(133)
K0, (134)
K = − (K+ +K− +K−1) . (135)
In Eqs. (129)–(135) we have taken σ = 1 as the length unit. Note that Eq.
(135) guarantees that c(0) = finite, while Eq. (131) yields c(σ+) − c(σ−) =
L(2)/2πα = g(σ+). The latter equation proves the continuity of the indirect
correlation function γ(r) ≡ h(r) − c(r) at r = σ. With the above results,
Eqs. (122) and (128), one may immediately write the function γ(r). Finally,
we note that the bridge function B(r) is linked to γ(r) and to the cavity
(or background) function y(r) ≡ eφ(r)/kBT g(r), where φ(r) is the interaction
potential, through
B(r) = ln y(r) − γ(r), (136)
and so, within the RFA method, the bridge function is also completely speci-
fied analytically for r > σ once Zs is prescribed.
If one wants to have B(r) also for 0 ≤ r ≤ σ, then an expression for the
cavity function is required in that region. Here we propose such an expression
using a limited number of constraints. First, since the cavity function and its
first derivative are continuous at r = σ, we have
Alternative Approaches to Hard-Sphere Liquids 37
0.0 0.2 0.4 0.6 0.8 1.0
1000 3=0.3
3=0.5
3=0.7
Fig. 9. Cavity function of a single component HS fluid in the overlap region for
ρσ3 = 0.3, 0.5, and 0.7. The solid lines represent our proposal (140) with Zs = Z
while the symbols represent Monte Carlo simulation results [84].
y(1) = gs,
y′(1)
− 1, (137)
where Eqs. (116) and (120) have been used and again σ = 1 has been taken.
Next, we consider the following exact zero-separation theorems [83]:
ln y(0) = Zs(η)− 1 +
′)− 1
, (138)
y′(0)
= −6ηy(1). (139)
The four conditions (137)–(139) can be enforced by assuming a cubic poly-
nomial form for ln y(r) inside the core, namely
y(r) = exp
Y0 + Y1r + Y2r
2 + Y3r
, (0 ≤ r ≤ 1), (140)
where
Y0 = Zs(η)− 1 +
′)− 1
, (141)
Y1 = −6ηy(1), (142)
38 M. López de Haro, S. B. Yuste and A. Santos
0 2 4 6 8 10
 RFA, =0.3
 RFA, =0.49
Fig. 10. Parametric plot of the bridge function B(r) versus the indirect correlation
function γ(r). The dashed line refers to the RFA for η = 0.3, while the solid line
refers to the RFA for η = 0.49. In each case, the branch of the curve to the right
of the circle corresponds to r ≤ 1, while that to the left corresponds to r ≥ 1. For
comparison, the PY closure B(r) = ln[1 + γ(r)]− γ(r) is also plotted (dash-dotted
line).
Y2 = 3 ln y(1)−
y′(1)
− 3Y0 − 2Y1, (143)
Y3 = −2 ln y(1) +
y′(1)
+ 2Y0 + Y1. (144)
The proposal (140) is compared with available Monte Carlo data [84] in Fig.
9, where an excellent agreement can be observed.
Once the cavity function y(r) provided by the RFA method is comple-
mented by (140), the bridge function B(r) can be obtained at any distance.
Figure 10 presents a parametric plot of the bridge function versus the indirect
correlation function as given by the RFA method for two different packing
fractions, as well as the result associated with the PY closure. The fact that
one gets a smooth curve means that within the RFA the oscillations in γ(r)
are highly correlated to those of B(r). Further, the effective closure relation in
the RFA turns out to be density dependent, in contrast with what occurs for
the PY theory. Note that the absolute value |B(r)| for a given value of γ(r)
is smaller in the RFA than the PY value and that the RFA and PY curves
become paradoxically closer for larger densities. Since the PY theory is known
Alternative Approaches to Hard-Sphere Liquids 39
to yield rather poor values of the cavity function inside the core [85, 86], it
seems likely that the present differences may represent yet another manifes-
tation of the superiority of the RFA method, a point that certainly deserves
to be further explored.
3.2 The Multicomponent HS Fluid
The method outlined in the preceding subsection will be now extended to an
N -component mixture of additive HS. Note that in a multicomponent system
the isothermal compressibility χ is given by
χ−1 =
T,{xj}
T,{xj}
= 1− ρ
i,j=1
xixj c̃ij(0), (145)
where c̃ij(q) is the Fourier transform of the direct correlation function cij(r),
which is defined by the OZ equation
h̃ij(q) = c̃ij(q) +
ρkh̃ik(q)c̃kj(q), (146)
where hij(r) ≡ gij(r) − 1. Equations (145) and (146) are the multicompo-
nent extensions of Eqs. (96) and (95), respectively. Introducing the quantities
ĥij(q) ≡
ρiρj h̃ij(q) and ĉij(q) ≡
ρiρj c̃ij(q), the OZ relation (146) be-
comes, in matrix notation,
ĉ(q) = ĥ(q) · [I+ ĥ(q)]−1, (147)
where I is the N ×N identity matrix. Thus, Eq. (145) can be rewritten as
χ−1 =
i,j=1
xixj [δij − ĉij(0)] =
i,j=1
I+ ĥ(0)
. (148)
Similarly to what we did in the single component case, we introduce the
Laplace transforms of rgij(r):
Gij(s) =
dr e−srrgij(r). (149)
The counterparts of Eqs. (100) and (101) are
gij(r) = Θ(r − σij)
gij(σ
ij) + g
ij)(r − σij) + · · ·
, (150)
40 M. López de Haro, S. B. Yuste and A. Santos
eσijssGij(s) = σijgij(σ
ij) +
gij(σ
ij) + σijg
s−1 +O(s−2). (151)
Moreover, the condition of a finite compressibility implies that h̃ij(0) = finite.
As a consequence, for small s,
s2Gij(s) = 1 +H
3 + · · · (152)
with H
ij = finite and H
ij = −h̃ij(0)/4π = finite, where
dr (−r)nrhij(r). (153)
We are now in the position to generalize the approximation (113) to the
N -component case [87]. While such a generalization may be approached in a
variety of ways, two motivations are apparent. On the one hand, we want to
recover the PY result [8] as a particular case in much the same fashion as in
the single component system. On the other hand, we want to maintain the
development as simple as possible. Taking all of this into account, we propose
Gij(s) =
e−σijs
L(s) · [(1 + αs)I− A(s)]−1
, (154)
where L(s) and A(s) are the matrices
Lij(s) = L
ij + L
ij s+ L
2, (155)
Aij(s) = ρi
ϕ2(σis)σ
ij + ϕ1(σis)σ
ij + ϕ0(σis)σiL
, (156)
the functions ϕn(x) being defined by Eq. (109). We note that, by construc-
tion, Eq. (154) complies with the requirement lims→∞ e
σijssGij(s) = finite.
Further, in view of Eq. (152), the coefficients of s0 and s in the power series
expansion of s2Gij(s) must be 1 and 0, respectively. This yields 2N
2 condi-
tions that allow us to express L(0) and L(1) in terms of L(2) and α. The solution
is [87]
ij = ϑ1 + ϑ2σj + 2ϑ2α− ϑ1
ρkσkL
kj , (157)
ij = ϑ1σij +
ϑ2σiσj + (ϑ1 + ϑ2σi)α−
ρkσkL
kj , (158)
where ϑ1 ≡ 2π/(1− η) and ϑ2 ≡ 6π(M2/M3)η/(1− η)2.
In parallel with the development of the single component case, L(2) and α
can be chosen arbitrarily. Again, the choice L
ij = α = 0 gives the PY solution
[8, 88]. Since we want to go beyond this approximation, we will determine
those coefficients by taking prescribed values for gij(σij), which in turn, via
Eq. (2), give the EOS of the mixture. This also leads to the required value of
Alternative Approaches to Hard-Sphere Liquids 41
χ−1 = ∂(ρZ)/∂ρ, thus making the theory thermodynamically consistent. In
particular, according to Eq. (151),
ij = 2πασijgij(σ
ij). (159)
The condition related to χ is more involved. Making use of Eq. (152), one can
get h̃ij(0) = −4πH(1)ij in terms of L(2) and α and then insert it into Eq. (148).
Finally, elimination of L
ij in favor of α from Eq. (159) produces an algebraic
equation of degree 2N , whose physical root is determined by the requirement
that Gij(s) is positive definite for positive real s. It turns out that the physical
solution corresponds to the smallest of the real roots. Once α is known, upon
substitution into Eqs. (154), (157), (158), and (159), the scheme is complete.
Also, using Eq. (151), one can easily derive the result
g′ij(σ
ij) =
2πασij
ij − L
. (160)
It is straightforward to check that the results of the preceding subsection are
recovered by setting σi = σ, regardless of the values of the mole factions.
Once Gij(s) has been determined, inverse Laplace transformation directly
yields rgij(r). Although in principle this can be done analytically, it is more
practical to use one of the efficient methods discussed by Abate and Whitt [89]
to numerically invert Laplace transforms [90].
In Fig. 11 we present a comparison between the results of the RFA method
with the PY theory and simulation data [50] for the RDF of a ternary mixture.
In the case of the RFA, we have used the eCS2 contact values and the cor-
responding isothermal compressibility. The improvement of the RFA over the
PY prediction, particularly in the region near contact, is noticeable. Although
the RFA accounts nicely for the observed oscillations, it seems to somewhat
overestimate the depth of the first minimum.
Explicit knowledge of Gij(s) also allows us to determine the Fourier trans-
form h̃ij(q) through the relation
h̃ij(q) = −2π
Gij(s)−Gij(−s)
. (161)
The structure factor Sij(q) may be expressed in terms of h̃ij(q) as [4]
Sij(q) = xiδij + ρxixj h̃ij(q). (162)
In the particular case of a binary mixture, rather than the individual structure
factors Sij(q), it is some combination of them which may be easily associated
with fluctuations of the thermodynamic variables [91, 92]. Specifically, the
quantities [4]
Snn(q) = S11(q) + S22(q) + 2S12(q), (163)
42 M. López de Haro, S. B. Yuste and A. Santos
2 3 4 5 6 7
2 3 4 5 6 7 8
3 4 5 6 7 8
10.750.50.250
10.750.50.250
1 2 3 4 5 6 7
2 3 4 5 6 7 8
3 4 5 6 7 8 9
Fig. 11. Radial distribution functions gij(r) for a ternary mixture with diameters
σ1 = 1, σ2 = 2, and σ3 = 3 at a packing fraction η = 0.49 with mole fractions
x1 = 0.7, x2 = 0.2, and x3 = 0.1. The circles are simulation results [50], the solid
lines are the RFA predictions, and the dotted lines are the PY predictions.
Snc(q) = x2S11(q)− x1S22(q) + (x2 − x1)S12(q), (164)
Scc(q) = x
2S11(q) + x
1S22(q)− 2x1x2S12(q) (165)
are sometimes required.
After replacement of ĥij(q) =
ρiρj h̃ij(q) in Eq. (147), one easily gets
c̃ij(q). Subsequent inverse Fourier transformation yields cij(r). The result
gives cij(r) for r > σij as the superposition of N Yukawas [93], namely
Alternative Approaches to Hard-Sphere Liquids 43
cij(r) =
e−κℓr
, (166)
where q = ±iκℓ with ℓ = 1, . . . , N are the zeros of det
I+ ĥ(q)
and the
amplitudes K
ij are obtained by applying the residue theorem as
q→iκℓ
c̃ij(q)(q − iκℓ). (167)
The indirect correlation functions γij(r) ≡ hij(r) − cij(r) readily follow
from the previous results for the RDF and direct correlation functions. Finally,
in this case the bridge functions Bij(r) for r > σij are linked to gij(r) and
cij(r) through
Bij(r) = ln gij(r) − γij(r) (168)
and so once more we have a full set of analytical results for the structural
properties of a multicomponent fluid mixture of HS once the contact values
gij(σij) are specified.
4 Other Related Systems
The philosophy behind the RFA method to derive the structural properties
of three-dimensional HS systems can be adapted to deal with other related
systems. The main common features of the RFA can be summarized as follows.
First, one chooses to represent the RDF in Laplace space. Next, using as a
guide the low-density form of the Laplace transform, an auxiliary function is
defined which is approximated by a rational or a rational-like form. Finally,
the coefficients are determined by imposing some basic consistency conditions.
In this section we consider the cases of sticky-hard-sphere, square-well, and
hard-disk fluids. In the two former cases the RFA program is followed quite
literally, while in the latter case it is done more indirectly through the RFA
method as applied to hard rods (d = 1) and hard spheres (d = 3).
4.1 Sticky Hard Spheres
The sticky-hard-sphere (SHS) fluid model has received a lot of attention since
it was first introduced by Baxter in 1968 [94] and later extended to multi-
component mixtures by Perram and Smith [95] and, independently, by Bar-
boy [96]. In this model, the molecular interaction may be defined via square-
well (SW) potentials of infinite depth and vanishing width, thus embodying
the two essential characteristics of real molecular interactions, namely a harsh
repulsion and an attractive part. In spite of their known shortcomings [97], an
important feature of SHS systems is that they allow for an exact solution of the
44 M. López de Haro, S. B. Yuste and A. Santos
OZ equation in the PY approximation [94,95]. Furthermore, they are thought
to be appropriate for describing structural properties of colloidal systems, mi-
celles, and microemulsions, as well as some aspects of gas-liquid equilibrium,
ionic fluids and mixtures, solvent mediated forces, adsorption phenomena,
polydisperse systems, and fluids containing chainlike molecules [98–102].
Let us consider an N -component mixture of spherical particles interacting
according to the SW potential
φij(r) =
∞, r < σij ,
−ǫij , σij < r < Rij ,
0, r > Rij .
(169)
As in the case of additive HS, σij = (σi + σj)/2 is the distance between the
centers of a sphere of species i and a sphere of species j at contact. In addition,
ǫij is the well depth and Rij − σij indicates the well width. We now take the
SHS limit [94], namely
Rij → σij , ǫij → ∞, τij ≡
Rij − σij
e−ǫij/kBT = finite, (170)
where the τij are monotonically increasing functions of the temperature T and
their inverses measure the degree of “adhesiveness” of the interacting spheres i
and j. Even without strictly taking the mathematical limits (170), short-range
SW fluids can be well described in practice by the SHS model [103].
The virial EOS for the SHS mixture is given by
Z = 1 +
i,j=1
dr ryij(r)
e−φij(r)/kBT
= 1 +
i,j=1
xixjσ
ijyij(σij)
12τij
y′ij(σij)
yij(σij)
, (171)
where yij(r) ≡ gij(r)eφij(r)/kBT is the cavity function and y′ij(r) = dyij(r)/dr.
Since yij(r) must be continuous, it follows that
gij(r) = yij(r)
Θ(r − σij) +
12τij
δ(r − σij)
. (172)
The case of a HS system is recovered by taking the limit of vanishing adhe-
siveness τ−1ij → 0, in which case Eq. (171) reduces to the three-dimensional
version of Eq. (2). On the other hand, the compressibility EOS, Eq. (145), is
valid for any interaction potential, including SHS.
As in the case of HS, it is convenient to define the Laplace transform (149).
The condition yij(σij) = finite translates into the following large s behavior
of Gij(s):
eσijsGij(s) = σ
ijyij(σij)
12τij
+ σ−1ij s
+O(s−2), (173)
Alternative Approaches to Hard-Sphere Liquids 45
which differs from (151): while eσijsGij(s) ∼ s−1 for HS, eσijsGij(s) ∼ s0 for
SHS. However, the small s behavior is still given by Eq. (152), as a consequence
of the condition χ−1 = finite.
The RFA proposal for SHS mixtures [104] keeps the form (154), except
that now
Lij(s) = L
ij + L
ij s+ L
2 + L
3, (174)
Aij(s) = ρi
ϕ2(σis)σ
ij + ϕ1(σis)σ
ij + ϕ0(σis)σiL
ij − e
−σisL
(175)
instead of Eqs. (155) and (156). By construction, Eqs. (154), (174), and (175)
comply with the requirement lims→∞ e
σijsGij(s) = finite. Further, in view
of Eq. (152), the coefficients of s0 and s in the power series expansion of
s2Gij(s) must be 1 and 0, respectively. This yields 2N
2 conditions that allow
us to express L(0) and L(1) in terms of L(2), L(3), and α as [104]
ij = ϑ1+ϑ2σj+2ϑ2α−ϑ1
kj − L
ρkσkL
kj , (176)
ij = ϑ1σij +
ϑ2σiσj + (ϑ1 + ϑ2σi)α−
kj − L
(ϑ1 + ϑ2σi)
ρkσkL
kj , (177)
where ϑ1 and ϑ2 are defined below Eq. (158). We have the freedom to choose
(3) and α, but L(2) is constrained by the condition (173), i.e., the ratio
between the first and second terms in the expansion of eσijsGij(s) for large s
must be exactly equal to σij/12τij.
First-Order Approximation (PY Solution)
The simplest approximation consists of making α = 0. In view of the condition
eσijsGij(s) ∼ s0 for large s, this implies L(3)ij = 0. In that case, the large s
behavior that follows from Eq. (154) is
2πeσijsGij(s) = L
(2) · D
s−1 +O(s−2), (178)
where
Dij ≡ ρi
σ2i L
ij − σiL
ij + L
. (179)
Comparison with Eq. (173) yields
yij(σij) =
πσ2ij
ij , (180)
46 M. López de Haro, S. B. Yuste and A. Santos
12τijL
ik Dkj . (181)
Taking into account Eqs. (176) and (177) (with L
ij = L
ji and of course also
with α = 0 and L(3) = 0), Eq. (181) becomes a closed equation for L(2):
12τijL
= ϑ1σij+
ϑ2σiσj−
ki σj + L
kj σi
(182)
The physical root L(2) of Eq. (182) is the one vanishing in the HS limit τij →
∞. Once known, Eq. (180) gives the contact values.
This first-order approximation obtained from the RFA method turns out
to coincide with the exact solution of the PY theory for SHS [95].
Second-Order Approximation
As in the case of HS mixtures, a more flexible proposal is obtained by keeping
α (and, consequently, L
ij ) different from zero. In that case, instead of Eq.
(178), one has
2πeσijsGij(s) =
+O(s−2). (183)
This implies
πσ2ij
αyij(σij), (184)
12τijL
. (185)
If we fix yij(σij), Eqs. (176), (177), (184), and (185) allow one to express L
(1), L(2), and L(3) as linear functions of α. Thus, only the scalar parameter
α remains to be fixed, analogously to what happens in the HS case. As done
in the latter case, one possibility is to choose α in order to reproduce the
isothermal compressibility χ given by Eq. (148). To do so, one needs to find
the coefficients H
ij appearing in Eq. (152). The result is [104]
(0) = C(0) ·
I− A(0)
, (186)
(1) = C(1) ·
I− A(0)
, (187)
where
Alternative Approaches to Hard-Sphere Liquids 47
αδkj −A(1)kj
δkj −A(0)kj
(188)
σ2ik +H
αδkj −A(1)kj
σ3ik + σikH
δkj −A(0)kj
, (189)
ij = (−1)
σn+3i
(n+ 3)!
σn+2i
(n+ 2)!
σn+1i
(n+ 1)!
(190)
Equation (187) gives H(1) in terms of α: H
ij = Pij(α)/[Q(α)]
2, where Pij(α)
denotes a polynomial in α of degree 2N and Q(α) denotes a polynomial of
degree N . It turns out then that, seen as a function of α, χ is the ratio of
two polynomials of degree 2N . Given a value of χ, one may solve for α. The
physical solution, which has to fulfill the requirement that Gij(s) is positive
definite for positive real s, corresponds to the smallest positive real root.
Once α is known, the scheme is complete: Eq. (184) gives L(3), then L(2) is
obtained from Eq. (185), and finally L(1) and L(0) are given by Eqs. (176) and
(177), respectively. Explicit knowledge of Gij(s) through Eqs. (154), (174),
and (175) allows one to determine the Fourier transform h̃ij(q) and the struc-
ture factor Sij(q) through Eqs. (161) and (162), respectively. Finally, inverse
Laplace transformation of Gij(s) yields gij(r) [90].
Single Component SHS Fluids
The special case of single component SHS fluids [105, 106] can be obtained
from the multicomponent one by taking σij = σ and τij = τ . Thus, the
Laplace transform of rg(r) in the RFA is
G(s) =
L(0) + L(1)s+ L(2)s2 + L(3)s3
1 + αs− ρ
ϕ2(s)L(0) + ϕ1(s)L(1) + ϕ0(s)L(2) − e−sL(3)
(191)
where we have taken σ = 1. Equations (176) and (177) become
L(0) = 2π
1 + 2η
(1− η)2
− L(2)
(1 − η)2
(1 − 4η)L(3), (192)
L(1) = 2π
1 + 1
(1− η)2
1 + 2η
α− 3ηL(2)
− 18η
(1− η)2
L(3). (193)
The choice α = L(3) = 0 makes Eq. (191) coincide with the exact solution
to the PY approximation for SHS [94], where L(2) is the physical root (i.e., the
one vanishing in the limit τ → ∞) of the quadratic equation [see Eq. (182)]
48 M. López de Haro, S. B. Yuste and A. Santos
12τL(2) = 2π
1 + 2η
(1− η)2
− 12η
L(2) +
ηL(2)
. (194)
We can go beyond the PY approximation by prescribing a contact value
y(1), so that, according to Eqs. (184) and (185),
L(3) =
y(1), (195)
L(2) =
12τ +
L(3). (196)
By prescribing the isothermal compressibility χ, the parameter α can be ob-
tained as the physical solution (namely, the one remaining finite in the limit
τ → ∞) of a quadratic equation [106]. Thus, given an EOS for the SHS
fluid, one can get the thermodynamically consistent values of y(1) and χ and
determine from them all the coefficients appearing in Eq. (191).
Figure 12 shows the cavity function for η = 0.164 and τ = 0.13 as obtained
from Monte Carlo simulations [101] and as predicted by the PY and RFA
theories, the latter making use of the EOS recently proposed by Miller and
Frenkel [102]. It can be observed that the RFA is not only more accurate than
the PY approximation near r = 1 but also near r = 2. On the other hand,
none of these two approximations account for the singularities (delta-peaks
and/or discontinuities) of y(r) at r =
8/3, 5/3,
3, 2, . . . [100, 101].
1.0 1.5 2.0 2.5
=0.164, =0.13 MC
Fig. 12. Cavity function of a single component SHS fluid for η = 0.164 and τ = 0.13.
The solid line represents simulation data [101]. The dotted and dashed lines represent
the PY and RFA approaches, respectively.
Alternative Approaches to Hard-Sphere Liquids 49
4.2 Single Component Square-Well Fluids
Now we consider again the SW interaction potential (169) but for a single
fluid, i.e., σij = σ, ǫij = ǫ, Rij = R. Since no exact solution of the PY theory
for the SW potential is known, the application of the RFA method is more
challenging in this case than for HS and SHS fluids.
As in the cases of HS and SHS, the key quantity is the Laplace transform
of rg(r) defined by Eq. (98). It is again convenient to introduce the auxiliary
function Ψ(s) through Eq. (99). As before, the conditions g(r) = finite and
χ = finite imply Eqs. (102) and (104), respectively. However, the important
difference between HS and SHS fluids is that in the latter case G(s) must
reflect the fact that g(r) is discontinuous at r = R as a consequence of the
discontinuity of the potential φ(r) and the continuity of the cavity function
y(r). This implies that G(s), and hence Ψ(s), must contain the exponential
term e−(R−σ)s. This manifests itself in the low-density limit, where the con-
dition limρ→0 y(r) = 1 yields
Ψ(s) =
(1 + s)− e−(R−1)s(e1/T∗ − 1)(1 +Rs)
, (197)
where T ∗ ≡ kBT/ǫ and we have taken σ = 1.
In the spirit of the RFA method, the simplest form that complies with Eq.
(102) and is consistent with Eq. (197) is [107]
Ψ(s) =
−12η + E1s+ E2s2 + E3s3
1 +Q0 +Q1s− e−(R−1)s (Q0 +Q2s)
, (198)
where the coefficients Q0, Q1, Q2, E1, E2, and E3 are functions of η, T
∗, and
R. The condition (104) allows one to express the parameters Q1, E1, E2, and
E3 as linear functions of Q0 and Q2 [107, 108]:
1 + 2η
+ 2η(R3 − 1)Q2 −
(R − 1)2(R2 + 2R+ 3)Q0
+Q2 − (R − 1)Q0, (199)
1 + 2η
3− 4(R3 − 1)Q2 + (R− 1)2(R2 + 2R+ 3)Q0
, (200)
1 + 2η
{1− η − 2(R− 1) [1− 2ηR(R+ 1)]Q2
+(R− 1)2
(1− η(R + 1)2
, (201)
1 + 2η
(1 − η)2 + 6η(R− 1)
R+ 1− 2ηR2
−η(R − 1)2[4 + 2R− η(3R2 + 2R+ 1)]Q0
. (202)
50 M. López de Haro, S. B. Yuste and A. Santos
From Eq. (102), we have
g(1+) =
. (203)
The complete RDF is given by Eq. (122), where now Eq. (198) must be used
in Eq. (123). In particular, ψ1(r) and ψ2(r) are
ψ1(r) = ψ10(r)Θ(r) + ψ11(r + 1−R)Θ(r + 1−R), (204)
ψ2(r) = ψ20(r)Θ(r)+ψ21(r+1−R)Θ(r+1−R)+ψ22(r+2−2R)Θ(r+2−2R),
(205)
where
ψ1k(r) = 2π
W1k(si)
E′(si)
six, (206)
ψ2k(r) = −4π2
rW2k(si) +W
2k(si)−W2k(si)
E′′(si)
E′(si)
[E′(si)]2
. (207)
Here, si are the three distinct roots of E(s) ≡ −12η+E1s+E2s2 +E3s3 and
W10(s) ≡ 1 +Q0 +Q1s, W11(s) ≡ −(Q0 +Q2s). (208)
W20(s) ≡ s[W10(s)]2, W21(s) ≡ 2sW10(s)W11(s), W22(s) ≡ s[W11(s)]2.
(209)
To close the proposal, we need to determine the parameters Q0 and Q2 by
imposing two new conditions. An obvious condition is the continuity of the
cavity function at r = R, what implies
g(R+) = e1/T
g(R−). (210)
This yields (
1− e−1/T
ψ10(R− 1) = −ψ11(0) = 2π
. (211)
As an extra condition, we could enforce the continuity of the first derivative
y′(r) at r = R [109]. However, this complicates the problem too much without
any relevant gain in accuracy. In principle, it might be possible to impose
consistency with a given EOS, via either the virial route, the compressibility
route, or the energy route. But this is not practical since no simple EOS for
SW fluids is at our disposal for wide values of density, temperature, and range.
As a compromise between simplicity and accuracy, we fix the parameter Q0
at its exact zero-density limit value, namely Q0 = e
1/T∗ − 1 [107]. Therefore,
Eq. (211) becomes a transcendental equation for Q2 that needs to be solved
numerically. For narrow SW potentials, however, it is possible to replace the
exact condition (210) by a simpler one allowing Q2 to be obtained analytically
[108], which is especially useful for determining the thermodynamic properties
[108, 110].
Alternative Approaches to Hard-Sphere Liquids 51
R=1.5, 3=0.4, T*=1.5
1.0 1.5 2.0 2.5 3.0
2.0 R=2, 3=0.4, T*=3
R=1.05, 3=0.8, T*=0.5
Fig. 13. Radial distribution function of a single component SW fluid for R = 1.05,
ρσ3 = 0.8, and T ∗ = 0.5 (top panel), for R = 1.5, ρσ3 = 0.4, and T ∗ = 1.5 (middle
panel), and for R = 2.0, ρσ3 = 0.4, and T ∗ = 3.0 (bottom panel). The circles
represent simulation data [111] and the solid lines refer to the results obtained from
the RFA method.
It can be proven that the RFA proposal (198) reduces to the exact solutions
of the PY equation [29, 94] in the HS limit, i.e., ǫ → 0 or R → 1, and in the
SHS limit, i.e., ǫ→ ∞ and R → 1 with (R− 1)e1/T∗ = finite [107, 108].
Comparison with computer simulations [107, 108, 110, 111] shows that the
RFA for SW fluids is rather accurate at any fluid density if the potential well
is sufficiently narrow (say R ≤ 1.2), as well as for any width if the density
is small enough (say ρσ3 ≤ 0.4). However, as the width and/or the density
increase, the RFA predictions worsen, especially at low temperatures. As an
52 M. López de Haro, S. B. Yuste and A. Santos
illustration, Fig. 13 compares the RDF provided by the RFA with Monte Carlo
data [111] for three representative cases.
4.3 Hard Disks
As is well known, the PY theory is exactly solvable for HS fluids with an
odd number of dimensions [112–114]. In particular, in the case of hard rods
(d = 1), the PY theory provides the exact RDF g(r) or, equivalently, the exact
cavity function y(r) outside the hard core (i.e., for r > σ). However, it does
not reproduce the exact y(r) in the overlapping region (i.e., for r < σ) [85].
The full exact one-dimensional cavity function is [85]
yHR(r|η) =
e−(r−1)η/(1−η)
ηn−1e−(r−n)η/(1−η)
(1 − η)n(n− 1)!
(r − n)n−1Θ(r − n),
(212)
where the subscript HR stands for hard rods and, as usual, σ = 1 has been
taken. Consequently, one has
gHR(1
+|η) = 1
dr rhHR(r|η) ≡ H(0)HR(η) = −
η2. (213)
When d is even, the PY equation is not analytically solvable for the HS
interaction. In particular, in the important case of hard disks (d = 2), one
must resort to numerical solutions of the PY equation [1, 115]. Alternatively,
a simple heuristic approach has proven to yield reasonably good results [116].
Such an approach is based on the näıve assumption that the structure and
spatial correlations of a hard-disk fluid share some features with those of a
hard-rod and a hard-sphere fluid. This fuzzy idea becomes a more specific one
by means of the following simple model [116]:
gHD(r|η) = ν(η)gHR(r|ω1(η)η) + [1− ν(η)]gHS(r|ω3(η)η). (214)
Here, the subscript HD stands for hard disks (d = 2) and the subscript HS
stands for hard spheres (d = 3). The parameter ν(η) is a density-dependent
mixing parameter, while ω1(η)η and ω3(η)η are the packing fractions in one
and three dimensions, respectively, which are “equivalent” to the packing
fraction η in two dimensions. In Eq. (214), it is natural to take for gHR(r|η)
the exact solution, Eq. (212). As for gHR(r|η), one might use the RFA recipe
described in Section 3. However, in order to keep the model (214) as simple
as possible, it is sufficient for practical purposes to take the PY solution, Eq.
(108). In the latter approximation,
gHS(1
+|η) = 1 + η/2
(1− η)2
dr rhHS(r|η) ≡ H(0)HS (η) = −
10− 2η + η2
20(1 + 2η)
(215)
Alternative Approaches to Hard-Sphere Liquids 53
In order to close the model (214), we still need to determine the parameters
ν(η), ω1(η), and ω3(η). To that end, we first impose the condition that Eq.
(214) must be consistent with a prescribed contact value gHD(1
+|η) or, equiv-
alently, with a prescribed compressibility factor ZHD(η) = 1 + 2ηgHD(1
+|η),
with independence of the choice of the mixing parameter ν(η). In other words,
gHD(1
+|η) = gHR(1+|ω1(η)η) = gHS(1+|ω3(η)η). (216)
Making use of Eqs. (213) and (215), this yields
ω1(η) =
gHD(1
+|η)− 1
ηgHD(1+|η)
, ω3(η) =
4gHD(1
+|η) + 1−
24gHD(1+|η) + 1
4ηgHD(1+|η)
(217)
Once ω1(η) and ω3(η) are known, we can determine ν(η) by imposing that
the model (214) reproduces the isothermal compressibility χHD(η) thermody-
namically consistent with the prescribed ZHD(η) [cf. Eq. (97)]. From Eqs. (96)
and (214) one has
χHD(η) = 1 + 8η
dr r {ν(η)hHR(r|ω1(η)η) + [1− ν(η)] hHS(r|ω3(η)η)} ,
(218)
so that
ν(η) =
[χHD(η)− 1] /8η −H(0)HS (ω3(η)η)
HR(ω1(η)η) −H
HS (ω3(η)η)
, (219)
where H
HR(η) and H
HS (η) are given by Eqs. (213) and (215), respectively.
Once a sensible EOS for hard disks is chosen [see, for instance, Table 1],
Eqs. (217) and (219) provide the parameters of the model (214). The results
show that the scaling factor ω1(η) is a decreasing function, while ω3(η) is
an increasing function [116]. As for the mixing parameter ν(η), it is hardly
dependent of density and takes values around ν(η) ≃ 0.35–0.40.
Comparison of the interpolation model (214) with computer simulation re-
sults shows a surprisingly good agreement, despite the crudeness of the model
and the absence of empirical fitting parameters, especially at low and mod-
erate densities [116]. The discrepancies become important only for distances
beyond the location of the second peak and for densities close to the stability
threshold.
5 Perturbation Theory
When one wants to deal with realistic intermolecular interactions, the prob-
lem of deriving the thermodynamic and structural properties of the system
becomes rather formidable. Thus, perturbation theories of liquids have been
devised since the mid twentieth century. In the case of single component flu-
ids, the use of an accurate and well characterized RDF for the HS fluid in a
54 M. López de Haro, S. B. Yuste and A. Santos
perturbation theory opens up the possibility of deriving a closed theoretical
scheme for the determination of the thermodynamic and structural proper-
ties of more realistic models, such as the Lennard–Jones (LJ) fluid. In this
section, we will consider this model system, which captures the basic physical
properties of real non-polar fluids, to illustrate the procedure.
In the application of the perturbation theory of liquids, the stepping stone
has been the use of the HS RDF obtained from the solution to the PY equa-
tion. Unfortunately, the absence of thermodynamic consistency present in the
PY approximation (as well as in other integral equation theories) may clearly
contaminate the results derived from its use within a perturbative treatment.
In what follows we will reanalyze the different theoretical schemes for the ther-
modynamics of LJ fluids that have been constructed with perturbation theory,
taking as the reference system the HS fluid. This includes the consideration
of the RDF as obtained with the RFA method, which embodies thermody-
namic consistency, as well as the proposal of a unifying framework in which
all schemes fit in. With our development, we will be able to present a formula-
tion which lends itself to relatively easy numerical calculations while retaining
the merits that analytical results provide, namely a detailed knowledge and
control of all the approximations involved.
Let us consider a three-dimensional fluid system defined by a pair inter-
action potential φ(r). The virial and energy EOS express the compressibility
factor Z and the excess part of the Helmholtz free energy per unit volume
f ex, respectively, in terms of the RDF of the system as
Z = 1− 2
∂φ(r)
g(r)r3, (220)
= 2πρβ
dr φ(r)g(r)r2 , (221)
where β ≡ 1/kBT . Let us now assume that φ(r) is split into a known (ref-
erence) part φ0(r) and a perturbation part φ1(r). The usual perturbative
expansion for the Helmholtz free energy to first order in β leads to [2]
+ 2πρβ
dr φ1(r)g0(r)r
, (222)
where f0 and g0(r) are the free energy and the RDF of the reference system,
respectively.
The LJ potential is
φLJ(r) = 4ǫ
r−12 − r−6
, (223)
where ǫ is the depth of the well and, for simplicity, we have taken the distance
at which the potential vanishes as the length unit, i.e., φLJ(r = 1) = 0. For
this potential the reference system may be forced to be a HS system, i.e., one
can set
Alternative Approaches to Hard-Sphere Liquids 55
φ0(r) = φHS(r) =
∞, r ≤ σ0,
0, r > σ0,
(224)
where σ0 is a conveniently chosen effective HS diameter. In this case the
Helmholtz free energy to this order is approximated by
≈ fHS
+ 2πρβ
dr φLJ(r)gHS(r/σ0)r
2. (225)
Note that Eq. (225) may be rewritten in terms of the Laplace transform G(s)
of (r/σ0)gHS(r/σ0) as
≈ fHS
+ 2πρβσ30
ds ΦLJ(s)G(s), (226)
where ΦLJ(s) satisfies
rφLJ(r) = σ0
ds e−rs/σ0ΦLJ(s), (227)
so that
ΦLJ(s) = 4ǫσ
(s/σ0)
− (s/σ0)
. (228)
Irrespective of the value of the diameter σ0 of the reference system, the
right hand side of Eq. (226) represents always an upper bound for the value of
the free energy of the real system. Therefore, it is natural to determine σ0 so
as to provide the least upper bound. This is precisely the variational scheme of
Mansoori and Canfield [117,118] and Rasaiah and Stell [119], usually referred
to as MC/RS, and originally implemented with the PY theory for G(s), Eq.
(108). In our case, however, we will considerG(s) as given by the RFA method,
Eq. (113). Therefore, at fixed ρ and β, the effective diameter σ0 in the MC/RS
scheme is obtained from the conditions
{∫ η0
ZHS(η)− 1
+ 48βǫσ−20
dsG(s|η0)
(s/σ0)
− (s/σ0)
= 0, (229)
{∫ η0
ZHS(η)− 1
+ 48βǫσ−20
dsG(s|η0)
(s/σ0)
− (s/σ0)
> 0. (230)
In these equations, use has been made of the thermodynamic relationship
between the free energy and the compressibility factor, Eq. (78). Moreover, we
have called η0 ≡ (π/6)ρσ30 and have made explicit with the notation G(s|η0)
the fact that the HS RDF depends on the packing fraction η0.
56 M. López de Haro, S. B. Yuste and A. Santos
Even if the reference system is not forced to be a HS fluid, one can still use
Eq. (226) provided an adequate choice for σ0 is made such that the expansion
involved in the right hand side of Eq. (222) yields the right hand side of Eq.
(226) to order β2. This is the idea of the Barker and Henderson [120] first
order perturbation scheme (BH1), where the effective HS diameter is
1− e−βφLJ(r)
. (231)
The same ideas may be carried out to higher order in the perturbation
expansion. The inclusion of the second order term in the expansion yields the
so-called macroscopic compressibility approximation [2] for the free energy,
namely
+ 2πρβ
dr φ1(r)g0(r)r
−πρβ2χ0
dr φ21(r)g0(r)r
, (232)
where χ0 is the (reduced) isothermal compressibility of the reference system
[121].
To implement a particular perturbation scheme in this approximation un-
der a unifying framework that eventually leads to easy numerical evaluation,
two further assumptions may prove convenient. First, the perturbation poten-
tial φ1(r) ≡ φLJ(r)−φ0(r) may be split into two parts using some “molecular
size” parameter ξ ≥ σ0 such that
φ1(r) =
φ1a(r), 0 ≤ r ≤ ξ,
φ1b(r), r > ξ.
(233)
Next, a choice for the RDF for the reference system is done in the form
g0 (r) ≈ θ(r)yHS(r/σ0), (234)
where yHS is the cavity (background) correlation function of the HS system
and θ(r) is a step function defined by
θ(r) =
θa(r), 0 ≤ r ≤ ξ,
θb(r), r > ξ,
(235)
in which the functions θa(r) and θb(r) depend on the scheme.
With these assumptions the integrals involved in Eq. (232) may be rewrit-
ten as
dr φn1 (r)g0(r)r
dr φn1a(r)θa(r)yHS(r/σ0)r
dr φn1a(r)θa(r)gHS(r/σ0)r
dr φn1b(r)θb(r)gHS(r/σ0)r
2, (236)
Alternative Approaches to Hard-Sphere Liquids 57
with n = 1, 2 and where the fact that yHS(r/σ0) = gHS(r/σ0) when r > σ0
has been used. Decomposing the last integral as
and applying
the same step as in Eq. (226), Eq. (236) becomes
In = σ
ds Φnb(s)G(s) +
dr φn1a(r)θa(r)yHS(r/σ0)r
dr [φn1a(r)θa(r) − φn1b(r)θb(r)] gHS(r/σ0)r2, (237)
where the functions Φ1b(s) and Φ2b(s) are defined by the relation
rφn1b(r)θb(r) = σ0
ds e−rs/σ0Φnb(s). (238)
In the Barker–Henderson second order perturbation scheme (BH2), one
takes
θa(r) = 0, θb(r) = 1, ξ = σ0, φ1a(r) = 0, φ1b(r) = 4ǫ
r−12 − r−6
(239)
and σ0 is computed according to Eq. (231). This choice ensures that
+ 2πρβ
dr φ1(r)gHS(r/σ0)r
−πρβ2χHS
dr φ21(r)gHS(r/σ0)r
. (240)
On the other hand, if one chooses
θa(r) = exp [−β (φLJ(r) + ǫ)] , θb(r) = 1, ξ = 21/6, (241)
φ1a(r) = −ǫ, φ1b(r) = 4ǫ
r−12 − r−6
, (242)
the scheme leads to the Weeks–Chandler–Andersen (WCA) theory [122] if one
determines the HS diameter through the condition χ0 = χHS [123], which in
turn implies
dr r2e−βφ0(r)yHS(r/σ0) =
∫ 21/6
dr r2gHS(r/σ0)
1− e−βφ0(r)
. (243)
To close the scheme, the HS cavity function has to be provided in the range 0 ≤
r ≤ σ0. Fortunately, relatively simple expressions for yHS(r/σ0) are available
in the literature [124–126], apart from our own proposal, Eq. (140).
Note that θb(r) and φ1b(r), and thus also Φnb(s), are the same functions in
the BH2 and WCA schemes. It is convenient, in order to have all the quantities
needed to evaluate fLJ in these schemes, to provide explicit expressions for
Φ1b(s) and Φ2b(s). These are given by [cf. Eq. (228)]
58 M. López de Haro, S. B. Yuste and A. Santos
Φ1b(s) = ΦLJ(s), (244)
Φ2b(s) = 16ǫ
2σ−20
(s/σ0)
− 2(s/σ0)
(s/σ0)
. (245)
Up to this point, we have embodied the most popular perturbation schemes
within a unified framework that requires as input only the EOS of the HS fluid
in order to compute the Helmholtz free energy of the LJ system and leads
to relatively easy numerical computations. It should be clear that a variety
of other possible schemes, requiring the same little input, fit in our unified
framework, which is based on the RFA method for gHS(r/σ0) and G(s). Once
fLJ has been determined, the compressibility factor of the LJ fluid at a given
order of the perturbation expansion readily follows from Eqs. (222) or (232)
through the thermodynamic relation
ZLJ = ρ
. (246)
Taking into account that the HS fluid presents a fluid-solid transition at
a freezing packing fraction ηf ≃ 0.494 [127] and a solid-fluid transition at
a melting packing fraction ηm ≃ 0.54 [127], the fluid-solid and solid-fluid
coexistence lines for the LJ system may be computed from the values (ρ, T )
determined from the conditions (π/6)ρσ30(ρ, T ) = ηf and (π/6)ρσ
0(ρ, T ) = ηm,
respectively, with the effective diameter σ0(ρ, T ) obtained using any of the
perturbative schemes. Similarly, admitting that there is a glass transition in
the HS fluid at the packing fraction ηg ≃ 0.56 [128], one can now determine
the location of the liquid-glass transition line for the LJ fluid in the (ρ, T )
plane from the simple relationship (π/6)ρσ30(ρ, T ) = ηg. With a proper choice
for ZHS, it has been shown [76,129,130] that the critical point, the structure,
and the phase diagram (including a glass transition) of the LJ fluid may be
adequately described with this approach.
6 Perspectives
In this chapter we have given a self-contained account of a simple (mostly
analytical) framework for the study of the thermodynamic and structural
properties of hard-core systems. Whenever possible, the developments have
attempted to cater for mixtures with an arbitrary number of components (in-
cluding polydisperse systems) and arbitrary dimensionality. We started con-
sidering the contact values of the RDF because they enter directly into the
EOS and are required as input in the RFA method to compute the structural
properties. With the aid of consistency conditions, we were able to devise var-
ious approximate proposals which, when used in conjunction with a sensible
choice for the contact value of the RDF of the single component fluid (re-
quired in the formulation but otherwise chosen at will), have been shown to
Alternative Approaches to Hard-Sphere Liquids 59
be in reasonably good agreement with simulation results and lead to accurate
EOS both for additive and non-additive mixtures. Some aspects of the results
that follow from the use of these EOS were illustrated by looking at demixing
problems in these mixtures, including the far from intuitive case of a binary
mixture of non-additive hard spheres in infinite dimensionality.
After that, restricting ourselves to three-dimensional systems, we described
the RFA method as applied to a single component hard-sphere fluid and to
a multicomponent mixture of HS. Using this approach, we have been able to
obtain explicit analytical results for the RDF, the direct correlation function,
the static structure factor, and the bridge function, in the end requiring as
input only the contact value of the RDF of the single component HS fluid
(or equivalently its compressibility factor). One of the nice assets of the RFA
approach is that it eliminates the thermodynamic consistency problem which
is present in most of the integral equation formulations for the computation
of structural quantities. Once again, when a sensible choice for the single
component EOS is made, we have shown, through the comparison between
the results of the RFA approach and simulation data for some illustrative
cases, the very good performance of our development. Also, the use of the
RFA approach in connection with some other related systems (sticky hard
spheres, square-well fluids, and hard disks) has been addressed.
The final part of the chapter concerns the use of HS results for more
realistic intermolecular potentials in the perturbation theory of liquids. In
this instance we have been able to provide a unifying scheme in which the
most popular perturbation theory formulations may be expressed and which
was devised to allow for easy computations. We illustrated this for a LJ fluid
but it should be clear that a similar approach might be followed for other fluids
and in fact it has recently been done in connection with the glass transition
of hard-core Yukawa fluids [131].
Finally, it should be clear that there are many facets of the equilibrium and
structural properties of hard-core systems that may be studied with a simi-
lar approach but that up to now have not been considered. For instance, the
generalizations of the RFA approach for systems such as hard hyperspheres,
non-additive hard spheres, square-well mixtures, penetrable spheres [132], or
the Jagla potential [133] appear as interesting challenges. Similarly, the ex-
tension of the perturbation theory scheme to the case of LJ mixtures seems
a worthwhile task. We hope to address some of these problems in the future
and would be very much rewarded if some others were taken up by researchers
who might find these developments also a valuable tool for their work.
References
1. J. A. Barker and D. Henderson, Rev. Mod. Phys. 48, 587 (1976).
2. D. A. McQuarrie, Statistical Mechanics (Harper & Row, N. Y., 1976).
3. H. L. Friedman, A Course in Statistical Mechanics (Prentice Hall, Englewood
Cliffs, 1985).
60 M. López de Haro, S. B. Yuste and A. Santos
4. J.-P. Hansen and I. R. McDonald, Theory of Simple Liquids, (Academic Press,
London, 1986).
5. J. L. Lebowitz and D. Zomick, J. Chem. Phys. 54, 3335 (1971).
6. J. T. Jenkins and F. Mancini, J. Appl. Mech. 54, 27 (1987).
7. C. Barrio and J. R. Solana, J. Chem. Phys. 115, 7123 (2001); 117, 2451(E)
(2002).
8. J. L. Lebowitz, Phys. Rev. A 133, 895 (1964).
9. H. Reiss, H. L. Frisch, and J. L. Lebowitz, J. Chem. Phys. 31, 369 (1959); E.
Helfand, H. L. Frisch, and J. L. Lebowitz, J. Chem. Phys. 34, 1037 (1961); J.
L. Lebowitz, E. Helfand, and E. Praestgaard, J. Chem. Phys. 43, 774 (1965).
10. M. J. Mandell and H. Reiss, J. Stat. Phys. 13, 113 (1975).
11. Y. Rosenfeld, J. Chem. Phys. 89, 4272 (1988).
12. M. Heying and D. S. Corti, J. Phys. Chem. B 108, 19756 (2004).
13. T. Boubĺık, J. Chem. Phys. 53, 471 (1970).
14. E. W. Grundke and D. Henderson, Mol. Phys. 24, 269 (1972).
15. L. L. Lee and D. Levesque, Mol. Phys. 26, 1351 (1973).
16. G. A. Mansoori, N. F. Carnahan, K. E. Starling, and J. T. W. Leland, J. Chem.
Phys. 54, 1523 (1971).
17. D. Henderson, A. Malijevský, S. Lab́ık, and K. Y. Chan, Mol. Phys. 87, 273
(1996).
18. D. H. L. Yau, K.-Y. Chan, and D. Henderson, Mol. Phys. 88, 1237 (1996); 91,
1137 (1997).
19. D. Henderson and K. Y. Chan, J. Chem. Phys. 108, 9946 (1998); Mol. Phys.
94, 253 (1998); 98, 1005 (2000).
20. D. Henderson, D. Boda, K. Y. Chan, and D. T. Wasan, Mol. Phys. 95, 131
(1998).
21. D. Matyushov, D. Henderson, and K.-Y. Chan, Mol. Phys. 96, 1813 (1999).
22. D. Cao, K.-Y. Chan, D. Henderson, and W. Wang, Mol. Phys. 98, 619 (2000).
23. D. V. Matyushov and B. M. Ladanyi, J. Chem. Phys. 107, 5815 (1997).
24. C. Barrio and J. R. Solana, J. Chem. Phys. 113, 10180 (2000).
25. D. Viduna and W. R. Smith, Mol. Phys. 100, 2903 (2002); J. Chem. Phys.
117, 1214 (2002).
26. D. Henderson, Mol. Phys. 30, 971 (1975).
27. A. Santos, M. López de Haro, and S. B. Yuste, J. Chem. Phys. 103, 4622
(1995); M. López de Haro, A. Santos, and S. B. Yuste, Eur. J. Phys. 19, 281
(1998).
28. S. Luding, Phys. Rev. E 63, 042201 (2001); S. Luding, Adv. Compl. Syst.
4, 379 (2002); S. Luding and O. Strauß, in Granular Gases, T. Pöschel and
S. Luding, eds. (LNP 564, Springer-Verlag, Berlin, 2001), pp. 389–409.
29. M. S. Wertheim, Phys. Rev. Lett. 10, 321 (1963); E. Thiele, J. Chem. Phys.
39, 474 (1963).
30. N. F. Carnahan and K. E. Starling, J. Chem. Phys. 51, 635 (1969).
31. M. Luban and J. P. J. Michels, Phys. Rev. A 41, 6796 (1990).
32. E. Hamad, J. Chem. Phys. 101, 10195 (1994).
33. C. Vega, J. Chem. Phys. 108, 3074 (1998).
34. N. M. Tukur, E. Z. Hamad, and G. A. Mansoori, J. Chem. Phys. 110, 3463
(1999).
35. A. Santos, S. B. Yuste, and M. López de Haro, Mol. Phys. 96, 1 (1999).
36. A. Malijevský and J. Veverka, Phys. Chem. Chem. Phys. 1, 4267 (1999).
Alternative Approaches to Hard-Sphere Liquids 61
37. A. Santos, S. B. Yuste, and M. López de Haro, Mol. Phys. 99, 1959 (2001).
38. M. González-Melchor, J. Alejandre, and M. López de Haro, J. Chem. Phys.
114, 4905 (2001).
39. M. López de Haro, S. B. Yuste, and A. Santos, Phys. Rev. E 66, 031202
(2002).
40. A. Santos, Mol. Phys. 96, 1185 (1999); 99, 617(E) (2001).
41. C. Regnaut, A. Dyan, and S. Amokrane, Mol. Phys. 99, 2055 (2001); 100,
2907(E) (2002).
42. A. Santos, S. B. Yuste, and M. López de Haro, J. Chem. Phys. 117, 5785
(2002).
43. A. Santos, S. B. Yuste and M. López de Haro, J. Chem. Phys. 123, 234512
(2005); M. López de Haro, S. B. Yuste, and A. Santos, Mol. Phys. 104, 3461
(2006).
44. S. Luding and A. Santos, J. Chem. Phys. 121, 8458 (2004).
45. M. Barošová, A. Malijevský, S. Lab́ık, and W. R. Smith, Mol. Phys. 87, 423
(1996).
46. H. Hansen-Goos and R. Roth, J. Chem. Phys. 124, 154506 (2006).
47. R. Evans, in Liquids and Interfaces, edited by J. Charvolin, J. F. Joanny, and
J. Zinn-Justin (North-Holland, Amsterdam, 1990).
48. Y. Rosenfeld, Phys. Rev. Lett. 63, 980 (1989).
49. A. Malijevský, M. Barošová, and W. R. Smith, Mol. Phys. 91, 65 (1997).
50. Al. Malijevský, A. Malijevský, S. B. Yuste, A. Santos, and M. López de Haro,
Phys. Rev. E 66, 061203 (2002).
51. M. Buzzacchi, I. Pagonabarraga, and N. B. Wilding, J. Chem. Phys. 121, 11362
(2004).
52. Al. Malijevský, S. B. Yuste, A. Santos, and M. López de Haro, preprint arXiv:
cond-mat/0702284.
53. F. Lado, Phys. Rev. E 54, 4411 (1996).
54. I. Prigogine and S. Lafleur, Bull. Classe Sci. Acad. Roy. Belg. 40, 484, 497
(1954).
55. S. Asakura and F. Oosawa, J. Chem. Phys. 22, 1255 (1954); J. Polym. Sci. 33,
183 (1958).
56. R. Kikuchi, J. Chem. Phys. 23, 2327 (1955).
57. P. Ballone, G. Pastore, G. Galli, and D. Gazzillo, Mol. Phys. 59, 275 (1986).
58. D. Gazzillo, G. Pastore, and S. Enzo, J. Phys.: Condens. Matter 1, 3469 (1989);
D. Gazzillo, G. Pastore, and R. Frattini, J. Phys.: Condens. Matter 2,8465
(1990).
59. A. Santos, M. López de Haro, and S. B. Yuste, J. Chem. Phys. 122, 024514
(2005).
60. E. Z. Hamad, J. Chem. Phys. 105, 3229 (1996).
61. E. Z. Hamad, J. Chem. Phys. 105, 3222 (1996).
62. H. Hammawa and E. Z. Hamad, J. Chem. Soc. Faraday Trans. 92, 4943 (1996).
63. M. Al-Naafa, J. B. El-Yakubu, and E. Z. Hamad, Fluid Phase Equilibria 154,
33 (1999).
64. J. Jung, M. S. Jhon, and F. H. Ree, J. Chem. Phys. 100, 528 (1994).
65. J. Jung, M. S. Jhon, and F. H. Ree, J. Chem. Phys. 100, 9064 (1994).
66. T. Coussaert and M. Baus, J. Chem. Phys. 109, 6012 (1998).
67. A. Yu. Vlasov and A. J. Masters, Fluid Phase Equilibria 212, 183 (2003).
68. M. López de Haro and C. F. Tejero, J. Chem. Phys. 121, 6918 (2004).
62 M. López de Haro, S. B. Yuste and A. Santos
69. S. B. Yuste, A. Santos, and M. López de Haro, Europhys. Lett. 52, 158 (2000).
70. H.-O. Carmesin, H. L. Frisch, and J. K. Percus, J. Stat. Phys. 63, 791 (1991).
71. A. Santos and M. López de Haro, Phys. Rev. E 72, 010501(R) (2005).
72. R. Roth, R. Evans, and A. A. Louis, Phys. Rev. E 64, 051202 (2001).
73. S. B. Yuste and A. Santos, Phys. Rev. A 43, 5418 (1991).
74. S. B. Yuste, M. López de Haro, and A. Santos, Phys. Rev. E 53, 4820 (1996).
75. M. Robles, M. López de Haro, A. Santos, and S. B. Yuste, J. Chem. Phys. 108,
1290 (1998).
76. M. Robles and M. López de Haro, Europhys. Lett. 62, 56 (2003).
77. M. Robles and M. López de Haro, J. Chem. Phys. 107, 4648 (1997).
78. E. Waisman, Mol. Phys. 25, 45 (1973); D. Henderson and L. Blum, Mol. Phys.
32, 1627 (1976); J. S. Høye and L. Blum, J. Stat. Phys. 16, 399 (1977).
79. A. Dı́ez, J. Largo, and J. R. Solana, J. Chem. Phys. 125, 074509 (2006).
80. J. Kolafa, S. Lab́ık, and A. Malijevský, Phys. Chem. Chem. Phys. 6, 2335
(2004). See also http://www.vscht.cz/fch/software/hsmd/ for molecular dy-
namics results of g(r).
81. A. Trokhymchuk, I. Nezbeda, J. Jirsák, and D. Henderson, J. Chem. Phys.
123, 024501 (2005).
82. M. López de Haro, A. Santos, and S. B. Yuste, J. Chem. Phys. 124, 236102
(2006).
83. L. L. Lee, J. Chem. Phys. 103, 9388 (1995); L. L. Lee, D. Ghonasgi, and E.
Lomba, J. Chem. Phys. 104, 8058 (1996); L. L. Lee and A. Malijevský, J.
Chem. Phys. 114, 7109 (2001).
84. S. Lab́ık and A. Malijevský, Mol. Phys. 53, 381 (1984).
85. Al. Malijevský and A. Santos, J. Chem. Phys. 124, 074508 (2006).
86. A. Santos and Al. Malijevský, Phys. Rev. E 75, 021201 (2007).
87. S. B. Yuste, A. Santos, and M. López de Haro, J. Chem. Phys. 108, 3683
(1998).
88. L. Blum and J. S. Høye, J. Phys. Chem. 81, 1311 (1977).
89. J. Abate and W. Whitt, Queuing Systems 10, 5 (1992).
90. A code using the Mathematica computer algebra system to obtain Gij(s)
and gij(r) with the present method is available from the web page
http://www.unex.es/eweb/fisteor/santos/filesRFA.html.
91. N. W. Ashcroft and D. C. Langreth, Phys. Rev. 156, 685 (1967).
92. A. B. Bathia and D. E. Thornton, Phys. Rev. B 8, 3004 (1970).
93. S. B. Yuste, A. Santos, and M. López de Haro, Mol. Phys. 98, 439 (2000).
94. R. J. Baxter, J. Chem. Phys. 49, 2270 (1968).
95. J. W. Perram and E. R. Smith, Chem. Phys. Lett. 35, 138 (1975).
96. B. Barboy, Chem. Phys. 11, 357 (1975); B. Barboy and R. Tenne, Chem. Phys.
38, 369 (1979).
97. G. Stell, J. Stat. Phys. 63, 1203 (1991); B. Borštnik, C. G. Jesudason, and G.
Stell, J. Chem. Phys. 106, 9762 (1997).
98. B. Barboy, J. Chem. Phys. 61, 3194 (1974).
99. J. W. Perram and E. R. Smith, Chem. Phys. Lett. 39, 328 (1975); P. T.
Cummings, J. W. Perram, and E. R. Smith, Mol. Phys. 31, 535 (1976); E. R.
Smith and J. W. Perram, J. Stat. Phys. 17, 47 (1977); J. W. Perram and E.
R. Smith, Proc. R. Soc. London A353, 193 (1977); W. G. T. Kranendonk and
D. Frenkel, Mol. Phys. 64, 403 (1988); C. Regnaut and J. C. Ravey, J. Chem.
Phys. 91, 1211 (1989); G. Stell and Y. Zhou, J. Chem. Phys. 91, 3618 (1989);
Alternative Approaches to Hard-Sphere Liquids 63
J. N. Herrera and L. Blum, J. Chem. Phys. 94, 6190 (1991); A. Jamnik, D.
Bratko, and D. J. Henderson, J. Chem. Phys. 94, 8210 (1991); S. V. G. Menon,
C. Manohar, and K. S. Rao, J. Chem. Phys.; 95, 9186 (1991); Y. Zhou and G.
Stell, J. Chem. Phys. 96, 1504 (1992); E. Dickinson, J. Chem. Soc. Faraday
Trans. 88, 3561 (1992); C. F. Tejero and M. Baus, Phys. Rev. E, 48, 3793
(1993); K. Shukla and R. Rajagopalan, Mol. Phys. 81, 1093 (1994); C. Regnaut,
S. Amokrane, and Y. Heno, J. Chem. Phys. 102, 6230 (1995); C. Regnaut, S.
Amokrane, and P. Bobola, Prog. Colloid Polym. Sci. 98, 151 (1995); Y. Zhou,
C. K. Hall, and G. Stell, Mol. Phys. 86, 1485 (1995); J. N. Herrera-Pacheco
and J. F. Rojas-Rodŕıguez, Mol. Phys. 86, 837 (1995); Y. Hu, H. Liu, and J. M.
Prausnitz, J. Chem. Phys. 104, 396 (1996); O. Bernard and L. Blum, J. Chem.
Phys. 104, 4746 (1996); L. Blum, M. F. Holovko, and I. A. Protsykevych, J.
Stat. Phys. 84, 191 (1996); S. Amokrane, P. Bobola and C. Regnaut, Prog.
Colloid Polym. Sci. 100, 186 (1996); S. Amokrane and C. Regnaut, J. Chem.
Phys. 106, 376 (1997); C. Tutschka, G. Kahl, and E. Riegler, Mol. Phys. 100,
1025 (2002); D. Gazzillo and A. Giacometti, Mol. Phys. 100, 3307 (2002); M.
A. Miller and D. Frenkel, Phys. Rev. Lett. 90, 135702 (2003); D. Gazzillo and
A. Giacometti, J. Chem. Phys. 120, 4742 (2004); R. Fantoni, D. Gazzillo, and
A. Giacometti, Phys. Rev. E 72, 011503 (2005); A. Jamnik, Chem. Phys. Lett.
423, 23 (2006).
100. A. J. Post and E. D. Glandt, J. Chem. Phys. 84, 4585 (1986); N. A. Seaton
and E. D. Glandt, J. Chem. Phys. 84, 4595 (1986); 86, 4668 (1986); 87, 1785
(1987).
101. M. A. Miller and D. Frenkel, J. Phys.: Condens. Matter 16, S4901 (2004).
102. M. A. Miller and D. Frenkel, J. Chem. Phys. 121, 535 (2004).
103. Al. Malijevský, S. B. Yuste, and A. Santos, J. Chem. Phys. 125, 074507 (2006).
104. A. Santos, S. B. Yuste, and M. López de Haro, J. Chem. Phys. 109, 6814
(1998).
105. S. B. Yuste and A. Santos, J. Stat. Phys. 72, 703 (1993).
106. S. B. Yuste and A. Santos, Phys. Rev. E 48, 4599 (1993).
107. S. B. Yuste and A. Santos, J. Chem. Phys. 101, 2355 (1994).
108. L. Acedo and A. Santos, J. Chem. Phys. 115, 2805 (2001).
109. L. Acedo, J. Stat. Phys. 99, 707 (2000).
110. J. Largo, J. R. Solana, L. Acedo, and A. Santos, Mol. Phys. 101, 2981 (2003).
111. J. Largo, J. R. Solana, S. B. Yuste, and A. Santos, J. Chem. Phys. 122, 084510
(2005).
112. C. Freasier and D. J. Isbister, Mol. Phys. 42, 927 (1981).
113. E. Leutheusser, Physica A 127, 667 (1984).
114. M. Robles, M. López de Haro, and A. Santos, J. Chem. Phys. 120, 9113 (2004).
115. D. G. Chae, F. H. Ree, and T. Ree, J. Chem. Phys. 50, 1581 (1976).
116. S. B. Yuste and A. Santos, J. Chem. Phys. 99, 2020 (1993).
117. G. A. Mansoori and F. B. Canfield, J. Chem. Phys. 51, 4958 (1969).
118. G. A. Mansoori, J. A. Provine, and F. B. Canfield, J. Chem. Phys. 51, 5295
(1969).
119. J. Rasaiah and G. Stell, Mol. Phys. 18, 249 (1970).
120. J. A. Barker and D. Henderson, J. Chem. Phys. 47, 2856 (1967).
121. The macroscopic compressibility approach is only one of the possibilities of ap-
proximation to the second order Barker–Henderson perturbation theory term.
Another successful approach is the local-compressibility approximation (see
64 M. López de Haro, S. B. Yuste and A. Santos
Ref. [2], p. 308). This expresses the free energy in terms of φ1(r) and HS
quantities.
122. J. D. Weeks, D. Chandler, and H. C. Andersen, J. Chem. Phys. 53, 149 (1971).
123. A simple algorithm to compute a rather accurate approximation for the HS
diameter σ0 in the WCA theory has been given in L. Verlet and J. J. Weis,
Phys. Rev. A 5, 939 (1972).
124. D. Henderson and E. W. Grundke, J. Chem. Phys. 63, 601 (1975).
125. J. A. Ballance and R. J. Speedy, Mol. Phys. 54, 1035 (1985).
126. Y. Zhou and G. Stell, J. Stat. Phys. 52, 1389 (1988).
127. J.-P. Hansen and L. Verlet, Phys. Rev. 184, 151 (1969).
128. R. J. Speedy, J. Chem. Phys. 100, 6684 (1994).
129. M. Robles and M. López de Haro, Phys. Chem. Chem. Phys. 3, 5528 (2001).
130. M. López de Haro and M. Robles, J. Phys.: Condens. Matt. 16, S2089 (2004).
131. M. López de Haro and M. Robles, Physica A 372, 307 (2006).
132. C. N. Likos, Phys. Rep. 348, 267 (2001).
133. E. A. Jagla, J. Chem. Phys. 111, 8980 (1999).

  An overview of some analytical approaches to the computation of the
structural and thermodynamic properties of single component and multicomponent
hard-sphere fluids is provided. For the structural properties, they yield a
thermodynamically consistent formulation, thus improving and extending the
known analytical results of the Percus-Yevick theory. Approximate expressions
for the contact values of the radial distribution functions and the
corresponding analytical equations of state are also discussed. Extensions of
this methodology to related systems, such as sticky hard spheres and
square-well fluids, as well as its use in connection with the perturbation
theory of fluids are briefly addressed.

Introduction
In the statistical thermodynamic approach to the theory of simple liquids,
there is a close connection between the thermodynamic and structural prop-
erties [1–4]. These properties depend on the intermolecular potential of the
system, which is generally assumed to be well represented by pair interactions.
The simplest model pair potential is that of a hard-core fluid (rods, disks,
spheres, hyperspheres) in which attractive forces are completely neglected. In
fact, it is a model that has been most studied and has rendered some analytical
results, although up to this day no general (exact) explicit expression for the
equation of state is available, except for the one-dimensional case. Something
similar applies to the structural properties. An interesting feature concerning
the thermodynamic properties is that in hard-core systems the equation of
state depends only on the contact values of the radial distribution functions.
In the absence of a completely analytical approach, the most popular methods
http://arxiv.org/abs/0704.0157v1
2 M. López de Haro, S. B. Yuste and A. Santos
to deal with both kinds of properties of these systems are integral equation
theories and computer simulations.
It is well known that in real gases and liquids at high temperatures the
state and thermodynamic properties are determined almost entirely by the
repulsive forces among molecules. At lower temperatures, attractive forces
become significant, but even in this case they affect very little the configu-
ration of the system at moderate and high densities. These facts are taken
into account in the application of the perturbation theory of fluids, where
hard-core fluids are used as the reference systems in the computation of the
thermodynamic and structural properties of real fluids. However, successful
results using perturbation theory are rather limited due to the fact that, as
mentioned above, there are in general no exact (analytical) expressions for
the thermodynamic and structural properties of the reference systems which
are in principle required in the calculations. On the other hand, in the realm
of soft condensed matter the use of the hard-sphere model in connection, for
instance, with sterically stabilized colloidal systems is quite common. This is
due to the fact that nowadays it is possible to prepare (almost) monodisperse
spherical colloidal particles with short-ranged harshly repulsive interparticle
forces that may be well described theoretically with the hard-sphere potential.
This chapter presents an overview of the efforts we have made over the
last few years to compute the thermodynamic and structural properties of
hard-core systems using relatively simple (approximate) analytical methods.
It is structured as follows. In Section 2 we describe our proposals to derive
the contact values of the radial distribution functions of a multicomponent
mixture (with an arbitrary size distribution, either discrete or continuous) of
d-dimensional hard spheres from the use of some consistency conditions and
the knowledge of the contact value of the radial distribution function of the
corresponding single component system. In turn, these contact values lead to
equations of state both for additive and non-additive hard spheres. Some con-
sequences of such equations of state, in particular the demixing transition, are
briefly analyzed. This is followed in Section 3 by the description of the Ratio-
nal Function Approximation method to obtain analytical expressions for the
structural quantities of three-dimensional single component and multicompo-
nent fluids. The only required inputs in this approach are the contact values of
the radial distribution functions and so the connection with the work of the
previous section follows naturally. Structural properties of related systems,
like sticky hard spheres or square-well fluids, that may also be tackled with
the same philosophy are also discussed in Section 4. Section 5 provides an
account of the reformulation of the perturbation theory of liquids using the
results of the Rational Function Approximation method for a single compo-
nent hard-sphere fluid and its illustration in the case of the Lennard–Jones
fluid. In the final section, we provide some perspectives of the achievements
obtained so far and of the challenges that remain ahead.
Alternative Approaches to Hard-Sphere Liquids 3
2 Contact Values and Equations of State for Mixtures
As stated in the Introduction, a nice feature of hard-core fluids is that the
expressions of all their thermodynamic properties in terms of the radial distri-
bution functions (RDF) are particularly simple. In fact, for these systems the
internal energy reduces to that of the ideal gas and in the pressure equation
it is only the contact values rather than the full RDF which appear explicitly.
In this section we present our approach to the derivation of the contact values
of hard-core fluid mixtures in d dimensions.
2.1 Additive Systems in d Dimensions
If σij denotes the distance of separation at contact between the centers of
two interacting fluid particles, one of species i and the other of species j,
the mixture is said to be additive if σij is just the arithmetic mean of the
hard-core diameters of each species. Otherwise, the system is non-additive.
We deal in this subsection and in Subsection 2.2 with additive systems, while
non-additive hard-core mixtures will be treated in Subsection 2.3.
Definitions
Let us consider an additive mixture of hard spheres (HS) in d dimensions with
an arbitrary number N of components. In fact, our discussion will remain valid
for N → ∞, i.e., for polydisperse mixtures with a continuous distribution of
sizes.
The additive hard core of the interaction between a sphere of species i and
a sphere of species j is σij =
(σi + σj), where the diameter of a sphere of
species i is σii = σi. Let the number density of the mixture be ρ and the mole
fraction of species i be xi = ρi/ρ, where ρi is the number density of species i.
From these quantities one can define the packing fraction η = vdρMd, where
vd = (π/4)
d/2/Γ (1 + d/2) is the volume of a d-dimensional sphere of unit
diameter and
Mn ≡ 〈σn〉 =
i (1)
denotes the nth moment of the diameter distribution.
In a HS mixture, the knowledge of the contact values gij(σij) of the RDF
gij(r), where r is the distance, is important for a number of reasons. For
example, the availability of gij(σij) is sufficient to get the equation of state
(EOS) of the mixture via the virial expression
Z(η) = 1 +
i,j=1
xixjσ
ijgij(σij), (2)
4 M. López de Haro, S. B. Yuste and A. Santos
where Z = p/ρkBT is the compressibility factor of the mixture, p being the
pressure, kB the Boltzmann constant, and T the absolute temperature.
The exact form of gij(σij) as functions of the packing fraction η, the set
of diameters {σk}, and the set of mole fractions {xk} is only known in the
one-dimensional case, where one simply has [5]
gij(σij) =
, (d = 1). (3)
Consequently, for d ≥ 2 one has to resort to approximate theories or empirical
expressions. For hard-disk mixtures, an accurate expression is provided by
Jenkins and Mancini’s (JM) approximation [6, 7],
gJMij (σij) =
(1− η)2
σiσjM1
σijM2
, (d = 2). (4)
The associated compressibility factor is
ZJM(η) =
1 + η/8
(1 − η)2
, (d = 2). (5)
In the case of three-dimensional systems, some important analytical expres-
sions for the contact values and the corresponding compressibility factor also
exist. For instance, the expressions which follow from the solution of the
Percus–Yevick (PY) equation of additive HS mixtures by Lebowitz [8] are
gPYij (σij) =
(1− η)2
σiσjM2
σijM3
, (d = 3), (6)
ZPY(η) =
(1− η)2
(1 − η)2
, (d = 3). (7)
Also analytical are the results obtained from the Scaled Particle Theory (SPT)
[9–12],
gSPTij (σij) =
(1 − η)2
σiσjM2
σijM3
(1 − η)3
σiσjM2
σijM3
, (d = 3),
ZSPT(η) =
(1− η)2
(1− η)3
, (d = 3). (9)
Neither the PY nor the SPT lead to particularly accurate values and so
Boubĺık [13] and, independently, Grundke and Henderson [14] and Lee and
Levesque [15] proposed an interpolation between the PY and the SPT contact
values, that we will refer to as the BGHLL values:
gBGHLLij (σij) =
(1− η)2
σiσjM2
σijM3
(1− η)3
σiσjM2
σijM3
, (d = 3).
Alternative Approaches to Hard-Sphere Liquids 5
This leads through Eq. (2) to the widely used and rather accurate Boubĺık–
Mansoori–Carnahan–Starling–Leland (BMCSL) EOS [13,16] for HS mixtures:
ZBMCSL(η) =
(1− η)2
η2(3− η)
(1 − η)3
, (d = 3). (11)
Refinements of the BGHLL values have been subsequently introduced, among
others, by Henderson et al. [17–22], Matyushov and Ladanyi [23], and Barrio
and Solana [24] to eliminate some drawbacks of the BMCSL EOS in the so-
called colloidal limit of binary HS mixtures. On a different path, but also
having to do with the colloidal limit, Viduna and Smith [25] have proposed
a method to obtain contact values of the RDF of HS mixtures from a given
EOS. However, none of these proposals may be easily generalized so as to be
valid for any dimensionality and any number of components. Therefore, if one
wants to have a more general framework able to deal with arbitrary d and N
an alternative strategy is called for.
Universality Ansatz
In order to follow our alternative strategy, it is useful to make use of exact
limit results that can help one in the construction of approximate expressions
for gij(σij). Let us consider first the limit in which one of the species, say i,
is made of point particles, i.e., σi → 0. In that case, gii(σi) takes the ideal
gas value, except that one has to take into account that the available volume
fraction is 1− η. Thus,
gii(σi) =
. (12)
An even simpler situation occurs when all the species have the same size,
{σk} → σ, so that the system becomes equivalent to a single component
system. Therefore,
{σk}→σ
gij(σij) = gs, (13)
where gs is the contact value of the RDF of the single component fluid at
the same packing fraction η as that of the mixture. Table 1 lists some of
the most widely used proposals for the contact value gs and the associated
compressibility factor
Zs = 1 + 2
d−1ηgs (14)
in the case of the single component HS fluid.
Equations (12) and (13) represent the simplest and most basic conditions
that gij(σij) must satisfy. There is a number of other less trivial consistency
conditions [11, 17, 19, 20, 23, 24,32–34], some of which will be used later on.
In order to proceed, in line with a property shared by earlier proposals
[see, in particular, Eqs. (4), (6), (8), and (10)], we assume that, at a given
packing fraction η, the dependence of gij(σij) on the parameters {σk} and
{xk} takes place only through the scaled quantity
6 M. López de Haro, S. B. Yuste and A. Santos
Table 1. Some expressions of gs and Zs for the single component HS fluid. In the
SHY proposal, ηcp = (
3/6)π is the crystalline close-packing fraction for hard disks.
In the LM proposal, b3 and b4 are the (reduced) third and fourth virial coefficients,
ζ(η) = 1.2973(59)−0.062(13)η/ηcp for d = 4, and ζ(η) = 1.074(16)+0.163(45)η/ηcp
for d = 5, where the values of the close-packing fractions are ηcp = π
2/16 ≃ 0.617
and ηcp = π
2/30 ≃ 0.465 for d = 4 and d = 5, respectively.
d gs Zs Label Ref.
1− 7η/16
(1− η)2
1 + η2/8
(1− η)2
H [26]
1− η(2ηcp − 1)/2η2cp
1− 2η + η2(ηcp − 1)/2η2cp
1− 2η + η2(ηcp − 1)/2η2cp
SHY [27]
2 gHs −
27(1− η)4
ZHs −
26(1− η)4
L [28]
1 + η/2
(1− η)2
1 + 2η + 3η2
(1− η)2
PY [29]
1− η/2 + η2/4
(1− η)3
1 + η + η2
(1− η)3
SPT [9]
1− η/2
(1− η)3
1 + η + η2 − η3
(1− η)3
CS [30]
1 + [21−db3 − ζ(η)b4/b3]η
1− ζ(η)(b4/b3)η + [ζ(η)− 1] 21−db4η2
1 + 2d−1ηgLMs LM [31]
zij ≡
. (15)
More specifically, we assume
gij(σij) = G(η, zij), (16)
where the function G(η, z) is universal in the sense that it is a common func-
tion for all the pairs (i, j), regardless of the composition and number of compo-
nents of the mixture. Of course, the function G(η, z) is in principle different for
each dimensionality d. To clarify the implications of this universality ansatz,
let us imagine two mixtures M and M′ having the same packing fraction
η but strongly differing in the set of mole fractions, the sizes of the parti-
cles, and even the number of components. Suppose now that there exists a
pair (i, j) in mixture M and another pair (i′, j′) in mixture M′ such that
zij = zi′j′ . Then, according to Eq. (16), the contact value of the RDF for the
pair (i, j) in mixture M is the same as that for the pair (i′, j′) in mixture
M′, i.e., gij(σij) = gi′j′ (σi′j′ ). In order to ascribe a physical meaning to the
parameter zij , note that the ratio Md−1/Md can be understood as a “typ-
ical” inverse diameter (or curvature) of the particles of the mixture. Thus,
z−1ij =
(σ−1i + σ
j )/(Md−1/Md) represents the arithmetic mean curvature,
in units of Md−1/Md, of a particle of species i and a particle of species j.
Alternative Approaches to Hard-Sphere Liquids 7
Once the ansatz (16) is adopted, one may use the limits in (12) and (13)
to get G(η, z) at z = 0 and z = 1, respectively. Since zii → 0 in the limit
σi → 0, insertion of Eq. (12) into (16) yields
G(η, 0) = 1
≡ G0(η). (17)
Next, if all the diameters are equal, zij → 1, so that Eq. (13) implies that
G(η, 1) = gs. (18)
Linear Approximation
As the simplest approximation [35], one may assume a linear dependence of
G on z that satisfies the basic requirements (17) and (18), namely
G(η, z) = 1
z. (19)
Inserting this into Eq. (16), one has
ge1ij (σij) =
. (20)
Here, the label “e1” is meant to indicate that (i) the contact values used are
an extension of the single component contact value gs and that (ii) G(η, z) is
a linear polynomial in z. This notation will become handy below. Although
the proposal (20) is rather crude and does not produce especially accurate
results for gij(σij) when d ≥ 3, it nevertheless leads to an EOS that exhibits
an excellent agreement with simulations in 2, 3, 4, and 5 dimensions, provided
that an accurate gs is used as input [35–39]. This EOS may be written as
Ze1(η) = 1 +
2d−1(Ω0 −Ω1) + [Zs(η)− 1]Ω1, (21)
where the coefficients Ωm depend only on the composition of the mixture and
are defined by
Ωm = 2
−(d−m)
Mmd−1
Mm+1d
Mn+mMd−n. (22)
In particular, for d = 2 and d = 3,
Ze1(η) =
Zs(η) −
, (d = 2), (23)
8 M. López de Haro, S. B. Yuste and A. Santos
Ze1(η) =
Zs(η)−
, (d = 3). (24)
As an extra asset, from Eq. (21) one may write the virial coefficients of
the mixture Bn, defined by
Z = 1 +
Bn+1ρ
n, (25)
in terms of the (reduced) virial coefficients of the single component fluid bn
defined by
Zs = 1 +
bn+1η
n. (26)
The result is
Bn = v
Ω1bn + 2
d−1(Ω0 −Ω1)
. (27)
In the case of binary mixtures, these coefficients are in very good agreement
with the available exact and simulation results [35,37], except when the mix-
ture involves components of very disparate sizes, especially for high dimen-
sionalities. One may perform a slight modification such that this deficiency is
avoided and thus get a modified EOS [37, 40]. For d = 2 and d = 3 it reads
Z(η) = Zs(η) + x1
1− η2
1− η2
− Zs(η)
σ2 − σ1
1− η1
1− η1
− Zs(η)
σ1 − σ2
, (d = 2, 3),
where ηi = vdρiσ
i is the partial volume packing fraction due to species i. In
contrast to most of the approaches (PY, SPT, BMCSL, e1, . . . ), the proposal
(28) expresses Z(η) in terms not only of Zs(η) but also involves Zs
and Zs
. Equation (28) should in principle be useful in particular for
binary mixtures involving components of very disparate sizes. However, it is
slightly less accurate than the one given in Eq. (21) for ordinary mixtures [37].
Quadratic Approximation
In order to improve the proposal contained in Eq. (20), in addition to the
consistency requirements (12) and (13), one may consider the condition stem-
ming from a binary mixture in which one of the species (say i = 1) is much
Alternative Approaches to Hard-Sphere Liquids 9
larger than the other one (i.e., σ1/σ2 → ∞), but occupies a negligible volume
(i.e., x1(σ1/σ2)
d → 0). In that case, a sphere of species 1 is felt as a wall by
particles of species 2, so that [17, 20, 41]
σ1/σ2→∞
x1(σ1/σ2)d→0
g12(σ12)− 2d−1ηg22(σ2)
= 1. (29)
Hence, in the limit considered in Eq. (29), we have z22 → 1, z12 → 2. Conse-
quently, under the universality ansatz (16), one may rewrite Eq. (29) as
G(η, 2) = 1 + 2d−1ηG(η, 1). (30)
Thus, Eqs. (17), (18), and (30) provide complete information on the function
G at z = 0, z = 1, and z = 2, respectively, in terms of the contact value gs of
the single component RDF.
The simplest functional form of G that complies with the above consistency
conditions is a quadratic function of z [42]:
G(η, z) = G0(η) + G1(η)z + G2(η)z2, (31)
where the coefficients G1(η) and G2(η) are explicitly given by
G1(η) = (2− 2d−2η)gs −
2− η/2
1− η , (32)
G2(η) =
1− η/2
− (1− 2d−2η)gs. (33)
Therefore, the explicit expression for the contact values is
ge2ij (σij) =
(2− 2d−2η)gs −
2− η/2
1− η/2
− (1− 2d−2η)gs
. (34)
Following the same criterion as the one used in connection with Eq. (20), the
label “e2” is meant to indicate that (i) the resulting contact values represent
an extension of the single component contact value gs and that (ii) G(η, z)
is a quadratic polynomial in z. Of course, the quadratic form (31) is not the
only choice compatible with conditions (17), (18), and (30). For instance, a
rational function was also considered in Ref. [42]. However, although it is
rather accurate, it does not lead to a closed form for the EOS. In contrast,
when Eq. (34) is inserted into Eq. (2), one gets a closed expression for the
compressibility factor in terms of the packing fraction η and the first few
moments Mn, n ≤ d. The result is
Ze2(η) = 1 + 2
d−2 η
1− η [2(Ω0 − 2Ω1 +Ω2) + (Ω1 −Ω2)η]
+ [Zs(η)− 1]
2Ω1 −Ω2 + 2d−2(Ω2 −Ω1)η
, (35)
10 M. López de Haro, S. B. Yuste and A. Santos
where the quantities Ωm are defined in Eq. (22). Quite interestingly, in the
two-dimensional case Eq. (35) reduces to Eq. (23), i.e.,
Ze1(η) = Ze2(η), (d = 2). (36)
This illustrates the fact that two different proposals for the contact values
gij(σij) can yield the same EOS when inserted into Eq. (2). On the other
hand, for three-dimensional mixtures Eq. (35) becomes
Ze2(η) =
1− η + M
Zs(η)−
, (d = 3),
which differs from Eq. (24). In fact,
Ze1(η)− Ze2(η) =
1 + η
− (1 − 2η)Zs(η)
, (d = 3).
Specific Examples
In this subsection, rather than carrying out an exhaustive comparison with
the wealth of results available in the literature, we will consider only a few
representative examples. In particular, for d = 3, we will restrict ourselves
to a comparison with classical proposals (say BGHLL, PY, and SPT for the
contact values). The comparison with more recent ones may be found in Refs.
[35, 42, 43].
Thus far the development has been rather general since gs remains free in
Eqs. (20) and (34). In order to get specific results, it is necessary to fix gs [cf.
Table 1]. In the one-dimensional case, one has gs = 1/(1− η) and so one gets
the exact result (3) after substitution into Eq. (20). Similarly Eqs. (32) and
(33) lead to G1 = G2 = 0 and so we recover again the exact result.
If in the two-dimensional case we take Henderson’s value [26] gs = g
then the linear approximation (20) reduces to the JM approximation, Eq. (4).
This equivalence can be symbolically represented as geH1ij = g
ij , where the
label “eH1” refers to the extension of Henderson’s single component value in
the linear approximation. While gJMij is very accurate, even better results are
provided by the quadratic form (34), especially if Luding’s value [28] gs = g
is used [44].
In the three-dimensional case, Eq. (20) is of the form of the solution of the
PY equation [8]. In fact, insertion of gs = g
s leads to Eq. (6), i.e., g
gPYij . Similarly, if the SPT expression [9] gs = g
s is used for the single
component contact value in the quadratic approximation (34), we reobtain
the SPT expression for the mixture, Eq. (8). In other words, geSPT2ij = g
On the other hand, if the much more accurate CS [30] expression gs = g
used as input, we arrive at the following expression:
Alternative Approaches to Hard-Sphere Liquids 11
geCS2ij =
η(1 − η/3)
(1 − η)2
σiσjM2
σijM3
η2(1− η/2)
(1− η)3
σiσjM2
σijM3
, (d = 3),
which is different from the BGHLL one, Eq. (10), improves the latter for
zij > 1, and leads to similar results for zij < 1, as comparison with computer
simulations shows [42]. The four approximations (6), (8), (10), and (39) are
consistent with conditions (12) and (13), but only the SPT and eCS2 are also
consistent with condition (29). It should also be noted that if one considers
a binary mixture in the infinite solute dilution limit, namely x1 → 0, so that
z12 → 2/(1 + σ2/σ1), Eq. (39) yields the same result for g12(σ12) as the one
proposed by Matyushov and Ladanyi [23] for this quantity on the basis of
exact geometrical relations. However, the extension that the same authors
propose when there is a non-vanishing solute concentration, i.e., for x1 6= 0,
is different from Eq. (39).
Equation (34) can also be used in the case of hyperspheres (d ≥ 4) [42]. In
particular, a very good agreement with available computer simulations [38] is
obtained for d = 4 and d = 5 by using Luban and Michels [31] value gs = g
0.0 0.1 0.2 0.3 0.4 0.5
-0.02
 eCS1
 eCS2
Fig. 1. Deviation of the compressibility factor from the BMCSL value, as a function
of the packing fraction η for an equimolar three-dimensional binary mixture with
σ2/σ1 = 0.6. The open (Ref. [18]) and closed (Ref. [45]) circles are simulation data.
The lines are the PY EOS (– · · –), the SPT EOS (– · – ·), the eCS1 EOS (· · · ), and
the eCS2 EOS (– – –).
12 M. López de Haro, S. B. Yuste and A. Santos
0.00 0.05 0.10 0.15 0.20 0.25
 4D, / =1/2
 4D, / =1/3
 5D, / =2/5
Fig. 2. Compressibility factor for three equimolar mixtures in 4D and 5D systems.
Lines are the eLM1 predictions, while symbols are simulation data [38].
Now we turn to the compressibility factors (21) and (35), which are ob-
tained from the contact values (20) and (34), respectively. Since they depend
on the details of the composition through the d first moments, they are mean-
ingful even for continuous polydisperse mixtures.
As said above, in the two-dimensional case both Eqs. (21) and (35) reduce
to Eq. (23), which yield very accurate results when a good Zs is used as
input [39, 42, 44]. For three-dimensional mixtures, insertion of Zs = Z
Eqs. (24) and (37) yields
ZeCS1(η) = ZBMCSL(η) +
(1 − η)3M23
M1M3 −M22
, (d = 3), (40)
ZeCS2(η) = ZBMCSL(η)−
(1 − η)2M23
M1M3 −M22
, (d = 3), (41)
where ZBMCSL(η) is given by Eq. (11). Note that ZeCS1(η) > ZBMCSL(η) >
ZeCS2(η). Since simulation data indicate that the BMCSL EOS tends to un-
derestimate the compressibility factor, it turns out that, as illustrated in Fig.
1 for an equimolar binary mixture with σ2/σ1 = 0.6, the performance of ZeCS1
is, paradoxically, better than that of ZeCS2 [42], despite the fact that the un-
derlying linear approximation for the contact values is much less accurate than
the quadratic approximation. This shows that a rather crude approximation
such as Eq. (20) may lead to an extremely good EOS [35, 37–39], which, as
Alternative Approaches to Hard-Sphere Liquids 13
clearly seen in Fig. 1, represents a substantial improvement over the classical
proposals. Interestingly, the EOS corresponding to ZeCS1 has recently been
independently derived as the second order approximation of the Fundamental
Measure Theory for the HS fluid by Hansen-Goos and Roth [46].
In the case of d = 4 and d = 5, use of Zs(η) = Z
s (η) in Eq. (21)
produces a simple extended EOS of a mixture of hard additive hyperspheres
in these dimensionalities. The accuracy of these two EOS for hard hypersphere
mixtures in the fluid region has been confirmed by simulation data [38] for a
wide range of compositions and size ratios. In Fig. 2, this accuracy is explicitly
exhibited in the case of three equimolar mixtures, two in 4D and one in 5D.
2.2 A More Consistent Approximation for Three-Dimensional
Additive Mixtures
Up to this point, we have considered an arbitrary dimensionality d and have
constructed, under the universality assumption (16), the acurate quadratic
approximation (34), which fulfills the consistency conditions (12), (13), and
(29). However, there exist extra consistency conditions that are not necessarily
satisfied by (34). In particular, when the mixture is in contact with a hard wall,
the state of equilibrium imposes that the pressure evaluated near the wall by
considering the impacts with the wall must be the same as the pressure in the
bulk evaluated from the particle-particle collisions. This consistency condition
is especially important if one is interested in deriving accurate expressions for
the contact values of the particle-wall correlation functions.
Since a hard wall can be seen as a sphere of infinite diameter, the contact
value gwj of the correlation function of a sphere of diameter σj with the wall
can be obtained from gij(σij) as
gwj = lim
gij(σij). (42)
Note that gwj provides the ratio between the density of particles of species
j adjacent to the wall and the density of those particles far away from the
wall. The sum rule connecting the pressure of the fluid and the above contact
values is [47]
Zw(η) =
xjgwj, (43)
where the subscript w in Zw has been used to emphasize that Eq. (43) repre-
sents a route alternative to the virial one, Eq. (2), to get the EOS of the HS
mixture. The condition Z = Zw is equivalent to (29) in the special case where
one has a single fluid in the presence of the wall. However, in the general case
of a mixture plus a wall, the condition Z = Zw is stronger than Eq. (29).
In the two-dimensional case, it turns out that the quadratic approximation
(34) already satisfies the requirement Z = Zw, regardless of the density and
composition of the mixture [44]. However, this is not the case for d ≥ 3.
14 M. López de Haro, S. B. Yuste and A. Santos
Our problem now consists of computing gij(σij) and the associated gwj for
the HS mixture in the presence of a hard wall, so that the condition Z = Zw
is satisfied for an arbitrary mixture [43]. Due to the mathematical complexity
of the problem, here we will restrict ourselves to three-dimensional systems
(d = 3). Similarly to what we did in the preceding subsection, we consider
a class of approximations of the universal type (16), so that conditions (12)
and (13) lead again to Eqs. (17) and (18), respectively. Notice that Eq. (16)
implies in particular that
gwj = G(η, zwj), zwj = 2σj
. (44)
Assuming that z = 0 is a regular point and taking into account condition
(17), G(η, z) can be expanded in a power series in z:
G(η, z) = G0(η) +
Gn(η)zn. (45)
After simple algebra, using the ansatz (16) and Eq. (45) in Eqs. (2) (with
d = 3) and (43) one gets
Z = G0 + 3η
G0 + 4η
Mn+13
i,j=1
xixjσ
ij , (46)
Zw = G0 +
Mn. (47)
Notice that if the series (45) is truncated after a given order n ≥ 3, Zw is given
by the first n moments of the size distribution only. On the other hand, Z still
involves an infinite number of moments if the truncation is made after n ≥ 4
due to the presence of terms like
i,j xixjσ
4/σij ,
i,j xixjσ
5/σ2ij , . . . .
Therefore, if we want the consistency condition Z = Zw to be satisfied for any
discrete or continuous polydisperse mixture, either the whole infinite series
(45) needs to be considered or it must be truncated after n = 3. The latter is
of course the simplest possibility and thus we make the approximation
G(η, z) = G0(η) + G1(η)z + G2(η)z2 + G3(η)z3. (48)
As a consequence, Z and Zw depend functionally on the size distribution of
the mixture only through the first three moments (which is in the spirit of
Rosenfeld’s Fundamental Measure Theory [48]).
Using the approximation (48) in Eqs. (46) and (47) we are led to
Z = G0 + η
(3G0 + 2G1) + 2
(G1 + 2G2 + 2G3)
, (49)
Alternative Approaches to Hard-Sphere Liquids 15
Zw = G0 + 2
G1 + 4
(G2 + 2G3) . (50)
Thus far, the dependence of both Z and Zw on the momentsM1,M2, andM3
is explicit and we only lack the packing-fraction dependence of G1, G2, and
G3. From Eqs. (49) and (50) it follows that the difference between Z and Zw
is given by
Z−Zw =
[3ηG0 − 2(1− η)G1]+2
[ηG1 − 2(1− η)G2 − 2(2− η)G3] .
Therefore, Z = Zw for any dispersity provided that
G1(η) =
2 (1− η)2
, (52)
G2(η) =
4 (1− η)3
− 2− η
G3(η), (53)
where use has been made of the definition of G0, Eq. (17). To close the problem,
we use the equal size limit given in Eq. (18), which yields G0+G1+G2+G3 = gs.
After a little algebra we are led to
G2(η) = (2− η)gs −
2 + η2/4
(1− η)2
, (54)
G3(η) = (1− η)
gSPTs − gs
. (55)
This completes the derivation of our improved approximation, which we will
call “e3”, following the same criterion as the one used to call “e1” and “e2” to
the approximations (20) and (34), respectively. In Eq. (55), gSPTs is the SPT
contact value for a single fluid, whose expression appears in Table 1. From
Eq. (55) it is obvious that the choice gs = g
s makes our e3 approximation
to become the e2 approximation, both reducing to the SPT for mixtures, Eq.
(8). This means that the SPT is fully internally consistent with the require-
ment Z = Zw, although it has the shortcoming of not being too accurate in
the single component case. The e3 proposal, on the other hand, satisfies the
condition Z = Zw and has the flexibility of accommodating any desired gs.
For the sake of concreteness, let us write explicitly the contact values in
the e3 aproximation:
ge3ij (σij) =
2 (1− η)2
(2 − η)gs −
2 + η2/4
(1− η)2
+ (1− η)
gSPTs − gs
, (56)
16 M. López de Haro, S. B. Yuste and A. Santos
ge3wj =
1− η +
(1− η)2
σj + 4
(2 − η)gs −
2 + η2/4
(1− η)2
+8(1− η)
gSPTs − gs
. (57)
With the above results the compressibility factor may be finally written in
terms of Zs as
Ze3(η) =
(1− η)
(1− η)2
Zs(η)−
. (58)
A few comments are in order at this stage. First, from Eq. (49) we can
observe that, for the class of approximations (48), the compressibility factor
Z does not depend on the individual values of the coefficients G2 and G3,
but only on their sum. As a consequence, two different approximations of
the form (48) sharing the same density dependence of G1 and G2 + G3 also
share the same virial EOS. For instance, if one makes the choice gs = g
then ZePY3 = ZPY, even though g
ij (σij) 6= gPYij (σij). Furthermore, if one
makes the more accurate choice gs = g
s , then ZeCS3 = ZBMCSL, but again
geCS3ij (σij) 6= gBGHLLij (σij). The eCS3 contact values are
geCS3ij (σij) =
2 (1− η)2
η2(1 + η)
4(1− η)3
4(1− η)2
, (59)
geCS3wj =
(1− η)2
η2(1 + η)
(1− η)3
(1 − η)2
. (60)
In Figs. 3 and 4 we display the performance of the contact values as given
by Eqs. (59) and (60), respectively, by comparison with results of computer
simulations for both discrete and polydisperse mixtures. In both figures we
have also included the results that follow from the classical proposals as well
as those of the eCS1 and eCS2 approximations. It is clear that for the wall-
particle contact values the eCS3 approximation yields the best performance,
while for the particle-particle contact values both the eCS2 and eCS3 are of
comparable accuracy. A further feature to be pointed out is that the practical
collapse on a common curve of the simulation data in Figs. 3 and 4 provide a
posteriori support for the universality ansatz made in Eq. (16).
As mentioned earlier, there exist extra consistency conditions (see for in-
stance Ref. [12]) that one might use as well within our approach. Assuming
Alternative Approaches to Hard-Sphere Liquids 17
0.0 0.5 1.0 1.5 2.0
 eCS1
 eCS2 
 eCS3 
Fig. 3. Plot of the difference gij(σij)− gBGHLLij (σij) as a function of the parameter
zij = (σiσj/σij)M2/M3 for hard spheres (d = 3) at a packing fraction η = 0.49.
The symbols are simulation data for the single fluid (circle, Ref. [36]), three binary
mixtures (squares, Ref. [49]) with σ2/σ1 = 0.3 and x1 = 0.0625, 0.125, and 0.25, and
a ternary mixture (triangles, Ref. [50]) with σ2/σ1 =
, σ3/σ1 =
, and x1 = 0.1,
x2 = 0.2. The lines are the PY approximation (– · · –), the SPT approximation
(– · – ·), the eCS1 approximation (· · · ), the eCS2 approximation (– – –), and the
eCS3 approximation (—).
that the ansatz (16) still holds, some of these conditions are related to the
derivatives of G with respect to z, namely
∂G(η, z)
2(1− η)2
, (61)
∂2G(η, z)
gPYs −
, (62)
∂3G(η, z)
= 0. (63)
Interestingly enough, as shown by Eq. (52), condition (61) is already satisfied
by our e3 approximation without having to be imposed. On the other hand,
condition (63) implies G3 = 0 in the e3 scheme and thus it is only satisfied if
gs = g
s , in which case we recover the SPT. Condition (62) is not fulfilled
either by the SPT or by the e3 approximation (except for a particular ex-
pression of gs which is otherwise not very accurate). Thus, fulfilling the extra
18 M. López de Haro, S. B. Yuste and A. Santos
0.0 0.5 1.0 1.5 2.0 2.5
 eCS1
 eCS2
 eCS3
Fig. 4. Plot of the difference gwj − gBGHLLwj as a function of the parameter zwj/2 =
σjM2/M3 for hard spheres (d = 3) at a packing fraction η = 0.4. The symbols are
simulation data for a polydisperse mixture with a narrow top-hat distribution (open
squares, Ref. [51]), a polydisperse mixture with a wide top-hat distribution (open
circles, Ref. [51]), a polydisperse mixture with a Schulz distribution (open triangles,
Ref. [51]), and a binary mixture (closed circles, Ref. [52]). The lines are the PY
approximation (– · · –), the SPT approximation (– · – ·), the eCS1 approximation
(· · · ), the eCS2 approximation (– – –), and the eCS3 approximation (—).
conditions (62) and (63) with a free gs requires either considering a higher
order polynomial in z (in which case the consistency condition Z = Zw can-
not be satisfied for arbitrary mixtures, as discussed before) or not using the
universality ansatz at all. In the first case, we have checked that a quartic or
even a quintic polynomial does not improve matters, whereas giving up the
universality assumption increases significantly the number of parameters to
be determined and seems not to be adequate in view of the behavior observed
in the simulation data.
An additional comment has to do with the restriction to d = 3 in this
subsection. As noted before, the approximation e1 reduces to the exact result
(3) for d = 1. For d = 2, the approximation e2 already fulfills the condition
Z = Zw and so there is no real need to go further in that case. Since we
have needed the approximation e3 to satisfy Z = Zw for d = 3, it is tempting
to speculate that a polynomial form for G(z) of degree d could be found to
be consistent with the condition Z = Zw for d ≥ 4. However, a detailed
analysis shows that this is not the case for an arbitrary mixture, since the
Alternative Approaches to Hard-Sphere Liquids 19
number of conditions exceeds the number of unknowns, unless the universality
assumption is partially relaxed.
As a final comment, let us stress that, although the discussion in this
section has referred, for the sake of simplicity, to discrete mixtures, all the
dependence on the details of the composition occurs through a finite number
of moments, so that the results remain meaningful even for continuous poly-
disperse mixtures [53]. In that case, instead of a set of mole fractions {xi} and
a set of diameters {σi}, one has to deal with a distribution function w(σ) such
that w(σ)dσ is the fraction of particles with a diameter comprised between σ
and σ + dσ. Therefore, the moments (1) are now defined as
dσ σnw(σ), (64)
and with such a change the results we have derived for discrete mixtures also
hold for polydisperse systems.
2.3 Non-Additive Systems
Non-additive hard-core mixtures, where the distance of closest approach be-
tween particles of different species is no longer the arithmetic mean of the
diameters of both particles, have received much less attention than additive
mixtures, in spite of their in principle more versatility to deal with interesting
aspects occurring in real systems (such as fluid-fluid phase separation) and
of their potential use as reference systems in perturbation calculations on the
thermodynamic and structural properties of, say, Lennard–Jones mixtures.
Nevertheless, the study of non-additive systems goes back fifty years [54–56]
and is still a rapidly developing and challenging problem.
As mentioned in the paper by Ballone et al. [57], where the relevant
references may be found, experimental work on alloys, aqueous electrolyte
solutions, and molten salts suggests that hetero-coordination and homo-
coordination may be interpreted in terms of excluded volume effects due to
non-additivity of the repulsive part of the intermolecular potential. In particu-
lar, positive non-additivity leads naturally to demixing in HS mixtures, so that
some of the experimental findings of phase separation in the above mentioned
(real) systems may be accounted for by using a model of a binary mixture of
(positive) non-additive HS. On the other hand, negative non-additivity seems
to account well for chemical short-range order in amorphous and liquid binary
mixtures with preferred hetero-coordination [58].
Some Preliminary Definitions
Let us consider an N -component mixture of non-additive HS in d dimensions.
In this case, σij =
(σi + σj)(1 + ∆ij), where ∆ij ≥ −1 is a symmetric
matrix with zero diagonal elements (∆ii = 0) that characterizes the degree
20 M. López de Haro, S. B. Yuste and A. Santos
of non-additivity of the interactions. If ∆ij > 0 the non-additivity character
of the ij interaction is said to be positive, while it is negative if ∆ij < 0. In
the case of a binary mixture (N = 2), the only non-additivity parameter is
∆ ≡ ∆12 = ∆21. The virial EOS (2) remains being valid in the non-additive
case.
The contact values gij(σij) can be expanded in a power series in density
gij(σij) = 1 + vdρ
xkck;ij + (vdρ)
k,ℓ=1
xkxℓckℓ;ij +O(ρ3). (65)
The coefficients ck;ij , ckℓ;ij , . . . are independent of the composition of the mix-
ture, but they are in general complicated nonlinear functions of the diameters
σij , σik, σjk, σkℓ, . . . . Insertion of the expansion (65) into Eq. (2) yields the
virial expansion of Z, namely
Z(ρ) = 1 +
Bn(vdρ)
= 1 + vdρ
i,j=1
Bijxixj + (vdρ)
i,j,k=1
Bijkxixjxk
+(vdρ)
i,j,k,ℓ=1
Bijkℓxixjxkxℓ +O(ρ4). (66)
Note that, for further convenience, we have introduced the coefficients Bn ≡
−(n−1)
d Bn, where Bn are the usual virial coefficients [cf. Eq. (25)]. The
composition-independent second, third, and fourth (barred) virial coefficients
are given by
Bij = 2
d−1σdij , (67)
Bijk =
ck;ijσ
ij + cj;ikσ
ik + ci;jkσ
, (68)
Bijkℓ =
ckℓ;ijσ
ij + cjℓ;ikσ
ik + ciℓ;jkσ
jk + cjk,iℓσ
iℓ + cik,jℓσ
+cij;kℓσ
. (69)
A Simple Proposal for the Equation of State of d-Dimensional
Non-Additive Mixtures
Our goal now is to generalize the e1 proposal given by Eq. (20) to the non-
additive case [59]. We will not try to extend the e2 and e3 proposals, Eqs. (34)
and (56), because of two reasons. First, given the inherent complexity of non-
additive systems, we want to keep the approach as simple as possible. Second,
Alternative Approaches to Hard-Sphere Liquids 21
we are more interested in the EOS than in the contact values themselves and,
as mentioned earlier, the e1 proposal provides excellent EOS, at least in the
additive case, despite the simplicity of the corresponding contact values.
As the simplest possible extension, we impose again the point particle and
equal size consistency conditions, Eqs. (12) and (13), and thus keep in this
case also the ansatz (16) and the linear structure of Eq. (19). However, instead
of using Eq. (15), we determine the parameters zij as to reproduce Eq. (65)
to first order in the density. The result is readily found to be [59]
zij =
)−1(∑
k xkck;ij
. (70)
Here b2 = 2
d−1 and b3 are the second and third virial coefficients for the
single component fluid, as defined by Eq. (26). The proposal of Eq. (19) sup-
plemented by Eq. (70) is, by construction, accurate for densities low enough
as to justify the truncated approximation gij(σij) ≈ 1 + vdρ
k xkck;ij . On
the other hand, the limitations of this truncated expansion for moderate and
large densities may be compensated by the use of gs. When Eqs. (16), (19),
and (70) are inserted into Eq. (2) one gets
Z(η) = 1 +
b3MdB2 − b2B3
(b3 − b2)M2d
+ [Zs(η) − 1]
B3 −MdB2
(b3 − b2)M2d
. (71)
Equation (71) is the sought generalization of Eq. (21) to non-additive hard-
core systems. As in the additive case, the the density dependence in the EOS
of the mixture is rather simple: Z(η)− 1 is expressed as a linear combination
of η/(1 − η) and Zs(η) − 1, with coefficients such that the second and third
virial coefficients are reproduced. Again, Eq. (71) is bound to be accurate for
sufficiently low densities, while the limitations of the truncated expansion for
moderate and large densities are compensated by the use of the EOS of the
pure fluid.
The exact second virial coefficient B2 is known from Eq. (67). In principle,
one should use the exact coefficients ck;ij to compute B3. However, to the
best of our knowledge they are only known for d ≤ 3. Since our objective is to
have a proposal which is explicit for any d, we can make use of a reasonable
approximation for them [59], as described below.
An Approximate Proposal for ck;ij
The values of the coefficients ck;ij are exactly known for d = 1 and d = 3 and
from these results one may approximate them in d dimensions as [59]
ck;ij = σ
k;ij +
σd−1k;ij
σi;jkσj;ik, (72)
where we have called
22 M. López de Haro, S. B. Yuste and A. Santos
σk;ij ≡ σik + σjk − σij (73)
and it is understood that σk;ij ≥ 0 for all sets ijk. Clearly, σi;ij = σi. For a
binary mixture Eq. (72) yields
c1;11 = (b3/b2)σ
c2;11 = (2σ12 − σ1)d + (b3/b2 − 1)σ1(2σ12 − σ1)d−1,
c1;12 = σ
1 + (b3/b2 − 1) (2σ12 − σ1)σd1/σ12.
Of course, Eqs. (72) and (74) reduce to the exact results for d = 1 (b2 = b3 = 1)
and for d = 3 (b2 = 4, b3 = 10).
The quantities σk;ij may be given a simple geometrical interpretation.
Assume that we have three spheres of species i, j, and k aligned in the sequence
ikj. In such a case, the distance of closest approach between the centers of
spheres i and j is σik + σjk. If the sphere of species k were not there, that
distance would of course be σij . Therefore σk;ij as given by Eq. (73) represents
a kind of effective diameter of sphere k, as seen from the point of view of the
interaction between spheres i and j.
Inserting Eq. (72) into Eq. (70), one gets
zij =
)−1(∑
k xkσ
k xkσ
k;ij σi;jkσj;ik
Mdσij
. (75)
It can be easily checked that in the additive case (σk;ij → σk), Eq. (75) reduces
to Eq. (15).
Equations (72) and (74) are restricted to the situation σk;ij ≥ 0 for any
choice of i, j, and k, i.e., 2σ12 ≥ max(σ1, σ2) in the binary case. This excludes
the possibility of dealing with mixtures with extremely high negative non-
additivity in which one sphere of species k might “fit in” between two spheres
of species i and j in contact. Since for d = 3 and N = 2 the coefficients ck;ij
are also known for such mixtures [60], we may extend our proposal to deal
with these cases:
c1;11 = (b3/b2)σ
c2;11 = σ̂
2 + (b3/b2 − 1)σ1σ̂d−12 ,
c1;12 = (2σ12 − σ̂2)d + (b3/b2 − 1) σ̂2σd1/σ12,
where we have defined
σ̂2 = max (2σ12 − σ1, 0) . (77)
With such an extension, we recover the exact values of ck;ij for a binary
mixture of hard spheres (d = 3), even if σ1 > 2σ12 or σ2 > 2σ12.
The EOS (71) becomes explicit when B3 is obtained from Eq. (68) by
using the approximation (72). The resulting virial coefficient is the exact one
for d = 1 and d = 3. For hard disks (d = 2), it turns out that the approximate
third virial coefficient is practically indistinguishable from the exact one [59].
Alternative Approaches to Hard-Sphere Liquids 23
-0.6 -0.4 -0.2 0.0 0.2 0.4 0.6
Fig. 5. Plot of the compressibility factor versus the non-additivity parameter ∆ for
a symmetric binary mixture of non-additive hard spheres (d = 3) at η = π/30 and
two different compositions. The solid lines are our proposal, Eq. (71), with Zs = Z
while the dashed lines are Hamad’s proposal (Refs. [61–63]). The symbols are results
from Monte Carlo simulations (Refs. [64,65]).
When the approximate B3 is used, Eq. (71) reduces to Eq. (21) in the additive
case.
From the comparison with simulation results, both for the compressibility
factor and higher order virial coefficients, we find that the EOS (71) does a
good job for non-additive mixtures, thus representing a reasonable compro-
mise between simplicity and accuracy, provided that Zs is accurate enough.
This is illustrated in Fig. 5, where the proposal (71) with Zs = Z
s and a
similar proposal by Hamad [61–63] are compared with simulation data [64,65]
for some three-dimensional symmetric mixtures. A more extensive compari-
son [59] shows that Eq. (71) seems to work better (especially as the density
is increased) in the case of positive non-additivities, at least for d = 1, d = 2,
and d = 3, but its performance is also reasonably good in highly asymmetric
mixtures, even for negative ∆. Of course the full assessment of this proposal
is still pending since it involves many facets (non-additivity parameters, size
ratios, density, and composition). Without this full assessment and given its
rather satisfactory performance so far, going beyond the approximation given
by Eq. (19) (taking similar steps to the ones described in Subsections 2.1 and
2.2 for additive systems) does not seem to be necessary at this stage, although
it is in principle feasible.
24 M. López de Haro, S. B. Yuste and A. Santos
2.4 Demixing
Demixing is a common phase transition in fluid mixtures usually originated
on the asymmetry of the interactions (e.g., their strength and/or range) be-
tween the different components in the mixture. In the case of athermal systems
such as HS mixtures in d dimensions, if fluid-fluid separation occurs, it would
represent a neat example of an entropy-driven phase transition, i.e., a phase
separation based only on the size asymmetry of the components. The exis-
tence of demixing in binary additive three dimensional HS mixtures has been
studied theoretically since decades, and the issue is still controversial. In this
subsection we will present our results following different but related routes
that attempt to clarify some aspects of this problem.
Binary Mixtures of Additive d-Dimensional Spheres (d = 3, d = 4
and d = 5)
Now we look at the possible instability of a binary fluid mixture of HS of
diameters σ1 and σ2 (σ1 > σ2) in d dimensions by looking at the Helmholtz
free energy per unit volume, f , which is given by
= −1 +
xi ln
Z(η′)− 1
, (78)
where λi is the thermal de Broglie wavelength of species i. We locate the
spinodals through the condition f11f22−f212 = 0, with fij ≡ ∂2f/∂ρi∂ρj . Due
to the spinodal instability, the mixture separates into two phases of different
composition. The coexistence conditions are determined through the equality
of the pressure p and the two chemical potentials µ1 and µ2 in both phases
(µi = ∂f/∂ρi), leading to binodal (or coexistence) curves.
We begin with the case d = 3. It is well known that the BMCSL EOS, Eq.
(11), does not lead to demixing. However, other EOS for HS mixtures have
been shown to predict demixing [41, 66], including the EOS that is obtained
by truncating the virial series after a certain number of terms [67, 68]. In
particular, it turns out that both Z = ZeCS1, Eq. (40), and Z = ZeCS2,
Eq. (41), lead to demixing for certain values of the parameter γ ≡ σ2/σ1
that measures the size asymmetry. The critical values of the pressure, the
composition, and the packing fraction are presented in Table 2 for a few
values of γ.
As discussed earlier, the eCS1 EOS and, to a lesser extent, the eCS2 EOS
are both in reasonably good agreement with the available simulation results
for the compressibility factor [18, 36, 45] and lead to the exact second and
third virial coefficients but differ in the predictions for Bn with n ≥ 4. The
scatter in the values for the critical constants shown in Table 2 is evident and
so there is no indication as to whether one should prefer one equation over
the other in connection with this problem. Notice, for instance, that the eCS2
Alternative Approaches to Hard-Sphere Liquids 25
Table 2. Critical constants pcσ
1/kBT , x1c, and ηc for different γ-values as obtained
from the two extended CS equations (40) and (41).
eCS1 eCS2
γ pcσ
1/kBT x1c ηc pcσ
1/kBT x1c ηc
0.05 3599 0.0093 0.822 1096 0.0004 0.204
0.1 1307 0.0203 0.757 832.0 0.0008 0.290
0.2 653.4 0.0537 0.725 — — —
0.3 581.9 0.0998 0.738 — — —
0.4 663.4 0.1532 0.766 — — —
does not predict demixing for γ ≥ 0.2, while both the values of the critical
pressures and packing fractions for which it occurs according to the eCS1 EOS
suggest that the transition might be metastable with respect to a fluid-solid
transition.
Now we turn to the cases d = 4 and d = 5. Here we use the extended
Luban–Michels equation (eLM1) described in Subsection 2.1 [see Eq. (21)
and Table 1]. As seen in Fig. 6, the location of the critical point tends to go
down and to the right in the η2 vs η1 plane as γ decreases for d = 4 [69].
On the other hand, while it also tends to go down as γ decreases if d = 5,
its behavior in the η2 vs η1 plane is rather more erratic in this case. Also,
the value of the critical pressure pc (in units of kBT/σ
1) is not a monotonic
function of γ; its minimum value lies between γ = 1/3 and γ = 1/2 when
d = 4, and it is around γ = 3/5 for d = 5. This non-monotonic behavior is
also observed for three-dimensional HS [66, 68].
It is conceivable that the demixing transition in binary mixtures of hard
hyperspheres in four and five dimensions described above may be metastable
with respect to a fluid-solid transition, as it may also be the case of 3D HS.
In fact, the value of the pressure at the freezing transition for the single
component fluid is [31] pfσ
d/kBT ≃ 12.7 (d = 3), 11.5 (d = 4), and 12.2
(d = 5), i.e., pfσ
d/kBT does not change appreciably with the dimensionality
but is clearly very small in comparison with the critical pressures pcσ
1/kBT
we obtain for the mixture; for instance, pcσ
1/kBT ≃ 600 (d = 3, γ = 3/10),
300 (d = 4, γ = 1/3) and 123 (d = 5, γ = 3/5). However, one should also
bear in mind that, if the concentration x1 of the bigger spheres decreases, the
value of the pressure at which the solid-fluid transition in the mixture occurs
in 3D is also considerably increased with respect to pf [cf. Fig. 6 of Ref. [66]].
Thus, for concentrations x1 ≃ 0.01 corresponding to the critical point of the
fluid-fluid transition, the maximum pressure of the fluid phase greatly exceeds
pf. If a similar trend with composition also holds in 4D and 5D, and given that
the critical pressures become smaller as the dimensionality d is increased, it is
not clear whether the competition between the fluid-solid and the fluid-fluid
transitions in these dimensionalities will always be won by the former. The
point clearly deserves further investigation.
26 M. López de Haro, S. B. Yuste and A. Santos
0.0 0.1 0.2 0.3 0.4 0.5 0.6
1400  2/ 1=3/4
4 1/k
0.0 0.1 0.2 0.3 0.4
5 1/k
Fig. 6. Spinodal curves (upper panels: lines) and binodal curves (upper panels: open
symbols; lower panels: lines) in a 4D system (left panels) and in a 5D system (right
panels). The closed symbols are the critical consolute points.
An interesting feature must be mentioned. There is a remarkable similarity
between the binodal curves represented in the pσdi –η1 and in the µi–η1 planes
[69]. By eliminating η1 as if it were a parameter, one can represent the binodal
curves in a µi vs pσ
i plane. Provided the origin of the chemical potentials
is such as to make λi = σi, the binodals in the µi–pσ
i plane practically
collapse into a single curve (which is in fact almost a straight line) for each
dimensionality (d = 3, d = 4, and d = 5) [69]. A closer analysis of this
Alternative Approaches to Hard-Sphere Liquids 27
phenomenon shows, however, that it is mainly due to the influence on µi
of terms which are quantitatively dominant but otherwise irrelevant to the
coexistence conditions.
Binary Mixtures of Non-Additive Hard Hyperspheres in the Limit
of High Dimensionality
Let us now consider a binary mixture of non-additive HS of diameters σ1
and σ2 in d dimensions. Thus in this case σ12 ≡ 12 (σ1 + σ2)(1 +∆) where as
before ∆ may be either positive or negative. Further assume (something that
will become exact in the limit d → ∞ [70]) that the EOS of the mixture is
described by the second virial coefficient only, namely
p = ρkBT [1 +B2(x1)ρ] , (79)
where, according to Eq. (67),
B2(x1) = vd2
1 + x
2 + 2x1x2σ
. (80)
The Helmholtz free energy per unit volume is given by f/ρkBT = −1 +∑2
i=1 xi ln
+B2ρ, where Eq. (78) has been used. The Gibbs free energy
per particle is
g = (f + p)/ρ =
xi ln
+ 2B2(x1)ρ, (81)
where without loss of generality we have set kBT = 1. Given a size ratio γ, a
value of ∆, and a dimensionality d, the consolute critical point (x1c, pc) is the
solution to
∂2g/∂x21
∂3g/∂x31
= 0, provided of course it exists. Then,
one can get the critical density ρc from Eq. (79).
We now introduce the scaled quantities [71]
p̃ ≡ 2d−1vdd−2pσd1/kBT, u ≡ d−1B2ρ. (82)
Consequently, Eqs. (79) and (81) can be rewritten as
p̃ = u
u+ d−1
/B̃2, (83)
xi ln (xiΛi) + ln
Adu/B̃2
+ 2du, (84)
where B̃2 ≡ B2/2d−1vdσd1 , Λi ≡ (λi/σ1)d, and Ad ≡ d/2d−1vd. Next we take
the limit d → ∞ and assume that the volume ratio γ̃ ≡ γd is kept fixed and
that there is a (slight) non-additivity ∆ = d−2∆̃ such that the scaled non-
additivity parameter ∆̃ is also kept fixed in this limit. Thus, the second virial
coefficient can be approximated by
28 M. López de Haro, S. B. Yuste and A. Santos
B̃2 = B̃
2 +B̃
−1+O(d−2), B̃(0)2 =
x1 + x2γ̃
2 = x1x2γ̃
1/2J,
J ≡ 1
(ln γ̃)
+ 2∆̃. (86)
Let us remark that, in order to find a consolute critical point, it is essential
to keep the term of order d−1 if ∆̃ ≤ 0. The EOS (83) can then be inverted
to yield
u = u(0)+u(1)d−1+O(d−2), u(0) =
2 , u
(1) = −1
1− u(0) B̃
In turn, the Gibbs free energy (84) becomes
g = g(0)d+ g(1) +O(d−1),
g(0) = 2u(0), g(1) =
i=1 xi ln (xiΛi) + ln
(0)/B̃
+ 2u(1),
while the chemical potentials µ1 = g+x2 (∂g/∂x1)p and µ2 = g−x1 (∂g/∂x1)p
are given by
µi = µ
i d+ µ
i +O(d−1), µ
1 = 2p̃
1 = ln
Adx1Λ1
p̃/B̃
2 + (x2/x1)(γ̃p̃)
1/2B̃
2 /B̃
where µ2 is obtained from µ1 by the changes x1 ↔ x2, Λ1 → Λ2/γ̃, γ̃ → 1/γ̃,
p̃→ p̃γ̃, B̃2 → B̃2/γ̃.
The coordinates of the critical point are readily found to be
x1c =
γ̃3/4
1 + γ̃3/4
, p̃c =
1 + γ̃1/4
4γ̃J2
. (90)
Note that x1c is independent of ∆̃. The coexistence curve, which has to be
obtained numerically, follows from the conditions µ
i (xA, p̃) = µ
i (xB , p̃)
(i = 1, 2) where x1 = xA and x1 = xB are the mole fractions of the co-
existing phases. Once the critical consolute point has been identified in the
pressure/concentration plane, we can obtain the critical density. The domi-
nant behaviors of B̃2 and u at the critical point are
2 (x1c) =
1− γ̃1/4 + γ̃1/2
)2 , u
1 + γ̃1/4
1− γ̃1/4 + γ̃1/2
. (91)
Hence, the critical density readily follows after substitution in the scaling
relation given in Eq. (82). It is also convenient to consider the scaled version
η̃ ≡ d−12dη of the packing fraction η = vdρσd1 (x1 + x2γ̃). At the critical point,
it takes the nice expression
Alternative Approaches to Hard-Sphere Liquids 29
0.00 0.04 0.08 0.12 0.16
p  =-0.1
 =0.1
=0.01
Fig. 7. Binodal curves in the planes ep vs x1 and eη vs x1 corresponding to eγ = 0.01
and e∆ = −0.1, e∆ = 0, and e∆ = 0.1.
η̃c =
γ̃1/8 + γ̃−1/8
. (92)
The previous results clearly indicate that a demixing transition is possible
not only for additive or positively non-additive mixtures but even for negative
non-additivities. The only requirement is J > 0, i.e., ∆̃ > − 1
(ln γ̃)
equivalently, ∆ > − 1
(ln γ)
. Figure 7 shows the binodal curves corresponding
to γ̃ = 0.01 and ∆̃ = −0.1 (negative non-additivity), ∆̃ = 0 (additivity), and
∆̃ = 0.1 (positive non-additivity).
While the high dimensionality limit has allowed us to address the prob-
lem in a mathematically simple and clear-cut way, the possibility of demixing
with negative non-additivity is not an artifact of that limit. As said before,
demixing is known to occur for positive non-additive binary mixtures of HS
in three dimensions and there is compelling evidence on the existence of this
phenomenon in the additive case, at least in the metastable fluid region. Even
though in a three-dimensional mixture the EOS is certainly more complicated
than Eq. (79) and the demixing transition that we have just discussed for neg-
ative non-additivity is possibly metastable with respect to the freezing transi-
tion, the main effects at work (namely the competition between depletion due
to size asymmetry and hetero-coordination due to negative non-additivity)
are also present. In fact, it is interesting to point out that Roth et al. [72],
using the approximation of an effective single component fluid with pair inter-
30 M. López de Haro, S. B. Yuste and A. Santos
actions to describe a binary mixture of non-additive 3D HS and employing an
empirical rule based on the effective second virial coefficient, have also sug-
gested that demixing is possible for small negative non-additivity and high
size asymmetry. Our exact results lend support to this suggestion and con-
firm that, in some cases, the limit d→ ∞ highlights features already present
in real systems.
3 The Rational Function Approximation (RFA) Method
for the Structure of Hard-Sphere Fluids
The RDF g(r) and its close relative the (static) structure factor S(q) are the
basic quantities used to discuss the structure of a single component fluid [1–4].
The latter quantity is defined as
S(q) = 1 + ρh̃(q), (93)
where
h̃(q) =
dr e−iq·rh(r) (94)
is the Fourier transform of the total correlation function h(r) ≡ g(r)−1, i being
the imaginary unit. An important related quantity is the direct correlation
function c(r), which is defined in Fourier space through the Ornstein–Zernike
(OZ) relation [1–4]
c̃(q) =
h̃(q)
1 + ρh̃(q)
, (95)
where c̃(q) is the Fourier transform of c(r)
The usual approach to obtain g(r) is through one of the integral equa-
tion theories, where the OZ equation is complemented by a closure relation
between c(r) and h(r) [1]. However, apart from requiring in general hard nu-
merical labor, a disappointing aspect is that the substitution of the (necessar-
ily) approximate values of g(r) obtained from them in the (exact) statistical
mechanical formulae may lead to the thermodynamic inconsistency problem.
The two basic routes to obtain the EOS of a single component fluid of HS
are the virial route, Eq. (14), and the compressibility route
χs ≡ kBT
= [1− ρc̃(0)]−1 = S(0)
= 1 + 2ddησ−d
dr rd−1h(r). (96)
Thermodynamic consistency implies that
χ−1s (η) =
[ηZs(η)], (97)
Alternative Approaches to Hard-Sphere Liquids 31
but, in general, this condition is not satisfied by an approximate RDF. In
the case of a HS mixture, the virial route is given by Eq. (2), while the
compressibility route is indicated below [cf. Eq. (145)].
In this section we describe the RFA method, which is an alternative to
the integral equation approach and in particular leads by construction to
thermodynamic consistency.
3.1 The Single Component HS Fluid
We begin with the case of a single component fluid of HS of diameter σ.
The following presentation is equivalent to the one given in Refs. [73, 74],
where all details can be found, but more suitable than the former for direct
generalization to the case of mixtures.
The starting point will be the Laplace transform
G(s) =
dr e−srrg(r) (98)
and the auxiliary function Ψ(s) defined through
G(s) =
[ρ+ esσΨ(s)]
. (99)
The choice of G(s) as the Laplace transform of rg(r) and the definition of
Ψ(s) from Eq. (99) are suggested by the exact form of g(r) to first order in
density [73].
Since g(r) = 0 for r < σ while g(σ+) = finite, one has
g(r) = Θ(r − σ)
g(σ+) + g′(σ+)(r − σ) + · · ·
, (100)
where g′(r) ≡ dg(r)/dr. This property imposes a constraint on the large s
behavior of G(s), namely
eσssG(s) = σg(σ+) +
g(σ+) + σg′(σ+)
s−1 +O(s−2). (101)
Therefore, lims→∞ e
sσsG(s) = σg(σ+) = finite or, equivalently,
s−2Ψ(s) =
2πσg(σ+)
= finite. (102)
On the other hand, according to Eq. (96) with d = 3,
χs = 1− 24ησ−3 lim
dr e−srr [g(r) − 1]
= 1− 24ησ−3 lim
G(s)− s−2
. (103)
Since the (reduced) isothermal compressibility χs is also finite, one has∫∞
dr r2 [g(r) − 1] = finite, so that the weaker condition
dr r [g(r)− 1] =
lims→0[G(s)− s−2] = finite must hold. This in turn implies
32 M. López de Haro, S. B. Yuste and A. Santos
Ψ(s) = −ρ+ ρσs−
ρσ2s2+
ρσ3 +
ρσ3 +
σs4+O(s5).
(104)
First-Order Approximation (PY Solution)
An interesting aspect to be remarked is that the minimal input we have just
described on the physical requirements related to the structure and thermo-
dynamics of the system is enough to determine the small and large s limits of
Ψ(s), Eqs. (102) and (104), respectively. While infinite choices for Ψ(s) would
comply with such limits, a particularly simple form is a rational function. In
particular, the rational function having the least number of coefficients to be
determined is
Ψ(s) =
E(0) + E(1)s+ E(2)s2 + E(3)s3
L(0) + L(1)s
, (105)
where one of the coefficients can be given an arbitrary non-zero value. We
choose E(3) = 1. With such a choice and in view of Eq. (104), one finds
E(0) = −ρL(0), E(1) = −ρ(L(1) − σL(0)), E(2) = ρ(σL(1) − 1
σ2L(0)), and
L(0) = 2π
1 + 2η
(1− η)2
, (106)
L(1) = 2πσ
1 + η/2
(1− η)2
. (107)
Upon substitution of these results into Eqs. (99) and (105), we get
G(s) =
L(0) + L(1)s
ϕ2(σs)σ3L(0) + ϕ1(σs)σ2L(1)
] , (108)
where
ϕn(x) ≡ x−(n+1)
(−x)m
− e−x
. (109)
In particular,
ϕ0(x) =
1− e−x
, ϕ1(x) =
1− x− e−x
, ϕ2(x) =
1− x+ x2/2− e−x
(110)
Note that limx→0 ϕn(x) = (−1)n/(n+ 1)!.
It is remarkable that Eq. (108), which has been derived here as the sim-
plest rational form for Ψ(s) complying with the requirements (102) and (104),
coincides with the solution to the PY closure, c(r) = 0 for r > σ, of the OZ
equation [29]. Application of Eq. (102) yields the PY contact value gPYs and
compressibility factor ZPYs shown in Table 1. Analogously, Eq. (103) yields
χPYs =
(1− η)4
(1 + 2η)2
. (111)
It can be easily checked that the thermodynamic relation (97) is not satisfied
by the PY theory.
Alternative Approaches to Hard-Sphere Liquids 33
Second-Order Approximation
In the spirit of the RFA, the simplest extension of the rational approximation
(105) involves two new terms, namely αs4 in the numerator and L(2)s2 in the
denominator, both of them necessary in order to satisfy Eq. (102). Such an
addition leads to
Ψ(s) =
E(0) + E(1)s+ E(2)s2 + E(3)s3 + αs4
L(0) + L(1)s+ L(2)s2
. (112)
Applying Eq. (104), it is possible to express E(0), E(1), E(2), E(3), L(0), and
L(1) in terms of α and L(2). This leads to
G(s) =
L(0) + L(1)s+ L(2)s2
1 + αs− ρ
ϕ2(σs)σ3L(0) + ϕ1(σs)σ2L(1) + ϕ0(σs)σL(2)
(113)
where
L(0) = 2π
1 + 2η
(1− η)2
, (114)
L(1) = 2πσ
1 + 1
(1 − η)2
1 + 2η
α− 3ηL
. (115)
Thus far, irrespective of the values of the coefficients L(2) and α, the condi-
tions lims→∞ e
sσsG(s) = finite and lims→0[G(s) − s−2] = finite are satisfied.
Of course, if L(2) = α = 0, one recovers the PY approximation. More gen-
erally, we may determine these coefficients by prescribing the compressibility
factor Zs (or equivalently the contact value gs) and then, in order to ensure
thermodynamic consistency, compute from it the isothermal compressibility
χs by means of Eq. (97). From Eqs. (102) and (103) one gets
L(2) = 2πασgs, (116)
1− 12η
1 + 2
αL(2)
. (117)
Clearly, upon substitution of Eqs. (114) and (116) into Eq. (117) a quadratic
algebraic equation for α is obtained. The physical root is
α = − 12η(1 + 2η)E4
(1− η)2 + 36η [1 + η − Zs(1− η)]E4
, (118)
where
Zs − 13
Zs − 13
Zs − ZPYs
)]1/2}
. (119)
The other root must be discarded because it corresponds to a negative value
of α, which, according to Eq. (116), yields a negative value of L(2). This would
34 M. López de Haro, S. B. Yuste and A. Santos
imply the existence of a positive real value of s at which G(s) = 0 [73, 74],
which is not compatible with a positive definite RDF. However, according to
the form of Eq. (119) it may well happen that, once Zs has been chosen, there
exists a certain packing fraction ηg above which α is no longer positive. This
may be interpreted as an indication that, at the packing fraction ηg where α
vanishes, the system ceases to be a fluid and a glass transition in the HS fluid
occurs [74–76].
Expanding (113) in powers of s and using Eq. (101) one can obtain the
derivatives of the RDF at r = σ+ [77]. In particular, the first derivative is
g′(σ+) =
L(1) − L(2)
, (120)
which may have some use in connection with perturbation theory [15].
It is worthwhile to point out that the structure implied by Eq. (113) coin-
cides in this single component case with the solution of the Generalized Mean
Spherical Approximation (GMSA) [78], where the OZ relation is solved under
the ansatz that the direct correlation function has a Yukawa form outside the
core.
For a given Zs, once G(s) has been determined, inverse Laplace trans-
formation yields rg(r). First, note that Eq. (99) can be formally rewritten
G(s) = −
ρn−1 [−Ψ(s)]−n e−nsσ. (121)
Thus, the RDF is then given by
g (r) =
ρn−1ψn (r − nσ)Θ (r − nσ) , (122)
with Θ (x) denoting the Heaviside step function and
ψn (r) = −L−1
s [−Ψ (s)]−n
, (123)
L−1 denoting the inverse Laplace transform. Explicitly, using the residue the-
orem,
ψn (r) = −
(n−m)!(m− 1)!
rn−m, (124)
where
a(i)mn = lim
s [−Ψ (s) /(s− si)]−n , (125)
si (i = 1, . . . , 4) being the poles of 1/Ψ(s), i.e., the roots of E
(0) + E(1)s +
E(2)s2 + E(3)s3 + αs4 = 0. Explicit expressions of g(r) up to the second
coordination shell σ ≤ r ≤ 3σ can be found in Ref. [79].
Alternative Approaches to Hard-Sphere Liquids 35
On the other hand, the static structure factor S(q) [cf. Eq. (93)] and the
Fourier transform h̃(q) may be related to G(s) by noting that
h̃(q) =
dr r sin(qr)h(r) = −2π G(s)−G(−s)
. (126)
Therefore, the basic structural quantities of the single component HS fluid,
namely the RDF and the static structure factor, may be analytically deter-
mined within the RFA method once the compressibility factor Zs, or equiva-
lently the contact value gs, is specified. In Fig. 8 we compare simulation data
of g(r) for a density ρσ3 = 0.9 [80] with the RFA prediction and a recent
approach by Trokhymchuk et al. [81], where Zs = Z
s [cf. Table 1] and the
associated compressibility
χCSs =
(1− η)4
1 + 4η + 4η2 − 4η3 + η4
(127)
are taken in both cases. Both theories are rather accurate, but the RFA cap-
tures better the maxima and minima of g(r) [82].
It is also possible to obtain within the RFA method the direct correlation
function c(r). Using Eqs. (95) and (126), and applying the residue theorem,
one gets, after some algebra,
1 2 3 4
1 2 3 4
Fig. 8. Radial distribution function of a single component HS fluid for ρσ3 = 0.9.
The solid lines represent simulation data [80]. The dashed lines represent the results
of the approach of Ref. [81], while the dotted lines refer to those of the RFA method.
The inset shows the oscillations of g(r) in more detail.
36 M. López de Haro, S. B. Yuste and A. Santos
c(r) =
+K0 +K1r +K3r
Θ(1 − r) +K
(128)
where
12αηL(2)/π + 1− 12α(1 + 2α)η/(1− η), (129)
4α2(1 − η)4κ6
2 [1 + 2(1 + 3α)η]± [2 + η + 2α(1 + 2η)]κ
+(1− η)
κ2 − η (12 + (κ± 6)κ)
L(2)/π
12η [1 + 2(1 + 3α)η]
±6η [3η − 2α(1− 4η)]κ− 6η(1 + 2α)(1− η)κ2 − (1 − η)2κ3(ακ∓ 1)
+6η(1− η)
κ2 − η (12 + (κ± 6)κ)
L(2)/π
, (130)
K−1 = −
κ +K−e
−κ +K0 +K1 +K3
, (131)
K0 = −
1 + 2 (1 + 3α) η − 6η (1− η)L(2)/π
ακ (1− η)2
, (132)
2α2κ2 (1− η)4
[2 + η + 2α(1 + 2η)]
2 − 4 (1− η) [1 + η
×(7 + η + 6α (2 + η))]L(2)/π + 12η (2 + η) (1− η)2L(2)
,(133)
K0, (134)
K = − (K+ +K− +K−1) . (135)
In Eqs. (129)–(135) we have taken σ = 1 as the length unit. Note that Eq.
(135) guarantees that c(0) = finite, while Eq. (131) yields c(σ+) − c(σ−) =
L(2)/2πα = g(σ+). The latter equation proves the continuity of the indirect
correlation function γ(r) ≡ h(r) − c(r) at r = σ. With the above results,
Eqs. (122) and (128), one may immediately write the function γ(r). Finally,
we note that the bridge function B(r) is linked to γ(r) and to the cavity
(or background) function y(r) ≡ eφ(r)/kBT g(r), where φ(r) is the interaction
potential, through
B(r) = ln y(r) − γ(r), (136)
and so, within the RFA method, the bridge function is also completely speci-
fied analytically for r > σ once Zs is prescribed.
If one wants to have B(r) also for 0 ≤ r ≤ σ, then an expression for the
cavity function is required in that region. Here we propose such an expression
using a limited number of constraints. First, since the cavity function and its
first derivative are continuous at r = σ, we have
Alternative Approaches to Hard-Sphere Liquids 37
0.0 0.2 0.4 0.6 0.8 1.0
1000 3=0.3
3=0.5
3=0.7
Fig. 9. Cavity function of a single component HS fluid in the overlap region for
ρσ3 = 0.3, 0.5, and 0.7. The solid lines represent our proposal (140) with Zs = Z
while the symbols represent Monte Carlo simulation results [84].
y(1) = gs,
y′(1)
− 1, (137)
where Eqs. (116) and (120) have been used and again σ = 1 has been taken.
Next, we consider the following exact zero-separation theorems [83]:
ln y(0) = Zs(η)− 1 +
′)− 1
, (138)
y′(0)
= −6ηy(1). (139)
The four conditions (137)–(139) can be enforced by assuming a cubic poly-
nomial form for ln y(r) inside the core, namely
y(r) = exp
Y0 + Y1r + Y2r
2 + Y3r
, (0 ≤ r ≤ 1), (140)
where
Y0 = Zs(η)− 1 +
′)− 1
, (141)
Y1 = −6ηy(1), (142)
38 M. López de Haro, S. B. Yuste and A. Santos
0 2 4 6 8 10
 RFA, =0.3
 RFA, =0.49
Fig. 10. Parametric plot of the bridge function B(r) versus the indirect correlation
function γ(r). The dashed line refers to the RFA for η = 0.3, while the solid line
refers to the RFA for η = 0.49. In each case, the branch of the curve to the right
of the circle corresponds to r ≤ 1, while that to the left corresponds to r ≥ 1. For
comparison, the PY closure B(r) = ln[1 + γ(r)]− γ(r) is also plotted (dash-dotted
line).
Y2 = 3 ln y(1)−
y′(1)
− 3Y0 − 2Y1, (143)
Y3 = −2 ln y(1) +
y′(1)
+ 2Y0 + Y1. (144)
The proposal (140) is compared with available Monte Carlo data [84] in Fig.
9, where an excellent agreement can be observed.
Once the cavity function y(r) provided by the RFA method is comple-
mented by (140), the bridge function B(r) can be obtained at any distance.
Figure 10 presents a parametric plot of the bridge function versus the indirect
correlation function as given by the RFA method for two different packing
fractions, as well as the result associated with the PY closure. The fact that
one gets a smooth curve means that within the RFA the oscillations in γ(r)
are highly correlated to those of B(r). Further, the effective closure relation in
the RFA turns out to be density dependent, in contrast with what occurs for
the PY theory. Note that the absolute value |B(r)| for a given value of γ(r)
is smaller in the RFA than the PY value and that the RFA and PY curves
become paradoxically closer for larger densities. Since the PY theory is known
Alternative Approaches to Hard-Sphere Liquids 39
to yield rather poor values of the cavity function inside the core [85, 86], it
seems likely that the present differences may represent yet another manifes-
tation of the superiority of the RFA method, a point that certainly deserves
to be further explored.
3.2 The Multicomponent HS Fluid
The method outlined in the preceding subsection will be now extended to an
N -component mixture of additive HS. Note that in a multicomponent system
the isothermal compressibility χ is given by
χ−1 =
T,{xj}
T,{xj}
= 1− ρ
i,j=1
xixj c̃ij(0), (145)
where c̃ij(q) is the Fourier transform of the direct correlation function cij(r),
which is defined by the OZ equation
h̃ij(q) = c̃ij(q) +
ρkh̃ik(q)c̃kj(q), (146)
where hij(r) ≡ gij(r) − 1. Equations (145) and (146) are the multicompo-
nent extensions of Eqs. (96) and (95), respectively. Introducing the quantities
ĥij(q) ≡
ρiρj h̃ij(q) and ĉij(q) ≡
ρiρj c̃ij(q), the OZ relation (146) be-
comes, in matrix notation,
ĉ(q) = ĥ(q) · [I+ ĥ(q)]−1, (147)
where I is the N ×N identity matrix. Thus, Eq. (145) can be rewritten as
χ−1 =
i,j=1
xixj [δij − ĉij(0)] =
i,j=1
I+ ĥ(0)
. (148)
Similarly to what we did in the single component case, we introduce the
Laplace transforms of rgij(r):
Gij(s) =
dr e−srrgij(r). (149)
The counterparts of Eqs. (100) and (101) are
gij(r) = Θ(r − σij)
gij(σ
ij) + g
ij)(r − σij) + · · ·
, (150)
40 M. López de Haro, S. B. Yuste and A. Santos
eσijssGij(s) = σijgij(σ
ij) +
gij(σ
ij) + σijg
s−1 +O(s−2). (151)
Moreover, the condition of a finite compressibility implies that h̃ij(0) = finite.
As a consequence, for small s,
s2Gij(s) = 1 +H
3 + · · · (152)
with H
ij = finite and H
ij = −h̃ij(0)/4π = finite, where
dr (−r)nrhij(r). (153)
We are now in the position to generalize the approximation (113) to the
N -component case [87]. While such a generalization may be approached in a
variety of ways, two motivations are apparent. On the one hand, we want to
recover the PY result [8] as a particular case in much the same fashion as in
the single component system. On the other hand, we want to maintain the
development as simple as possible. Taking all of this into account, we propose
Gij(s) =
e−σijs
L(s) · [(1 + αs)I− A(s)]−1
, (154)
where L(s) and A(s) are the matrices
Lij(s) = L
ij + L
ij s+ L
2, (155)
Aij(s) = ρi
ϕ2(σis)σ
ij + ϕ1(σis)σ
ij + ϕ0(σis)σiL
, (156)
the functions ϕn(x) being defined by Eq. (109). We note that, by construc-
tion, Eq. (154) complies with the requirement lims→∞ e
σijssGij(s) = finite.
Further, in view of Eq. (152), the coefficients of s0 and s in the power series
expansion of s2Gij(s) must be 1 and 0, respectively. This yields 2N
2 condi-
tions that allow us to express L(0) and L(1) in terms of L(2) and α. The solution
is [87]
ij = ϑ1 + ϑ2σj + 2ϑ2α− ϑ1
ρkσkL
kj , (157)
ij = ϑ1σij +
ϑ2σiσj + (ϑ1 + ϑ2σi)α−
ρkσkL
kj , (158)
where ϑ1 ≡ 2π/(1− η) and ϑ2 ≡ 6π(M2/M3)η/(1− η)2.
In parallel with the development of the single component case, L(2) and α
can be chosen arbitrarily. Again, the choice L
ij = α = 0 gives the PY solution
[8, 88]. Since we want to go beyond this approximation, we will determine
those coefficients by taking prescribed values for gij(σij), which in turn, via
Eq. (2), give the EOS of the mixture. This also leads to the required value of
Alternative Approaches to Hard-Sphere Liquids 41
χ−1 = ∂(ρZ)/∂ρ, thus making the theory thermodynamically consistent. In
particular, according to Eq. (151),
ij = 2πασijgij(σ
ij). (159)
The condition related to χ is more involved. Making use of Eq. (152), one can
get h̃ij(0) = −4πH(1)ij in terms of L(2) and α and then insert it into Eq. (148).
Finally, elimination of L
ij in favor of α from Eq. (159) produces an algebraic
equation of degree 2N , whose physical root is determined by the requirement
that Gij(s) is positive definite for positive real s. It turns out that the physical
solution corresponds to the smallest of the real roots. Once α is known, upon
substitution into Eqs. (154), (157), (158), and (159), the scheme is complete.
Also, using Eq. (151), one can easily derive the result
g′ij(σ
ij) =
2πασij
ij − L
. (160)
It is straightforward to check that the results of the preceding subsection are
recovered by setting σi = σ, regardless of the values of the mole factions.
Once Gij(s) has been determined, inverse Laplace transformation directly
yields rgij(r). Although in principle this can be done analytically, it is more
practical to use one of the efficient methods discussed by Abate and Whitt [89]
to numerically invert Laplace transforms [90].
In Fig. 11 we present a comparison between the results of the RFA method
with the PY theory and simulation data [50] for the RDF of a ternary mixture.
In the case of the RFA, we have used the eCS2 contact values and the cor-
responding isothermal compressibility. The improvement of the RFA over the
PY prediction, particularly in the region near contact, is noticeable. Although
the RFA accounts nicely for the observed oscillations, it seems to somewhat
overestimate the depth of the first minimum.
Explicit knowledge of Gij(s) also allows us to determine the Fourier trans-
form h̃ij(q) through the relation
h̃ij(q) = −2π
Gij(s)−Gij(−s)
. (161)
The structure factor Sij(q) may be expressed in terms of h̃ij(q) as [4]
Sij(q) = xiδij + ρxixj h̃ij(q). (162)
In the particular case of a binary mixture, rather than the individual structure
factors Sij(q), it is some combination of them which may be easily associated
with fluctuations of the thermodynamic variables [91, 92]. Specifically, the
quantities [4]
Snn(q) = S11(q) + S22(q) + 2S12(q), (163)
42 M. López de Haro, S. B. Yuste and A. Santos
2 3 4 5 6 7
2 3 4 5 6 7 8
3 4 5 6 7 8
10.750.50.250
10.750.50.250
1 2 3 4 5 6 7
2 3 4 5 6 7 8
3 4 5 6 7 8 9
Fig. 11. Radial distribution functions gij(r) for a ternary mixture with diameters
σ1 = 1, σ2 = 2, and σ3 = 3 at a packing fraction η = 0.49 with mole fractions
x1 = 0.7, x2 = 0.2, and x3 = 0.1. The circles are simulation results [50], the solid
lines are the RFA predictions, and the dotted lines are the PY predictions.
Snc(q) = x2S11(q)− x1S22(q) + (x2 − x1)S12(q), (164)
Scc(q) = x
2S11(q) + x
1S22(q)− 2x1x2S12(q) (165)
are sometimes required.
After replacement of ĥij(q) =
ρiρj h̃ij(q) in Eq. (147), one easily gets
c̃ij(q). Subsequent inverse Fourier transformation yields cij(r). The result
gives cij(r) for r > σij as the superposition of N Yukawas [93], namely
Alternative Approaches to Hard-Sphere Liquids 43
cij(r) =
e−κℓr
, (166)
where q = ±iκℓ with ℓ = 1, . . . , N are the zeros of det
I+ ĥ(q)
and the
amplitudes K
ij are obtained by applying the residue theorem as
q→iκℓ
c̃ij(q)(q − iκℓ). (167)
The indirect correlation functions γij(r) ≡ hij(r) − cij(r) readily follow
from the previous results for the RDF and direct correlation functions. Finally,
in this case the bridge functions Bij(r) for r > σij are linked to gij(r) and
cij(r) through
Bij(r) = ln gij(r) − γij(r) (168)
and so once more we have a full set of analytical results for the structural
properties of a multicomponent fluid mixture of HS once the contact values
gij(σij) are specified.
4 Other Related Systems
The philosophy behind the RFA method to derive the structural properties
of three-dimensional HS systems can be adapted to deal with other related
systems. The main common features of the RFA can be summarized as follows.
First, one chooses to represent the RDF in Laplace space. Next, using as a
guide the low-density form of the Laplace transform, an auxiliary function is
defined which is approximated by a rational or a rational-like form. Finally,
the coefficients are determined by imposing some basic consistency conditions.
In this section we consider the cases of sticky-hard-sphere, square-well, and
hard-disk fluids. In the two former cases the RFA program is followed quite
literally, while in the latter case it is done more indirectly through the RFA
method as applied to hard rods (d = 1) and hard spheres (d = 3).
4.1 Sticky Hard Spheres
The sticky-hard-sphere (SHS) fluid model has received a lot of attention since
it was first introduced by Baxter in 1968 [94] and later extended to multi-
component mixtures by Perram and Smith [95] and, independently, by Bar-
boy [96]. In this model, the molecular interaction may be defined via square-
well (SW) potentials of infinite depth and vanishing width, thus embodying
the two essential characteristics of real molecular interactions, namely a harsh
repulsion and an attractive part. In spite of their known shortcomings [97], an
important feature of SHS systems is that they allow for an exact solution of the
44 M. López de Haro, S. B. Yuste and A. Santos
OZ equation in the PY approximation [94,95]. Furthermore, they are thought
to be appropriate for describing structural properties of colloidal systems, mi-
celles, and microemulsions, as well as some aspects of gas-liquid equilibrium,
ionic fluids and mixtures, solvent mediated forces, adsorption phenomena,
polydisperse systems, and fluids containing chainlike molecules [98–102].
Let us consider an N -component mixture of spherical particles interacting
according to the SW potential
φij(r) =
∞, r < σij ,
−ǫij , σij < r < Rij ,
0, r > Rij .
(169)
As in the case of additive HS, σij = (σi + σj)/2 is the distance between the
centers of a sphere of species i and a sphere of species j at contact. In addition,
ǫij is the well depth and Rij − σij indicates the well width. We now take the
SHS limit [94], namely
Rij → σij , ǫij → ∞, τij ≡
Rij − σij
e−ǫij/kBT = finite, (170)
where the τij are monotonically increasing functions of the temperature T and
their inverses measure the degree of “adhesiveness” of the interacting spheres i
and j. Even without strictly taking the mathematical limits (170), short-range
SW fluids can be well described in practice by the SHS model [103].
The virial EOS for the SHS mixture is given by
Z = 1 +
i,j=1
dr ryij(r)
e−φij(r)/kBT
= 1 +
i,j=1
xixjσ
ijyij(σij)
12τij
y′ij(σij)
yij(σij)
, (171)
where yij(r) ≡ gij(r)eφij(r)/kBT is the cavity function and y′ij(r) = dyij(r)/dr.
Since yij(r) must be continuous, it follows that
gij(r) = yij(r)
Θ(r − σij) +
12τij
δ(r − σij)
. (172)
The case of a HS system is recovered by taking the limit of vanishing adhe-
siveness τ−1ij → 0, in which case Eq. (171) reduces to the three-dimensional
version of Eq. (2). On the other hand, the compressibility EOS, Eq. (145), is
valid for any interaction potential, including SHS.
As in the case of HS, it is convenient to define the Laplace transform (149).
The condition yij(σij) = finite translates into the following large s behavior
of Gij(s):
eσijsGij(s) = σ
ijyij(σij)
12τij
+ σ−1ij s
+O(s−2), (173)
Alternative Approaches to Hard-Sphere Liquids 45
which differs from (151): while eσijsGij(s) ∼ s−1 for HS, eσijsGij(s) ∼ s0 for
SHS. However, the small s behavior is still given by Eq. (152), as a consequence
of the condition χ−1 = finite.
The RFA proposal for SHS mixtures [104] keeps the form (154), except
that now
Lij(s) = L
ij + L
ij s+ L
2 + L
3, (174)
Aij(s) = ρi
ϕ2(σis)σ
ij + ϕ1(σis)σ
ij + ϕ0(σis)σiL
ij − e
−σisL
(175)
instead of Eqs. (155) and (156). By construction, Eqs. (154), (174), and (175)
comply with the requirement lims→∞ e
σijsGij(s) = finite. Further, in view
of Eq. (152), the coefficients of s0 and s in the power series expansion of
s2Gij(s) must be 1 and 0, respectively. This yields 2N
2 conditions that allow
us to express L(0) and L(1) in terms of L(2), L(3), and α as [104]
ij = ϑ1+ϑ2σj+2ϑ2α−ϑ1
kj − L
ρkσkL
kj , (176)
ij = ϑ1σij +
ϑ2σiσj + (ϑ1 + ϑ2σi)α−
kj − L
(ϑ1 + ϑ2σi)
ρkσkL
kj , (177)
where ϑ1 and ϑ2 are defined below Eq. (158). We have the freedom to choose
(3) and α, but L(2) is constrained by the condition (173), i.e., the ratio
between the first and second terms in the expansion of eσijsGij(s) for large s
must be exactly equal to σij/12τij.
First-Order Approximation (PY Solution)
The simplest approximation consists of making α = 0. In view of the condition
eσijsGij(s) ∼ s0 for large s, this implies L(3)ij = 0. In that case, the large s
behavior that follows from Eq. (154) is
2πeσijsGij(s) = L
(2) · D
s−1 +O(s−2), (178)
where
Dij ≡ ρi
σ2i L
ij − σiL
ij + L
. (179)
Comparison with Eq. (173) yields
yij(σij) =
πσ2ij
ij , (180)
46 M. López de Haro, S. B. Yuste and A. Santos
12τijL
ik Dkj . (181)
Taking into account Eqs. (176) and (177) (with L
ij = L
ji and of course also
with α = 0 and L(3) = 0), Eq. (181) becomes a closed equation for L(2):
12τijL
= ϑ1σij+
ϑ2σiσj−
ki σj + L
kj σi
(182)
The physical root L(2) of Eq. (182) is the one vanishing in the HS limit τij →
∞. Once known, Eq. (180) gives the contact values.
This first-order approximation obtained from the RFA method turns out
to coincide with the exact solution of the PY theory for SHS [95].
Second-Order Approximation
As in the case of HS mixtures, a more flexible proposal is obtained by keeping
α (and, consequently, L
ij ) different from zero. In that case, instead of Eq.
(178), one has
2πeσijsGij(s) =
+O(s−2). (183)
This implies
πσ2ij
αyij(σij), (184)
12τijL
. (185)
If we fix yij(σij), Eqs. (176), (177), (184), and (185) allow one to express L
(1), L(2), and L(3) as linear functions of α. Thus, only the scalar parameter
α remains to be fixed, analogously to what happens in the HS case. As done
in the latter case, one possibility is to choose α in order to reproduce the
isothermal compressibility χ given by Eq. (148). To do so, one needs to find
the coefficients H
ij appearing in Eq. (152). The result is [104]
(0) = C(0) ·
I− A(0)
, (186)
(1) = C(1) ·
I− A(0)
, (187)
where
Alternative Approaches to Hard-Sphere Liquids 47
αδkj −A(1)kj
δkj −A(0)kj
(188)
σ2ik +H
αδkj −A(1)kj
σ3ik + σikH
δkj −A(0)kj
, (189)
ij = (−1)
σn+3i
(n+ 3)!
σn+2i
(n+ 2)!
σn+1i
(n+ 1)!
(190)
Equation (187) gives H(1) in terms of α: H
ij = Pij(α)/[Q(α)]
2, where Pij(α)
denotes a polynomial in α of degree 2N and Q(α) denotes a polynomial of
degree N . It turns out then that, seen as a function of α, χ is the ratio of
two polynomials of degree 2N . Given a value of χ, one may solve for α. The
physical solution, which has to fulfill the requirement that Gij(s) is positive
definite for positive real s, corresponds to the smallest positive real root.
Once α is known, the scheme is complete: Eq. (184) gives L(3), then L(2) is
obtained from Eq. (185), and finally L(1) and L(0) are given by Eqs. (176) and
(177), respectively. Explicit knowledge of Gij(s) through Eqs. (154), (174),
and (175) allows one to determine the Fourier transform h̃ij(q) and the struc-
ture factor Sij(q) through Eqs. (161) and (162), respectively. Finally, inverse
Laplace transformation of Gij(s) yields gij(r) [90].
Single Component SHS Fluids
The special case of single component SHS fluids [105, 106] can be obtained
from the multicomponent one by taking σij = σ and τij = τ . Thus, the
Laplace transform of rg(r) in the RFA is
G(s) =
L(0) + L(1)s+ L(2)s2 + L(3)s3
1 + αs− ρ
ϕ2(s)L(0) + ϕ1(s)L(1) + ϕ0(s)L(2) − e−sL(3)
(191)
where we have taken σ = 1. Equations (176) and (177) become
L(0) = 2π
1 + 2η
(1− η)2
− L(2)
(1 − η)2
(1 − 4η)L(3), (192)
L(1) = 2π
1 + 1
(1− η)2
1 + 2η
α− 3ηL(2)
− 18η
(1− η)2
L(3). (193)
The choice α = L(3) = 0 makes Eq. (191) coincide with the exact solution
to the PY approximation for SHS [94], where L(2) is the physical root (i.e., the
one vanishing in the limit τ → ∞) of the quadratic equation [see Eq. (182)]
48 M. López de Haro, S. B. Yuste and A. Santos
12τL(2) = 2π
1 + 2η
(1− η)2
− 12η
L(2) +
ηL(2)
. (194)
We can go beyond the PY approximation by prescribing a contact value
y(1), so that, according to Eqs. (184) and (185),
L(3) =
y(1), (195)
L(2) =
12τ +
L(3). (196)
By prescribing the isothermal compressibility χ, the parameter α can be ob-
tained as the physical solution (namely, the one remaining finite in the limit
τ → ∞) of a quadratic equation [106]. Thus, given an EOS for the SHS
fluid, one can get the thermodynamically consistent values of y(1) and χ and
determine from them all the coefficients appearing in Eq. (191).
Figure 12 shows the cavity function for η = 0.164 and τ = 0.13 as obtained
from Monte Carlo simulations [101] and as predicted by the PY and RFA
theories, the latter making use of the EOS recently proposed by Miller and
Frenkel [102]. It can be observed that the RFA is not only more accurate than
the PY approximation near r = 1 but also near r = 2. On the other hand,
none of these two approximations account for the singularities (delta-peaks
and/or discontinuities) of y(r) at r =
8/3, 5/3,
3, 2, . . . [100, 101].
1.0 1.5 2.0 2.5
=0.164, =0.13 MC
Fig. 12. Cavity function of a single component SHS fluid for η = 0.164 and τ = 0.13.
The solid line represents simulation data [101]. The dotted and dashed lines represent
the PY and RFA approaches, respectively.
Alternative Approaches to Hard-Sphere Liquids 49
4.2 Single Component Square-Well Fluids
Now we consider again the SW interaction potential (169) but for a single
fluid, i.e., σij = σ, ǫij = ǫ, Rij = R. Since no exact solution of the PY theory
for the SW potential is known, the application of the RFA method is more
challenging in this case than for HS and SHS fluids.
As in the cases of HS and SHS, the key quantity is the Laplace transform
of rg(r) defined by Eq. (98). It is again convenient to introduce the auxiliary
function Ψ(s) through Eq. (99). As before, the conditions g(r) = finite and
χ = finite imply Eqs. (102) and (104), respectively. However, the important
difference between HS and SHS fluids is that in the latter case G(s) must
reflect the fact that g(r) is discontinuous at r = R as a consequence of the
discontinuity of the potential φ(r) and the continuity of the cavity function
y(r). This implies that G(s), and hence Ψ(s), must contain the exponential
term e−(R−σ)s. This manifests itself in the low-density limit, where the con-
dition limρ→0 y(r) = 1 yields
Ψ(s) =
(1 + s)− e−(R−1)s(e1/T∗ − 1)(1 +Rs)
, (197)
where T ∗ ≡ kBT/ǫ and we have taken σ = 1.
In the spirit of the RFA method, the simplest form that complies with Eq.
(102) and is consistent with Eq. (197) is [107]
Ψ(s) =
−12η + E1s+ E2s2 + E3s3
1 +Q0 +Q1s− e−(R−1)s (Q0 +Q2s)
, (198)
where the coefficients Q0, Q1, Q2, E1, E2, and E3 are functions of η, T
∗, and
R. The condition (104) allows one to express the parameters Q1, E1, E2, and
E3 as linear functions of Q0 and Q2 [107, 108]:
1 + 2η
+ 2η(R3 − 1)Q2 −
(R − 1)2(R2 + 2R+ 3)Q0
+Q2 − (R − 1)Q0, (199)
1 + 2η
3− 4(R3 − 1)Q2 + (R− 1)2(R2 + 2R+ 3)Q0
, (200)
1 + 2η
{1− η − 2(R− 1) [1− 2ηR(R+ 1)]Q2
+(R− 1)2
(1− η(R + 1)2
, (201)
1 + 2η
(1 − η)2 + 6η(R− 1)
R+ 1− 2ηR2
−η(R − 1)2[4 + 2R− η(3R2 + 2R+ 1)]Q0
. (202)
50 M. López de Haro, S. B. Yuste and A. Santos
From Eq. (102), we have
g(1+) =
. (203)
The complete RDF is given by Eq. (122), where now Eq. (198) must be used
in Eq. (123). In particular, ψ1(r) and ψ2(r) are
ψ1(r) = ψ10(r)Θ(r) + ψ11(r + 1−R)Θ(r + 1−R), (204)
ψ2(r) = ψ20(r)Θ(r)+ψ21(r+1−R)Θ(r+1−R)+ψ22(r+2−2R)Θ(r+2−2R),
(205)
where
ψ1k(r) = 2π
W1k(si)
E′(si)
six, (206)
ψ2k(r) = −4π2
rW2k(si) +W
2k(si)−W2k(si)
E′′(si)
E′(si)
[E′(si)]2
. (207)
Here, si are the three distinct roots of E(s) ≡ −12η+E1s+E2s2 +E3s3 and
W10(s) ≡ 1 +Q0 +Q1s, W11(s) ≡ −(Q0 +Q2s). (208)
W20(s) ≡ s[W10(s)]2, W21(s) ≡ 2sW10(s)W11(s), W22(s) ≡ s[W11(s)]2.
(209)
To close the proposal, we need to determine the parameters Q0 and Q2 by
imposing two new conditions. An obvious condition is the continuity of the
cavity function at r = R, what implies
g(R+) = e1/T
g(R−). (210)
This yields (
1− e−1/T
ψ10(R− 1) = −ψ11(0) = 2π
. (211)
As an extra condition, we could enforce the continuity of the first derivative
y′(r) at r = R [109]. However, this complicates the problem too much without
any relevant gain in accuracy. In principle, it might be possible to impose
consistency with a given EOS, via either the virial route, the compressibility
route, or the energy route. But this is not practical since no simple EOS for
SW fluids is at our disposal for wide values of density, temperature, and range.
As a compromise between simplicity and accuracy, we fix the parameter Q0
at its exact zero-density limit value, namely Q0 = e
1/T∗ − 1 [107]. Therefore,
Eq. (211) becomes a transcendental equation for Q2 that needs to be solved
numerically. For narrow SW potentials, however, it is possible to replace the
exact condition (210) by a simpler one allowing Q2 to be obtained analytically
[108], which is especially useful for determining the thermodynamic properties
[108, 110].
Alternative Approaches to Hard-Sphere Liquids 51
R=1.5, 3=0.4, T*=1.5
1.0 1.5 2.0 2.5 3.0
2.0 R=2, 3=0.4, T*=3
R=1.05, 3=0.8, T*=0.5
Fig. 13. Radial distribution function of a single component SW fluid for R = 1.05,
ρσ3 = 0.8, and T ∗ = 0.5 (top panel), for R = 1.5, ρσ3 = 0.4, and T ∗ = 1.5 (middle
panel), and for R = 2.0, ρσ3 = 0.4, and T ∗ = 3.0 (bottom panel). The circles
represent simulation data [111] and the solid lines refer to the results obtained from
the RFA method.
It can be proven that the RFA proposal (198) reduces to the exact solutions
of the PY equation [29, 94] in the HS limit, i.e., ǫ → 0 or R → 1, and in the
SHS limit, i.e., ǫ→ ∞ and R → 1 with (R− 1)e1/T∗ = finite [107, 108].
Comparison with computer simulations [107, 108, 110, 111] shows that the
RFA for SW fluids is rather accurate at any fluid density if the potential well
is sufficiently narrow (say R ≤ 1.2), as well as for any width if the density
is small enough (say ρσ3 ≤ 0.4). However, as the width and/or the density
increase, the RFA predictions worsen, especially at low temperatures. As an
52 M. López de Haro, S. B. Yuste and A. Santos
illustration, Fig. 13 compares the RDF provided by the RFA with Monte Carlo
data [111] for three representative cases.
4.3 Hard Disks
As is well known, the PY theory is exactly solvable for HS fluids with an
odd number of dimensions [112–114]. In particular, in the case of hard rods
(d = 1), the PY theory provides the exact RDF g(r) or, equivalently, the exact
cavity function y(r) outside the hard core (i.e., for r > σ). However, it does
not reproduce the exact y(r) in the overlapping region (i.e., for r < σ) [85].
The full exact one-dimensional cavity function is [85]
yHR(r|η) =
e−(r−1)η/(1−η)
ηn−1e−(r−n)η/(1−η)
(1 − η)n(n− 1)!
(r − n)n−1Θ(r − n),
(212)
where the subscript HR stands for hard rods and, as usual, σ = 1 has been
taken. Consequently, one has
gHR(1
+|η) = 1
dr rhHR(r|η) ≡ H(0)HR(η) = −
η2. (213)
When d is even, the PY equation is not analytically solvable for the HS
interaction. In particular, in the important case of hard disks (d = 2), one
must resort to numerical solutions of the PY equation [1, 115]. Alternatively,
a simple heuristic approach has proven to yield reasonably good results [116].
Such an approach is based on the näıve assumption that the structure and
spatial correlations of a hard-disk fluid share some features with those of a
hard-rod and a hard-sphere fluid. This fuzzy idea becomes a more specific one
by means of the following simple model [116]:
gHD(r|η) = ν(η)gHR(r|ω1(η)η) + [1− ν(η)]gHS(r|ω3(η)η). (214)
Here, the subscript HD stands for hard disks (d = 2) and the subscript HS
stands for hard spheres (d = 3). The parameter ν(η) is a density-dependent
mixing parameter, while ω1(η)η and ω3(η)η are the packing fractions in one
and three dimensions, respectively, which are “equivalent” to the packing
fraction η in two dimensions. In Eq. (214), it is natural to take for gHR(r|η)
the exact solution, Eq. (212). As for gHR(r|η), one might use the RFA recipe
described in Section 3. However, in order to keep the model (214) as simple
as possible, it is sufficient for practical purposes to take the PY solution, Eq.
(108). In the latter approximation,
gHS(1
+|η) = 1 + η/2
(1− η)2
dr rhHS(r|η) ≡ H(0)HS (η) = −
10− 2η + η2
20(1 + 2η)
(215)
Alternative Approaches to Hard-Sphere Liquids 53
In order to close the model (214), we still need to determine the parameters
ν(η), ω1(η), and ω3(η). To that end, we first impose the condition that Eq.
(214) must be consistent with a prescribed contact value gHD(1
+|η) or, equiv-
alently, with a prescribed compressibility factor ZHD(η) = 1 + 2ηgHD(1
+|η),
with independence of the choice of the mixing parameter ν(η). In other words,
gHD(1
+|η) = gHR(1+|ω1(η)η) = gHS(1+|ω3(η)η). (216)
Making use of Eqs. (213) and (215), this yields
ω1(η) =
gHD(1
+|η)− 1
ηgHD(1+|η)
, ω3(η) =
4gHD(1
+|η) + 1−
24gHD(1+|η) + 1
4ηgHD(1+|η)
(217)
Once ω1(η) and ω3(η) are known, we can determine ν(η) by imposing that
the model (214) reproduces the isothermal compressibility χHD(η) thermody-
namically consistent with the prescribed ZHD(η) [cf. Eq. (97)]. From Eqs. (96)
and (214) one has
χHD(η) = 1 + 8η
dr r {ν(η)hHR(r|ω1(η)η) + [1− ν(η)] hHS(r|ω3(η)η)} ,
(218)
so that
ν(η) =
[χHD(η)− 1] /8η −H(0)HS (ω3(η)η)
HR(ω1(η)η) −H
HS (ω3(η)η)
, (219)
where H
HR(η) and H
HS (η) are given by Eqs. (213) and (215), respectively.
Once a sensible EOS for hard disks is chosen [see, for instance, Table 1],
Eqs. (217) and (219) provide the parameters of the model (214). The results
show that the scaling factor ω1(η) is a decreasing function, while ω3(η) is
an increasing function [116]. As for the mixing parameter ν(η), it is hardly
dependent of density and takes values around ν(η) ≃ 0.35–0.40.
Comparison of the interpolation model (214) with computer simulation re-
sults shows a surprisingly good agreement, despite the crudeness of the model
and the absence of empirical fitting parameters, especially at low and mod-
erate densities [116]. The discrepancies become important only for distances
beyond the location of the second peak and for densities close to the stability
threshold.
5 Perturbation Theory
When one wants to deal with realistic intermolecular interactions, the prob-
lem of deriving the thermodynamic and structural properties of the system
becomes rather formidable. Thus, perturbation theories of liquids have been
devised since the mid twentieth century. In the case of single component flu-
ids, the use of an accurate and well characterized RDF for the HS fluid in a
54 M. López de Haro, S. B. Yuste and A. Santos
perturbation theory opens up the possibility of deriving a closed theoretical
scheme for the determination of the thermodynamic and structural proper-
ties of more realistic models, such as the Lennard–Jones (LJ) fluid. In this
section, we will consider this model system, which captures the basic physical
properties of real non-polar fluids, to illustrate the procedure.
In the application of the perturbation theory of liquids, the stepping stone
has been the use of the HS RDF obtained from the solution to the PY equa-
tion. Unfortunately, the absence of thermodynamic consistency present in the
PY approximation (as well as in other integral equation theories) may clearly
contaminate the results derived from its use within a perturbative treatment.
In what follows we will reanalyze the different theoretical schemes for the ther-
modynamics of LJ fluids that have been constructed with perturbation theory,
taking as the reference system the HS fluid. This includes the consideration
of the RDF as obtained with the RFA method, which embodies thermody-
namic consistency, as well as the proposal of a unifying framework in which
all schemes fit in. With our development, we will be able to present a formula-
tion which lends itself to relatively easy numerical calculations while retaining
the merits that analytical results provide, namely a detailed knowledge and
control of all the approximations involved.
Let us consider a three-dimensional fluid system defined by a pair inter-
action potential φ(r). The virial and energy EOS express the compressibility
factor Z and the excess part of the Helmholtz free energy per unit volume
f ex, respectively, in terms of the RDF of the system as
Z = 1− 2
∂φ(r)
g(r)r3, (220)
= 2πρβ
dr φ(r)g(r)r2 , (221)
where β ≡ 1/kBT . Let us now assume that φ(r) is split into a known (ref-
erence) part φ0(r) and a perturbation part φ1(r). The usual perturbative
expansion for the Helmholtz free energy to first order in β leads to [2]
+ 2πρβ
dr φ1(r)g0(r)r
, (222)
where f0 and g0(r) are the free energy and the RDF of the reference system,
respectively.
The LJ potential is
φLJ(r) = 4ǫ
r−12 − r−6
, (223)
where ǫ is the depth of the well and, for simplicity, we have taken the distance
at which the potential vanishes as the length unit, i.e., φLJ(r = 1) = 0. For
this potential the reference system may be forced to be a HS system, i.e., one
can set
Alternative Approaches to Hard-Sphere Liquids 55
φ0(r) = φHS(r) =
∞, r ≤ σ0,
0, r > σ0,
(224)
where σ0 is a conveniently chosen effective HS diameter. In this case the
Helmholtz free energy to this order is approximated by
≈ fHS
+ 2πρβ
dr φLJ(r)gHS(r/σ0)r
2. (225)
Note that Eq. (225) may be rewritten in terms of the Laplace transform G(s)
of (r/σ0)gHS(r/σ0) as
≈ fHS
+ 2πρβσ30
ds ΦLJ(s)G(s), (226)
where ΦLJ(s) satisfies
rφLJ(r) = σ0
ds e−rs/σ0ΦLJ(s), (227)
so that
ΦLJ(s) = 4ǫσ
(s/σ0)
− (s/σ0)
. (228)
Irrespective of the value of the diameter σ0 of the reference system, the
right hand side of Eq. (226) represents always an upper bound for the value of
the free energy of the real system. Therefore, it is natural to determine σ0 so
as to provide the least upper bound. This is precisely the variational scheme of
Mansoori and Canfield [117,118] and Rasaiah and Stell [119], usually referred
to as MC/RS, and originally implemented with the PY theory for G(s), Eq.
(108). In our case, however, we will considerG(s) as given by the RFA method,
Eq. (113). Therefore, at fixed ρ and β, the effective diameter σ0 in the MC/RS
scheme is obtained from the conditions
{∫ η0
ZHS(η)− 1
+ 48βǫσ−20
dsG(s|η0)
(s/σ0)
− (s/σ0)
= 0, (229)
{∫ η0
ZHS(η)− 1
+ 48βǫσ−20
dsG(s|η0)
(s/σ0)
− (s/σ0)
> 0. (230)
In these equations, use has been made of the thermodynamic relationship
between the free energy and the compressibility factor, Eq. (78). Moreover, we
have called η0 ≡ (π/6)ρσ30 and have made explicit with the notation G(s|η0)
the fact that the HS RDF depends on the packing fraction η0.
56 M. López de Haro, S. B. Yuste and A. Santos
Even if the reference system is not forced to be a HS fluid, one can still use
Eq. (226) provided an adequate choice for σ0 is made such that the expansion
involved in the right hand side of Eq. (222) yields the right hand side of Eq.
(226) to order β2. This is the idea of the Barker and Henderson [120] first
order perturbation scheme (BH1), where the effective HS diameter is
1− e−βφLJ(r)
. (231)
The same ideas may be carried out to higher order in the perturbation
expansion. The inclusion of the second order term in the expansion yields the
so-called macroscopic compressibility approximation [2] for the free energy,
namely
+ 2πρβ
dr φ1(r)g0(r)r
−πρβ2χ0
dr φ21(r)g0(r)r
, (232)
where χ0 is the (reduced) isothermal compressibility of the reference system
[121].
To implement a particular perturbation scheme in this approximation un-
der a unifying framework that eventually leads to easy numerical evaluation,
two further assumptions may prove convenient. First, the perturbation poten-
tial φ1(r) ≡ φLJ(r)−φ0(r) may be split into two parts using some “molecular
size” parameter ξ ≥ σ0 such that
φ1(r) =
φ1a(r), 0 ≤ r ≤ ξ,
φ1b(r), r > ξ.
(233)
Next, a choice for the RDF for the reference system is done in the form
g0 (r) ≈ θ(r)yHS(r/σ0), (234)
where yHS is the cavity (background) correlation function of the HS system
and θ(r) is a step function defined by
θ(r) =
θa(r), 0 ≤ r ≤ ξ,
θb(r), r > ξ,
(235)
in which the functions θa(r) and θb(r) depend on the scheme.
With these assumptions the integrals involved in Eq. (232) may be rewrit-
ten as
dr φn1 (r)g0(r)r
dr φn1a(r)θa(r)yHS(r/σ0)r
dr φn1a(r)θa(r)gHS(r/σ0)r
dr φn1b(r)θb(r)gHS(r/σ0)r
2, (236)
Alternative Approaches to Hard-Sphere Liquids 57
with n = 1, 2 and where the fact that yHS(r/σ0) = gHS(r/σ0) when r > σ0
has been used. Decomposing the last integral as
and applying
the same step as in Eq. (226), Eq. (236) becomes
In = σ
ds Φnb(s)G(s) +
dr φn1a(r)θa(r)yHS(r/σ0)r
dr [φn1a(r)θa(r) − φn1b(r)θb(r)] gHS(r/σ0)r2, (237)
where the functions Φ1b(s) and Φ2b(s) are defined by the relation
rφn1b(r)θb(r) = σ0
ds e−rs/σ0Φnb(s). (238)
In the Barker–Henderson second order perturbation scheme (BH2), one
takes
θa(r) = 0, θb(r) = 1, ξ = σ0, φ1a(r) = 0, φ1b(r) = 4ǫ
r−12 − r−6
(239)
and σ0 is computed according to Eq. (231). This choice ensures that
+ 2πρβ
dr φ1(r)gHS(r/σ0)r
−πρβ2χHS
dr φ21(r)gHS(r/σ0)r
. (240)
On the other hand, if one chooses
θa(r) = exp [−β (φLJ(r) + ǫ)] , θb(r) = 1, ξ = 21/6, (241)
φ1a(r) = −ǫ, φ1b(r) = 4ǫ
r−12 − r−6
, (242)
the scheme leads to the Weeks–Chandler–Andersen (WCA) theory [122] if one
determines the HS diameter through the condition χ0 = χHS [123], which in
turn implies
dr r2e−βφ0(r)yHS(r/σ0) =
∫ 21/6
dr r2gHS(r/σ0)
1− e−βφ0(r)
. (243)
To close the scheme, the HS cavity function has to be provided in the range 0 ≤
r ≤ σ0. Fortunately, relatively simple expressions for yHS(r/σ0) are available
in the literature [124–126], apart from our own proposal, Eq. (140).
Note that θb(r) and φ1b(r), and thus also Φnb(s), are the same functions in
the BH2 and WCA schemes. It is convenient, in order to have all the quantities
needed to evaluate fLJ in these schemes, to provide explicit expressions for
Φ1b(s) and Φ2b(s). These are given by [cf. Eq. (228)]
58 M. López de Haro, S. B. Yuste and A. Santos
Φ1b(s) = ΦLJ(s), (244)
Φ2b(s) = 16ǫ
2σ−20
(s/σ0)
− 2(s/σ0)
(s/σ0)
. (245)
Up to this point, we have embodied the most popular perturbation schemes
within a unified framework that requires as input only the EOS of the HS fluid
in order to compute the Helmholtz free energy of the LJ system and leads
to relatively easy numerical computations. It should be clear that a variety
of other possible schemes, requiring the same little input, fit in our unified
framework, which is based on the RFA method for gHS(r/σ0) and G(s). Once
fLJ has been determined, the compressibility factor of the LJ fluid at a given
order of the perturbation expansion readily follows from Eqs. (222) or (232)
through the thermodynamic relation
ZLJ = ρ
. (246)
Taking into account that the HS fluid presents a fluid-solid transition at
a freezing packing fraction ηf ≃ 0.494 [127] and a solid-fluid transition at
a melting packing fraction ηm ≃ 0.54 [127], the fluid-solid and solid-fluid
coexistence lines for the LJ system may be computed from the values (ρ, T )
determined from the conditions (π/6)ρσ30(ρ, T ) = ηf and (π/6)ρσ
0(ρ, T ) = ηm,
respectively, with the effective diameter σ0(ρ, T ) obtained using any of the
perturbative schemes. Similarly, admitting that there is a glass transition in
the HS fluid at the packing fraction ηg ≃ 0.56 [128], one can now determine
the location of the liquid-glass transition line for the LJ fluid in the (ρ, T )
plane from the simple relationship (π/6)ρσ30(ρ, T ) = ηg. With a proper choice
for ZHS, it has been shown [76,129,130] that the critical point, the structure,
and the phase diagram (including a glass transition) of the LJ fluid may be
adequately described with this approach.
6 Perspectives
In this chapter we have given a self-contained account of a simple (mostly
analytical) framework for the study of the thermodynamic and structural
properties of hard-core systems. Whenever possible, the developments have
attempted to cater for mixtures with an arbitrary number of components (in-
cluding polydisperse systems) and arbitrary dimensionality. We started con-
sidering the contact values of the RDF because they enter directly into the
EOS and are required as input in the RFA method to compute the structural
properties. With the aid of consistency conditions, we were able to devise var-
ious approximate proposals which, when used in conjunction with a sensible
choice for the contact value of the RDF of the single component fluid (re-
quired in the formulation but otherwise chosen at will), have been shown to
Alternative Approaches to Hard-Sphere Liquids 59
be in reasonably good agreement with simulation results and lead to accurate
EOS both for additive and non-additive mixtures. Some aspects of the results
that follow from the use of these EOS were illustrated by looking at demixing
problems in these mixtures, including the far from intuitive case of a binary
mixture of non-additive hard spheres in infinite dimensionality.
After that, restricting ourselves to three-dimensional systems, we described
the RFA method as applied to a single component hard-sphere fluid and to
a multicomponent mixture of HS. Using this approach, we have been able to
obtain explicit analytical results for the RDF, the direct correlation function,
the static structure factor, and the bridge function, in the end requiring as
input only the contact value of the RDF of the single component HS fluid
(or equivalently its compressibility factor). One of the nice assets of the RFA
approach is that it eliminates the thermodynamic consistency problem which
is present in most of the integral equation formulations for the computation
of structural quantities. Once again, when a sensible choice for the single
component EOS is made, we have shown, through the comparison between
the results of the RFA approach and simulation data for some illustrative
cases, the very good performance of our development. Also, the use of the
RFA approach in connection with some other related systems (sticky hard
spheres, square-well fluids, and hard disks) has been addressed.
The final part of the chapter concerns the use of HS results for more
realistic intermolecular potentials in the perturbation theory of liquids. In
this instance we have been able to provide a unifying scheme in which the
most popular perturbation theory formulations may be expressed and which
was devised to allow for easy computations. We illustrated this for a LJ fluid
but it should be clear that a similar approach might be followed for other fluids
and in fact it has recently been done in connection with the glass transition
of hard-core Yukawa fluids [131].
Finally, it should be clear that there are many facets of the equilibrium and
structural properties of hard-core systems that may be studied with a simi-
lar approach but that up to now have not been considered. For instance, the
generalizations of the RFA approach for systems such as hard hyperspheres,
non-additive hard spheres, square-well mixtures, penetrable spheres [132], or
the Jagla potential [133] appear as interesting challenges. Similarly, the ex-
tension of the perturbation theory scheme to the case of LJ mixtures seems
a worthwhile task. We hope to address some of these problems in the future
and would be very much rewarded if some others were taken up by researchers
who might find these developments also a valuable tool for their work.
References
1. J. A. Barker and D. Henderson, Rev. Mod. Phys. 48, 587 (1976).
2. D. A. McQuarrie, Statistical Mechanics (Harper & Row, N. Y., 1976).
3. H. L. Friedman, A Course in Statistical Mechanics (Prentice Hall, Englewood
Cliffs, 1985).
60 M. López de Haro, S. B. Yuste and A. Santos
4. J.-P. Hansen and I. R. McDonald, Theory of Simple Liquids, (Academic Press,
London, 1986).
5. J. L. Lebowitz and D. Zomick, J. Chem. Phys. 54, 3335 (1971).
6. J. T. Jenkins and F. Mancini, J. Appl. Mech. 54, 27 (1987).
7. C. Barrio and J. R. Solana, J. Chem. Phys. 115, 7123 (2001); 117, 2451(E)
(2002).
8. J. L. Lebowitz, Phys. Rev. A 133, 895 (1964).
9. H. Reiss, H. L. Frisch, and J. L. Lebowitz, J. Chem. Phys. 31, 369 (1959); E.
Helfand, H. L. Frisch, and J. L. Lebowitz, J. Chem. Phys. 34, 1037 (1961); J.
L. Lebowitz, E. Helfand, and E. Praestgaard, J. Chem. Phys. 43, 774 (1965).
10. M. J. Mandell and H. Reiss, J. Stat. Phys. 13, 113 (1975).
11. Y. Rosenfeld, J. Chem. Phys. 89, 4272 (1988).
12. M. Heying and D. S. Corti, J. Phys. Chem. B 108, 19756 (2004).
13. T. Boubĺık, J. Chem. Phys. 53, 471 (1970).
14. E. W. Grundke and D. Henderson, Mol. Phys. 24, 269 (1972).
15. L. L. Lee and D. Levesque, Mol. Phys. 26, 1351 (1973).
16. G. A. Mansoori, N. F. Carnahan, K. E. Starling, and J. T. W. Leland, J. Chem.
Phys. 54, 1523 (1971).
17. D. Henderson, A. Malijevský, S. Lab́ık, and K. Y. Chan, Mol. Phys. 87, 273
(1996).
18. D. H. L. Yau, K.-Y. Chan, and D. Henderson, Mol. Phys. 88, 1237 (1996); 91,
1137 (1997).
19. D. Henderson and K. Y. Chan, J. Chem. Phys. 108, 9946 (1998); Mol. Phys.
94, 253 (1998); 98, 1005 (2000).
20. D. Henderson, D. Boda, K. Y. Chan, and D. T. Wasan, Mol. Phys. 95, 131
(1998).
21. D. Matyushov, D. Henderson, and K.-Y. Chan, Mol. Phys. 96, 1813 (1999).
22. D. Cao, K.-Y. Chan, D. Henderson, and W. Wang, Mol. Phys. 98, 619 (2000).
23. D. V. Matyushov and B. M. Ladanyi, J. Chem. Phys. 107, 5815 (1997).
24. C. Barrio and J. R. Solana, J. Chem. Phys. 113, 10180 (2000).
25. D. Viduna and W. R. Smith, Mol. Phys. 100, 2903 (2002); J. Chem. Phys.
117, 1214 (2002).
26. D. Henderson, Mol. Phys. 30, 971 (1975).
27. A. Santos, M. López de Haro, and S. B. Yuste, J. Chem. Phys. 103, 4622
(1995); M. López de Haro, A. Santos, and S. B. Yuste, Eur. J. Phys. 19, 281
(1998).
28. S. Luding, Phys. Rev. E 63, 042201 (2001); S. Luding, Adv. Compl. Syst.
4, 379 (2002); S. Luding and O. Strauß, in Granular Gases, T. Pöschel and
S. Luding, eds. (LNP 564, Springer-Verlag, Berlin, 2001), pp. 389–409.
29. M. S. Wertheim, Phys. Rev. Lett. 10, 321 (1963); E. Thiele, J. Chem. Phys.
39, 474 (1963).
30. N. F. Carnahan and K. E. Starling, J. Chem. Phys. 51, 635 (1969).
31. M. Luban and J. P. J. Michels, Phys. Rev. A 41, 6796 (1990).
32. E. Hamad, J. Chem. Phys. 101, 10195 (1994).
33. C. Vega, J. Chem. Phys. 108, 3074 (1998).
34. N. M. Tukur, E. Z. Hamad, and G. A. Mansoori, J. Chem. Phys. 110, 3463
(1999).
35. A. Santos, S. B. Yuste, and M. López de Haro, Mol. Phys. 96, 1 (1999).
36. A. Malijevský and J. Veverka, Phys. Chem. Chem. Phys. 1, 4267 (1999).
Alternative Approaches to Hard-Sphere Liquids 61
37. A. Santos, S. B. Yuste, and M. López de Haro, Mol. Phys. 99, 1959 (2001).
38. M. González-Melchor, J. Alejandre, and M. López de Haro, J. Chem. Phys.
114, 4905 (2001).
39. M. López de Haro, S. B. Yuste, and A. Santos, Phys. Rev. E 66, 031202
(2002).
40. A. Santos, Mol. Phys. 96, 1185 (1999); 99, 617(E) (2001).
41. C. Regnaut, A. Dyan, and S. Amokrane, Mol. Phys. 99, 2055 (2001); 100,
2907(E) (2002).
42. A. Santos, S. B. Yuste, and M. López de Haro, J. Chem. Phys. 117, 5785
(2002).
43. A. Santos, S. B. Yuste and M. López de Haro, J. Chem. Phys. 123, 234512
(2005); M. López de Haro, S. B. Yuste, and A. Santos, Mol. Phys. 104, 3461
(2006).
44. S. Luding and A. Santos, J. Chem. Phys. 121, 8458 (2004).
45. M. Barošová, A. Malijevský, S. Lab́ık, and W. R. Smith, Mol. Phys. 87, 423
(1996).
46. H. Hansen-Goos and R. Roth, J. Chem. Phys. 124, 154506 (2006).
47. R. Evans, in Liquids and Interfaces, edited by J. Charvolin, J. F. Joanny, and
J. Zinn-Justin (North-Holland, Amsterdam, 1990).
48. Y. Rosenfeld, Phys. Rev. Lett. 63, 980 (1989).
49. A. Malijevský, M. Barošová, and W. R. Smith, Mol. Phys. 91, 65 (1997).
50. Al. Malijevský, A. Malijevský, S. B. Yuste, A. Santos, and M. López de Haro,
Phys. Rev. E 66, 061203 (2002).
51. M. Buzzacchi, I. Pagonabarraga, and N. B. Wilding, J. Chem. Phys. 121, 11362
(2004).
52. Al. Malijevský, S. B. Yuste, A. Santos, and M. López de Haro, preprint arXiv:
cond-mat/0702284.
53. F. Lado, Phys. Rev. E 54, 4411 (1996).
54. I. Prigogine and S. Lafleur, Bull. Classe Sci. Acad. Roy. Belg. 40, 484, 497
(1954).
55. S. Asakura and F. Oosawa, J. Chem. Phys. 22, 1255 (1954); J. Polym. Sci. 33,
183 (1958).
56. R. Kikuchi, J. Chem. Phys. 23, 2327 (1955).
57. P. Ballone, G. Pastore, G. Galli, and D. Gazzillo, Mol. Phys. 59, 275 (1986).
58. D. Gazzillo, G. Pastore, and S. Enzo, J. Phys.: Condens. Matter 1, 3469 (1989);
D. Gazzillo, G. Pastore, and R. Frattini, J. Phys.: Condens. Matter 2,8465
(1990).
59. A. Santos, M. López de Haro, and S. B. Yuste, J. Chem. Phys. 122, 024514
(2005).
60. E. Z. Hamad, J. Chem. Phys. 105, 3229 (1996).
61. E. Z. Hamad, J. Chem. Phys. 105, 3222 (1996).
62. H. Hammawa and E. Z. Hamad, J. Chem. Soc. Faraday Trans. 92, 4943 (1996).
63. M. Al-Naafa, J. B. El-Yakubu, and E. Z. Hamad, Fluid Phase Equilibria 154,
33 (1999).
64. J. Jung, M. S. Jhon, and F. H. Ree, J. Chem. Phys. 100, 528 (1994).
65. J. Jung, M. S. Jhon, and F. H. Ree, J. Chem. Phys. 100, 9064 (1994).
66. T. Coussaert and M. Baus, J. Chem. Phys. 109, 6012 (1998).
67. A. Yu. Vlasov and A. J. Masters, Fluid Phase Equilibria 212, 183 (2003).
68. M. López de Haro and C. F. Tejero, J. Chem. Phys. 121, 6918 (2004).
62 M. López de Haro, S. B. Yuste and A. Santos
69. S. B. Yuste, A. Santos, and M. López de Haro, Europhys. Lett. 52, 158 (2000).
70. H.-O. Carmesin, H. L. Frisch, and J. K. Percus, J. Stat. Phys. 63, 791 (1991).
71. A. Santos and M. López de Haro, Phys. Rev. E 72, 010501(R) (2005).
72. R. Roth, R. Evans, and A. A. Louis, Phys. Rev. E 64, 051202 (2001).
73. S. B. Yuste and A. Santos, Phys. Rev. A 43, 5418 (1991).
74. S. B. Yuste, M. López de Haro, and A. Santos, Phys. Rev. E 53, 4820 (1996).
75. M. Robles, M. López de Haro, A. Santos, and S. B. Yuste, J. Chem. Phys. 108,
1290 (1998).
76. M. Robles and M. López de Haro, Europhys. Lett. 62, 56 (2003).
77. M. Robles and M. López de Haro, J. Chem. Phys. 107, 4648 (1997).
78. E. Waisman, Mol. Phys. 25, 45 (1973); D. Henderson and L. Blum, Mol. Phys.
32, 1627 (1976); J. S. Høye and L. Blum, J. Stat. Phys. 16, 399 (1977).
79. A. Dı́ez, J. Largo, and J. R. Solana, J. Chem. Phys. 125, 074509 (2006).
80. J. Kolafa, S. Lab́ık, and A. Malijevský, Phys. Chem. Chem. Phys. 6, 2335
(2004). See also http://www.vscht.cz/fch/software/hsmd/ for molecular dy-
namics results of g(r).
81. A. Trokhymchuk, I. Nezbeda, J. Jirsák, and D. Henderson, J. Chem. Phys.
123, 024501 (2005).
82. M. López de Haro, A. Santos, and S. B. Yuste, J. Chem. Phys. 124, 236102
(2006).
83. L. L. Lee, J. Chem. Phys. 103, 9388 (1995); L. L. Lee, D. Ghonasgi, and E.
Lomba, J. Chem. Phys. 104, 8058 (1996); L. L. Lee and A. Malijevský, J.
Chem. Phys. 114, 7109 (2001).
84. S. Lab́ık and A. Malijevský, Mol. Phys. 53, 381 (1984).
85. Al. Malijevský and A. Santos, J. Chem. Phys. 124, 074508 (2006).
86. A. Santos and Al. Malijevský, Phys. Rev. E 75, 021201 (2007).
87. S. B. Yuste, A. Santos, and M. López de Haro, J. Chem. Phys. 108, 3683
(1998).
88. L. Blum and J. S. Høye, J. Phys. Chem. 81, 1311 (1977).
89. J. Abate and W. Whitt, Queuing Systems 10, 5 (1992).
90. A code using the Mathematica computer algebra system to obtain Gij(s)
and gij(r) with the present method is available from the web page
http://www.unex.es/eweb/fisteor/santos/filesRFA.html.
91. N. W. Ashcroft and D. C. Langreth, Phys. Rev. 156, 685 (1967).
92. A. B. Bathia and D. E. Thornton, Phys. Rev. B 8, 3004 (1970).
93. S. B. Yuste, A. Santos, and M. López de Haro, Mol. Phys. 98, 439 (2000).
94. R. J. Baxter, J. Chem. Phys. 49, 2270 (1968).
95. J. W. Perram and E. R. Smith, Chem. Phys. Lett. 35, 138 (1975).
96. B. Barboy, Chem. Phys. 11, 357 (1975); B. Barboy and R. Tenne, Chem. Phys.
38, 369 (1979).
97. G. Stell, J. Stat. Phys. 63, 1203 (1991); B. Borštnik, C. G. Jesudason, and G.
Stell, J. Chem. Phys. 106, 9762 (1997).
98. B. Barboy, J. Chem. Phys. 61, 3194 (1974).
99. J. W. Perram and E. R. Smith, Chem. Phys. Lett. 39, 328 (1975); P. T.
Cummings, J. W. Perram, and E. R. Smith, Mol. Phys. 31, 535 (1976); E. R.
Smith and J. W. Perram, J. Stat. Phys. 17, 47 (1977); J. W. Perram and E.
R. Smith, Proc. R. Soc. London A353, 193 (1977); W. G. T. Kranendonk and
D. Frenkel, Mol. Phys. 64, 403 (1988); C. Regnaut and J. C. Ravey, J. Chem.
Phys. 91, 1211 (1989); G. Stell and Y. Zhou, J. Chem. Phys. 91, 3618 (1989);
Alternative Approaches to Hard-Sphere Liquids 63
J. N. Herrera and L. Blum, J. Chem. Phys. 94, 6190 (1991); A. Jamnik, D.
Bratko, and D. J. Henderson, J. Chem. Phys. 94, 8210 (1991); S. V. G. Menon,
C. Manohar, and K. S. Rao, J. Chem. Phys.; 95, 9186 (1991); Y. Zhou and G.
Stell, J. Chem. Phys. 96, 1504 (1992); E. Dickinson, J. Chem. Soc. Faraday
Trans. 88, 3561 (1992); C. F. Tejero and M. Baus, Phys. Rev. E, 48, 3793
(1993); K. Shukla and R. Rajagopalan, Mol. Phys. 81, 1093 (1994); C. Regnaut,
S. Amokrane, and Y. Heno, J. Chem. Phys. 102, 6230 (1995); C. Regnaut, S.
Amokrane, and P. Bobola, Prog. Colloid Polym. Sci. 98, 151 (1995); Y. Zhou,
C. K. Hall, and G. Stell, Mol. Phys. 86, 1485 (1995); J. N. Herrera-Pacheco
and J. F. Rojas-Rodŕıguez, Mol. Phys. 86, 837 (1995); Y. Hu, H. Liu, and J. M.
Prausnitz, J. Chem. Phys. 104, 396 (1996); O. Bernard and L. Blum, J. Chem.
Phys. 104, 4746 (1996); L. Blum, M. F. Holovko, and I. A. Protsykevych, J.
Stat. Phys. 84, 191 (1996); S. Amokrane, P. Bobola and C. Regnaut, Prog.
Colloid Polym. Sci. 100, 186 (1996); S. Amokrane and C. Regnaut, J. Chem.
Phys. 106, 376 (1997); C. Tutschka, G. Kahl, and E. Riegler, Mol. Phys. 100,
1025 (2002); D. Gazzillo and A. Giacometti, Mol. Phys. 100, 3307 (2002); M.
A. Miller and D. Frenkel, Phys. Rev. Lett. 90, 135702 (2003); D. Gazzillo and
A. Giacometti, J. Chem. Phys. 120, 4742 (2004); R. Fantoni, D. Gazzillo, and
A. Giacometti, Phys. Rev. E 72, 011503 (2005); A. Jamnik, Chem. Phys. Lett.
423, 23 (2006).
100. A. J. Post and E. D. Glandt, J. Chem. Phys. 84, 4585 (1986); N. A. Seaton
and E. D. Glandt, J. Chem. Phys. 84, 4595 (1986); 86, 4668 (1986); 87, 1785
(1987).
101. M. A. Miller and D. Frenkel, J. Phys.: Condens. Matter 16, S4901 (2004).
102. M. A. Miller and D. Frenkel, J. Chem. Phys. 121, 535 (2004).
103. Al. Malijevský, S. B. Yuste, and A. Santos, J. Chem. Phys. 125, 074507 (2006).
104. A. Santos, S. B. Yuste, and M. López de Haro, J. Chem. Phys. 109, 6814
(1998).
105. S. B. Yuste and A. Santos, J. Stat. Phys. 72, 703 (1993).
106. S. B. Yuste and A. Santos, Phys. Rev. E 48, 4599 (1993).
107. S. B. Yuste and A. Santos, J. Chem. Phys. 101, 2355 (1994).
108. L. Acedo and A. Santos, J. Chem. Phys. 115, 2805 (2001).
109. L. Acedo, J. Stat. Phys. 99, 707 (2000).
110. J. Largo, J. R. Solana, L. Acedo, and A. Santos, Mol. Phys. 101, 2981 (2003).
111. J. Largo, J. R. Solana, S. B. Yuste, and A. Santos, J. Chem. Phys. 122, 084510
(2005).
112. C. Freasier and D. J. Isbister, Mol. Phys. 42, 927 (1981).
113. E. Leutheusser, Physica A 127, 667 (1984).
114. M. Robles, M. López de Haro, and A. Santos, J. Chem. Phys. 120, 9113 (2004).
115. D. G. Chae, F. H. Ree, and T. Ree, J. Chem. Phys. 50, 1581 (1976).
116. S. B. Yuste and A. Santos, J. Chem. Phys. 99, 2020 (1993).
117. G. A. Mansoori and F. B. Canfield, J. Chem. Phys. 51, 4958 (1969).
118. G. A. Mansoori, J. A. Provine, and F. B. Canfield, J. Chem. Phys. 51, 5295
(1969).
119. J. Rasaiah and G. Stell, Mol. Phys. 18, 249 (1970).
120. J. A. Barker and D. Henderson, J. Chem. Phys. 47, 2856 (1967).
121. The macroscopic compressibility approach is only one of the possibilities of ap-
proximation to the second order Barker–Henderson perturbation theory term.
Another successful approach is the local-compressibility approximation (see
64 M. López de Haro, S. B. Yuste and A. Santos
Ref. [2], p. 308). This expresses the free energy in terms of φ1(r) and HS
quantities.
122. J. D. Weeks, D. Chandler, and H. C. Andersen, J. Chem. Phys. 53, 149 (1971).
123. A simple algorithm to compute a rather accurate approximation for the HS
diameter σ0 in the WCA theory has been given in L. Verlet and J. J. Weis,
Phys. Rev. A 5, 939 (1972).
124. D. Henderson and E. W. Grundke, J. Chem. Phys. 63, 601 (1975).
125. J. A. Ballance and R. J. Speedy, Mol. Phys. 54, 1035 (1985).
126. Y. Zhou and G. Stell, J. Stat. Phys. 52, 1389 (1988).
127. J.-P. Hansen and L. Verlet, Phys. Rev. 184, 151 (1969).
128. R. J. Speedy, J. Chem. Phys. 100, 6684 (1994).
129. M. Robles and M. López de Haro, Phys. Chem. Chem. Phys. 3, 5528 (2001).
130. M. López de Haro and M. Robles, J. Phys.: Condens. Matt. 16, S2089 (2004).
131. M. López de Haro and M. Robles, Physica A 372, 307 (2006).
132. C. N. Likos, Phys. Rep. 348, 267 (2001).
133. E. A. Jagla, J. Chem. Phys. 111, 8980 (1999).
