0
dreiq·rul(r) , (116)
of the perturbational part ofu(r). Eq. (115) can be rewritten in the mnemonic form 1
S(q) ≈ 1
S0(q) +βρul(q) . (117)
The RPA is the most simple perturbation theory for fluid microstructures, usually suited only for small wave numbers. It has been successfully used, also in its multi-component extension, for the calculation of (partial) structure factors in dense systems with (ultra-) soft potentials where u0(r) = 0, such as polymer blends and star polymers [25, 26]. For an ideal gas as reference system, one hasc0(r) = 0and S0(q) = 1. The RPA reduces then to the MSA for point-like particles, referred to in the literature as a version of the Debye-H¨uckel approximation. One should notice that the perturbation must be sufficiently weak, or the density sufficiently low, to ensure thatβρul(q)S0(q) > −1 with S(q) > 0. Furthermore, the RPA does not ensure that g(r < σ) = 0in case of a hard-sphere reference system. This non-physical feature of the RPA is related to an ambiguity in the choice of the perturbation potentialul(r)in Eq. (116) forr < σ.
The trueg(r)should not depend on this choice.
In the optimized random phase approximation (ORPA), ul(r) is extended into the hard-core regimer < σsuch thatg(r < σ) = 0. While the ORPA is a considerable improvement of the RPA, there is a price to pay in form of a much larger numerical effort to calculateS(q).
3.4 Percus-Yevick solution for hard spheres
Hard spheres serve as a reference system in the theory of uncharged liquids, as an ideal gas does in the theory of dilute gases, and a harmonic solid in solid-state physics. The PY approxima-tion leads to an integral equaapproxima-tion for the hard-sphere cavity funcapproxima-tiony(r)which can be solved analytically.
The solution proceeds as follows. For hard spheres, c(r) =g(r)[
1−eβu(r)]
= 0, r > σ (118)
in PY approximation, i.e. the hard-sphere direct correlation is approximately set equal to zero for non-overlap distances. As a matter of fact, the truec(r)has a small but non-vanishing tail forr > σ. As can be noticed here, the PY closure is identical to the MSA closure (cf. Eq. (92)) in the case of a hard-sphere fluid.
The hard-sphere cavity function reads where the lower equality follows from the PY closure in Eq. (97). It follows thatc(r)andg(r) share a jump discontinuity of equal magnitude atr=σ, withg(r=σ+) =−c(r=σ−), since y(r)is continuous everywhere.
Upon inserting the hard-sphere potential into Eq. (100), one obtains a quadratic integral equa-tion fory(r)of the form
y(r) = 1 +ρ is known already. Following Wertheim [27], we use a third-order polynomial
c(r) = a0+a1r+a2r2+a3r3 (121) as a trial solution ofc(r < σ), with yet unknown density-dependent expansion coefficients{ai}. This ansatz is suggested from the low density form ofc(r), which is a third order polynomial in case of hard spheres. For a proof of this statement use Eq. (65) to show that to first order in density,y(r)is given by
y(r) = 1 +ρ
∫
dr′f(r′)f(|r−r′|) +O(ρ2) . (122) For hard spheres, the Mayer-f function is f(r) = −1 for r < σ and zero otherwise. The convolution integral in Eq. (122) is then equal to the volume of overlap of two spheres of equal radiiσwith centres separated byr. As a consequence
y(r) = 1 + 8 Φ with x = r/σ. The overlap volume is zero for x > 2 as expressed by the unit step function θ(2−x). Recall that the PY approximation describesg(r)exactly to first order inρ. In using the polynomial ansatz in Eq. (121), it is assumed that the functional form ofc(r)is valid for all volume fractions.
The four expansion coefficients,{ai}, are determined by employing the continuity ofy(r)and its first two derivatives at r = σ. Their continuity follows from Eq. (120) and its first two derivatives. A fourth condition follows from the PY integral equation (120) evaluated atr = 0:
y(0) = 1 +ρ
∫
r′<σ
dr′y(r′) . (124)
After inserting Eq. (121) in Eq. (120) and making use of the four boundary conditions to determine the{ai}, a lengthy calculation gives the following PY result for the hard-spherec(r):
c(r < σ) = −
with
λ1 = (1 + 2Φ)2
(1−Φ)4 , λ2 =−6Φ (1 + 0.5Φ)2
(1−Φ)4 . (126)
The PY result for −c(r) reduces, for small Φ, to the correct first order density form of y(r) given in Eq. (123).
r/ σ c(r)
r/ σ
0 1 2 3 4 5
g(r)
Fig. 13: Percus-Yevick direct correlation function (left) and radial distribution function (right) of hard spheres.
Fourier transformation of c(r) leads with Eq. (72) to an analytic expression for S(q). This expression reads explicitly [3]
S(y) = 1
X2(y) +Y2(y) (127)
with
X(y) = 1−12 Φ [A f1(y) +B f2(y) ] (128) Y(y) = −12 Φ [A f3(y) +B f4(y) ] , (129) where
A= 1 + 2Φ
(1−Φ)2 , B = 1 + 0.5Φ
(1−Φ)2 , (130)
and
f1(y) = y−sin(y)
y3 , f2(y) = cos(y)−1
y2 (131)
f3(y) = f2(y) y + 1
2y, f4(y) = −yf1(y) . (132) We have introduced here the reduced wave numbery=qσ.
The reduced isothermal compressibility follows as
qlim→0S(q) = (1−Φ)4
(1 + 2Φ)2 , (133)
which is a monotonically decreasing function inΦ. For given analyticalS(q), the hard-sphere g(r)can be determined in principle by numerical Fourier-inversion. However, to avoid prob-lems caused by the jump discontinuity ing(r), it is safer to calculate first the function γ(r) :=
r/ σ S(q)
q σ
Fig. 14: Percus-Yevick static structure factor of hard spheres.
2 3 4 5
r/a
0 1 2 3 4 5 6
g(r)
MD PY-VW PY RY φ = 0.49
Fig. 15: Hard-sphere g(r) for a volume fraction Φ = 0.49 close to the freezing transition.
Comparison between PY, RY, Verlet-Weis corrected PY and MD computer simulations (filled circles). From Ref. [28].
h(r)−c(r) by Fourier-inverting γ(q) = ρc(q)h(q) = [S(q)−1]2/(ρS(q)). The hard-sphere g(r)follows then in PY approximation from g(r) = y(r) = 1 +γ(r)forr > σ. Contrary to g(r), γ(r) and its first two derivatives are continuous also at r = σ as one can deduce from the OZ equation. Notice here that the identityy(r) = 1 +γ(r)holds true only within PY ap-proximation. The relation betweeny(r) andγ(r)is approximated in the HNC approximation by y(r) = exp{γ(r)}, which agrees with the corresponding PY relation for smallγ(r) only.
While the full PY-g(r)of hard spheres can not be represented analytically, one can derive closed expressions for the contact values ofg(r)and its first derivative from the continuity ofy(r)and its derivative atr =σ:
g(r=σ+) = 1 + 0.5 Φ
(1−Φ)2 , (134)
and
σdg
dr(r=σ+) =−4.5Φ (1 + Φ)
(1−Φ)3 . (135)
PY results for the hard-sphere c(r) and g(r), and for the static structure factor, are shown in Figs. 13 and 14, respectively, for various volume fractions. ForΦ = 0.1, we further show the
cavity function withy(r) = −c(r) for r < σ. The PY approximation provides a quite good representation of the true hard-sphereS(q)andg(r)for volume fractionsΦ ≤ 0.35. At larger values of Φit underestimates the contact value of g(r) and its oscillations are slightly out of phase, as can be seen from Fig. 15 in comparison with Molecular Dynamics (MD) computer simulation results. The PY approximation further fails to predict the liquid-solid freezing tran-sition which occurs for hard spheres atΦf = 0.494. This failure is not restricted to the PY approximation: The Ornstein-Zernike integral equations discussed in this lecture have been de-signed to describe the homogeneous fluid state, and are thus not suited to describe the density inhomogeneities and the symmetry change in a first-order liquid-crystal phase transition.
On the basis of the analytic PY solution, Verlet and Weis [29] provide a simple prescription to obtain results for the hard-sphereg(r)and S(q), which are in good agreement with computer simulation results even up to the freezing volume fraction (see also [30]).
For a givenσandΦ, the Verlet-Weis correctedg(r)is given by gV W(x; Φ) =gP Y(xσ
σ′; Φ′) +A(Φ′)e−(µ−1)x
x cos [µ(x−1)] , (136) withx=r/σ, and a rescaled volume fractionΦ′ ≤Φand diameterσ′ = (Φ′/Φ)1/3σ. The three parametersA, µ, andΦ′ are determined on demanding that gV W(σ+) = (pCS/pid−1)/4Φ = (1−0.5 Φ)/(1−ϕ)3 andSV W(q = 0) = kBT (∂ρ/∂pCS)T in agreement with the Carnahan-Starling equation of state for the hard-sphere pressure given in Eq. (140), and from minimizing the integral over|gM D(r/σ; Φ)−gP Y(r/σ′; Φ′)|for an interval ranging fromr/σ = 1.6to3.0.
Here, gM D(r) is the ’exact’ rdf obtained from Molecular Dynamics (MD) simulations. This leads toΦ′ = Φ (1−Φ/16), and toAandµexpressed as simple functions ofΦ′ [29]. Fig. 15 compares the rdf predictions of various integral equation schemes, including the VW-corrected PY and RY schemes, for a hard-sphere system at volume concentration Φ = 0.49 close to freezing. The principal peak height ofS(q)at wave numberqm, predicted by the VW-corrected PY scheme, is well described by the expression [7, 31]
S(qm) = 1 + 0.644 Φ(1−0.5Φ)
(1−Φ)3 , (137)
and conforms with the empirical Hansen-Verlet freezing criterion [32]. According to this cri-terion, a three-dimensional fluid freezes into a solid when the principal peak height,S(qm), of S(q)is located in between 2.8 to 3.1. The precise height depends to some extent on the range of the pair potential, and amounts toS(qm; Φf = 0.494) = 2.85 in the case of hard spheres.
For dispersions of strongly repulsive colloidal particles, this criterion becomes equivalent to a dynamic freezing criterion [33], as has been shown on the basis of a dynamic mode coupling theory [34]. For an alternativeanalyticparametrization of the rdf of hard-sphere fluids see [35].