Some simple results for the properties of polar fluids

The author's lecture notes concerning the correlation functions and the thermodynamics of a simple polar fluid are summarized. The emphasis is on the dipolar hard sphere fluid and the mean spherical approximation and on the relation of these results to the Clausius-Mossotti and Onsager formulae for the dielectric constant. Previous excerpts from these lecture notes, Condens. Matter Phys., 2009, 12, 127; ibid., 2010, 13, 13002, have contained results that were not widely known. It is hoped that this third, and likely final, excerpt will prove equally helpful by gathering several results together and making these more widely available and recording a few new results.


Introduction
This paper is dedicated to Yura Kalyuzhnyi on the occasion of his sixtieth birthday; it is the result of the beautiful work of Michael Wertheim on hard spheres and dipolar hard spheres that has inspired the author and many others, including Yura and his colleagues in Lviv. There is little in this paper, drawn from the author's lecture notes, that is not well-known to Michael but, perhaps, lesser mortals will benefit from this collection of results and those in two previous papers [1,2] taken from the author's lecture notes.
Dipolar hard spheres are a simple representative molecular fluid and polar fluid. For a canonical ensemble, the well known definition for the h-particle correlation function for a simple atomistic fluid of N molecules in a volume V is easily generalized for a molecular fluid, where is the energy of the system, β = 1/kT (T is the temperature and k is the Boltzmann constant), Q N is the configurational partition function given by and R ij = |r i − r j | is the distance between the centers of a pair of molecules, i and j, whose centers are located at r i and r j . The "volume" elements, dΩ, are normalized so that dΩ = 1. Thus, The function u(R ij , Ω i , Ω j ), that can be written as u(ij) for brevity, is the intermolecular potential between a pair of molecules.
The following notation is employed. A function, such as g(ij), that depends upon the orientation of each member of a pair of particles, i and j, is denoted by the presence of each of the indices of the two particles in the argument of the function. After integration over the orientations of the two particles, a function that depends only on the scalar separation, R, of the two molecules results. This spherically averaged function is denoted by the subscript 0. Thus, the radial distribution function (RDF) is given by g 0 (R) = g(12)dΩ 1 dΩ 2 = g (12) .
The Ornstein-Zernike (OZ) equation becomes h(12) = c(12) + ρ h(13)c(23)dr 3 dΩ 3 , where ρ = N/V , h(12) = g(12) − 1 and c (12) are the total and direct correlation functions. The common thermodynamic functions are given by E = 1 2 N ρ g (12)u(12)dr 2 dΩ 1 Ω 2 = 1 2 N ρ g (12)u(12) dr 2 , pV N kT = 1 − 1 6 βρ g (12)u (12)R 12 dr 2 , and In the above E is the energy in excess of the kinetic energy terms, p is the pressure, and h 0 (12) = g 0 (12) − 1. The functions, h 0 (12) and g 0 (12), are the total and pair correlation functions, respectively. This article gives only an outline of the field of molecular fluids. The discussion will be restricted to molecular fluids with a hard core. For convenience, molecular fluids can be divided further into two broad types, (1) fluids in which the hard core is spherical (the asymmetry comes from the attractive potential) and (2) fluids in which even the hard core is nonspherical. The first class is conceptually simpler and is considered here. Dipolar hard spheres will be considered as an example of this class. Liquid crystals are an example of the second class and may, perhaps, be considered in a future installment.

Dipolar hard spheres
As an example of a molecular fluid, we consider the dipolar hard sphere fluid where the intermolecular potential is given by where µ and σ are the magnitude of the dipole moment and diameter of the dipolar hard spheres, whereR 12 = R 12 /|R 12 | is a unit vector in the direction of R 12 ,ê 1 is a unit vector in the direction of dipole 1, and ∆(12) =ê 1 ·ê 2 .
The function ∆(12) does not appear in the intermolecular potential, except as part of D (12). However, D(12) and ∆(12) contribute independently to the correlation functions. The dipoles are assumed to be nonpolarizable.

33001-2
Barker [3] has proved the very useful theorem that is given in the following two equations, Note that b could beê j , i = j. We shall call these results Barker's theorem. Using Barker's theorem, it is easy to show that 1, D (12), and ∆(12) are orthogonal, and The normalization of 1, D (12), and ∆ (12) can also be obtained from Barker's theorem and is and This means that 1, D (12), and ∆(12) are part of an orthogonal basis set. Indeed, they are a subset of the spherical harmonics. A basis set is a linearly independent set of functions with the property that any function can be expressed as a linear combination of the members of the basis set. One basis set for three-dimensional Euclidean space is the set of vectors in the directions of the x, y, and z axes. The space of all functions for which the spherical harmonics are the basis set has an infinite dimension. As will be seen, in some special cases the functions 1, D(12), and ∆(12) form a complete basis set of finite (three) dimension but this is not usually the case. A basis set need not consist of orthogonal vectors or functions. However, it is convenient if they are orthogonal. A nonorthogonal basis set can be transformed into an orthogonal basis set by what is called the Schmidt orthogonalization procedure. Hence, we can expand where g 0 (R 12 ) is given by equation (5), and h ∆ (R 12 ) = 3 ∆(12)g(12)dΩ 1 dΩ 2 .

33001-3
The common thermodynamic functions are given by and where y(12) is the background, or cavity, function, y(12) = exp[βu (12)]g (12) and u(12) is the pair interaction. Note that y(12) is a continuous function even if u(12) is discontinuous. The functions y 0 (R) and h 0 (R) are the spherically averaged projections of y(12) and h(12), respectively. As we shall see shortly, the dielectric constant is also given by an integral involving h ∆ . This means that the dielectric constant and common thermodynamic functions can be obtained from g 0 , h D , and h ∆ even if the other projections are not known. Of course, in general, to obtain these three projections, the other projections must be calculated. In any case, these three projections can be called the active projections for the dipolar hard sphere fluid since they determine the common thermodynamic functions and the dielectric constant of this fluid.

Simple treatments of the dielectric constant
The simple treatments considered here are based on the concept of the local field, E loc , felt by a dipole. This is not equal to the applied field, E, because of the other dipoles. Let us carve out a sphere of volume a, centered at a dipole. Since the dipole-dipole interaction is long ranged, we may assume that the dielectric or polar fluid is a continuum outside this sphere. Choose the volume of this sphere to be equal to the volume per dipole, The average value of µ is related to E loc , As E loc is relatively small, the exponentials in equation (27) may be linearized. Thus, The integrals of the first term in the numerator and the second term in the denominator vanish. The result is Our task is to calculate E loc . We will consider two simple approaches first.

Clausius-Mossotti result for
The field inside the dielectric fluid is different from the applied field due to the polarization of this fluid.
where D and P are the electric displacement and polarization vectors. For an isotropic system, the vectors have the same direction. Thus, Lorentz [4,5] argued that there are four contributions to E loc : (1) the applied field, E; (2) the volume charge contribution of P, which is zero because P is a constant and ∇ · P = 0; (3) the surface charge contribution of P on the surface of the sphere of radius a; and (4) the field due to the dipole, which is independent of E and so does not contribute to µ . For a surface element of this sphere at a polar angle, θ, measured from the direction of P and E, the area of this element is dS = 2πa 2 sin θdθ. The surface charge density in the direction of the normal to the surface of the sphere at the polar angle θ due to the polarization is P n = P cos θ. Thus, the element of the field due to the surface polarization is dE = P n dS/a 2 and E is and Recalling that

33001-5
which is the Clausius-Mossotti result [6,7]. This is not a very good result because the CM diverges when y = 1, for which there is no experimental support. Sometimes this problem is called the polarization catastrophe. The CM result for is plotted and compared with some simulation results [8] in figure 1. There is no singularity in the simulation results. The constant y is not to be confused with the background function, y(12).

Onsager result for
To obtain E loc , Onsager [9] solved the boundary value problem for a sphere of radius a and dielectric constant equal to unity within an infinite dielectric medium whose dielectric constant is and with an applied field E. Denote the potential inside and outside the sphere by φ 1 and φ 2 , respectively. Thus, The potential and displacement are continuous across the surface of the sphere so that φ 1 (a) = φ 2 (a) and ∂φ 1 (a)/∂R = ∂φ 2 (a)/∂R. The potential φ 1 is finite inside the sphere (in particular at R = 0) and, far from the sphere, φ 2 = −ER cos θ. The solution of this boundary value problem is Hence, Using, the dielectric constant is given by Explicitly, = [1 + 9y + 3 1 + 2y + 9y 2 ]/4. This is Onsager's formula. It is plotted in figure 1. This result does not diverge and is much better than the Clausius-Mossotti result. Until Wertheim's result, this was the standard formula. Wertheim's results will now be considered. However, some preliminary formulae are needed.

Fourier transform of h ∆ (R)
As has been mentioned, h(12) can be expanded in spherical harmonics, A similar expansion can be made for c (12). We will want to substitute these expressions into the OZ equation, equation (6). To do this it is convenient to use the Fourier transform. The Fourier transforms of h 0 , c 0 , h ∆ , and c ∆ are straightforward. However, the Fourier transforms of h D (12) and c D (12) are more complicated because D(12) contains R and we must transform the combinations h D (R)D (12) and c D (R)D(12) as wholes. First, recall that the Fourier transform pair is where R = |r| and k = |k|. Choose the coordinate system so that r = R(sin θ sin φ, sin θ cos φ, cos θ) and k = (0, 0, k). For an easier notation, define T (12) = h D (R)D (12). After some algebra, where andh where The functionf (k) is sometimes called a Hankel transform. Note thath D (k) is noth D (k), the Fourier transform of h D (R), which is given bỹ The functions j 0 (x) and j 2 (x) are spherical Bessel functions. Equation (49) is a perfectly good result for the Fourier transform of h D (R)D(12) but it is a nuisance to have two kinds of transforms. Thus, it is useful to define an auxiliary function, because, as is seen by straightforward integration, the Fourier transform of F (R) is the Hankel transform of f (R). Thus, in summary, the Fourier transform of T (12) = h D (R)D(12) is given bỹ where D k (12) is given by equation (50) and The inverse of the last equation is This auxiliary function has another interesting property. then This property can be exploited to evaluate integrals of long ranged functions which would be difficult if direct integration were attempted.

Fourier transform of the OZ equation for dipolar hard spheres
It has been seen that and Hence,h 1 To take the transform of the convolution in the OZ relation, we must evaluate integrals of the form ∆(13)D(23)dΩ 3 .
To do this, Wertheim's "multiplication table", which is given in with similar equations for the transforms of the higher order terms. We know that and

33001-8
Equations (73)-(75) are fine but we are interested in H D (R) and C D (R) rather than in h D (R) and c D (R). It is easy to show that where The parameter, K, is independent of R but depends on T , ρ, , etc., and is not known until the problem is solved. We can establish an interesting result forC D (0). We know that Using equations (75) and (78) it follows that From this, we have

Some exact results for
Onsager's expression is a special case of the exact result [10,11] ( − 1)(2 + 1) 9 = 4πβρ 9 where M is the total dipole moment of the dielectric. Write this as The parameter g K is called the Kirkwood g K factor. The g K factor can be written as an integral, yielding so that ( − 1)(2 + 1) The Onsager approximation consists in neglecting the contribution of h ∆ (R). Some other interesting exact results for can be obtained using the OZ equation given above. We can use the truncated versions of the expressions for h 0 , h ∆ , and h D , namely, The missing terms do not contribute. Solving the truncated equation (89) forh ∆ gives Solving equation (90) forH D gives, for k = 0, where x = 1 − 1 3 ρc ∆ (0) and equations (82) and (91) have been used. Equation (91) can be rewritten as Thus, Using equation (92) yields The solution of this equation can be verified to be Hence, the Clausius-Mossotti result is obtained by neglecting c ∆ (R).
These three routes to may not be consistent for a given approximation. However, they will be consistent if the OZ relation is satisfied. Note also, that we have obtained exact expressions forh ∆ (0),H D (0),c ∆ (0), andC D (0)!

The mean spherical approximation for the dipolar hard sphere fluid
Because the MSA is a linearized approximation, 1, ∆, and D are a complete basis set for the MSA. Thus, equations (88)-(90) can be employed. The MSA is
The first thing to note is that h 0 and c 0 are decoupled from equations (89) and (90). Equations (100) and (101) are the Percus-Yevick (PY) approximation for a hard sphere (HS) fluid. Thus, Algorithms for calculating h PY HS (R) = g PY HS (R) − 1, that are based on the formulae of Thiele [13] and Wertheim [14,15], have been given previously [16,17]. The other two equations may be solved by introducing the new functions, and After a little algebra, the decoupled equations, follow. The MSA approximation consists of and Hence, and where h PY HS (R; ρ) are the PY hard sphere correlation functions. The equations for c + and c − are similar. The contact values of h + (R; ρ) and h ( R; ρ) are given by and where η = πρσ 3 /6. The algorithms of Smith et al. [16,17] can be used. These algorithms are robust and, with a small change, give sensible results, even for the negative densities required by equation (116). The required change is that the one line in the program where a cube root of a quantity involving η

33001-11
is taken, the instruction should be changed so that when η is negative, the absolute value of η is used and the resulting cube root is multiplied by −1. Thus, and The contact values of h ∆ (R) and H D (R) follow from equations (117) and (118) together with The function that we want is h D (R), not H D (R). This can be calculated from equation (57). In particular, the contact value of h D (R) is The parameter K is not yet specified but we are in a position to do so now. Using equation (82) and yields which specifies K, which has been renormalized so that it is dimensionless. Note that 0 < Kη < 1/2. When Kη = 0, y = 0 and when Kη = 1/2, y = ∞. The correction functions for the dipolar hard sphere fluid that follow from the MSA are plotted and compared with simulation results [18,19] in figure 2 for a representative case. The MSA gives fairly accurate results for g 0 (R). The simulation results for g 0 (R) for dipolar hard spheres are very nearly equal to those for hard spheres but are slightly larger. Hence, one prediction of the MSA is that g 0 (R) for dipolar hard spheres is independent of the magnitude of the dipole moment and is equal to the radial distribution function of a hard sphere fluid. This prediction is not exact but is  quite well satisfied by the simulation results. However, the MSA results for h D (R) and h ∆ (R) are rather poor. Interestingly, the approximations, called LEXP, and are much better.

MSA thermodynamic functions
Using the compressibility route, the thermodynamic functions are as follows: When the compressibility equation is used, this is a very poor result since the MSA incorrectly predicts, that there is no contribution from the dipolar part of the intermolecular potential.
Using the pressure route, which becomes pV N kT and, using the energy route, The points are simulation results [19]. The dashed curves marked MSA, 2, and 2+3 give the results of the MSA and perturbation theory when truncated after 2 and 3 terms, respectively. The solid curve gives the results of the Padé extrapolation of Rushbrooke et al. [24].
One interesting characteristic of the MSA is that the energy can be integrated analytically to give the (energy equation) free energy and this free energy can be differentiated to yield the (energy equation) pressure. After a little algebra, the results are 33001-13 where A HS and p HS are the hard sphere free energy and pressure, respectively. Note that K has been renormalized so that it is dimensionless. The free energy that results from the MSA is plotted and compared with simulation results [20] in figure 3. As is usually the case, the thermodynamics obtained from the energy equation are more accurate than those obtained from the compressibility or pressure equations.

MSA dielectric constant
The dielectric constant can be calculated by the three routes given above. All three routes yield the same expression for . For example, starting with equation (97), we obtain Solving for gives Some MSA results for are plotted in figure 1. The agreement of the MSA result with the simulation results is better than for the CM and Onsager theories but the MSA results are still too small. Expanding the MSA expression for gives which is correct to order y 3 . By contrast, the CM result is and the Onsager result is, on expansion, There is no term of order y 3 in the CM theory. The Onsager coefficient of the y 3 term is too negative.
Of course, approximations that are better than the MSA approximation can be used. For example, Fries and Patey [21] used the hypernetted chain (HNC) approximation. The HNC results are better than the MSA results but the calculations are lengthy and, in contrast to the theories considered here, do not yield analytic results. It is to be noted that other combinations of β, µ, and ρ, besides y, appear when the HNC approximation is employed. This is true of other more general theories, for example the perturbation theory that is considered below.

Perturbation theory for dipolar hard spheres
Perturbation theory has been found to be very successful for simple fluids. It is natural to wonder if perturbation theory might also be useful for a polar fluid. The answer is yes but some qualifications are necessary.
By expanding the free energy in powers of β, the following result is obtained

33001-14
where A HS is the hard sphere free energy. The quantities A 2 and A 3 are given by and where g HS (R) and g HS (123) are the pair and triplet distribution functions of the hard sphere fluid, and where the θ i are the three interior angles of the triangle formed by the three sides, R ij . The angle θ 1 is the angle opposite the side R 23 , etc. Barker's theorem has been used to perform/simplify the angular integrations. The term of order β vanishes on angular integration, as do some of the terms of order β 2 and β 3 that formally contribute. Barker et al. [22] and Tani et al. [23] have calculated I ddd by simulation and direct integration via the superposition approximation, g HS (123) = g HS (12)g HS (13)g HS (23). A numerical fit of their results is given by where ρ * = ρσ 3 . As is seen in figure 3, this truncated series gives poor results. However the Padé, that was proposed by Rushbrooke et al. [24], gives excellent agreement with the simulation results. A Padé tends to work best for alternating series. For example, the series 1−1+1−1+· · · is summed correctly to 1/2 by a Padé. Patey and Valleau [20] refer to the Padé results as "absurdly successful". This is meant as a positive comment and is a fair observation. Unfortunately, a Padé does not work well for the correlation functions. Some thoughts about the development of approximations that are consistent with equations (143) have been considered by Barker and Henderson [25]. However, nothing much has come of these efforts. The dielectric constant can be calculated from ( − 1)(2 + 1) 9 = y 1 + 9I dd∆ 16π 2 + · · · (144) where [22,23] I dd∆ = 3 cos 2 θ 3 − 1 (R 13 R 23 ) 3 g HS (123)dr 2 dr 3 = 17π 2 9 σ 6 1 − 0.93952ρ * + 0.36714ρ * 2 1 − 0.92398ρ * + 0.23323ρ * 2 .
This gives poor results, even with a Padé. However, the direct expansion, due to Tani et al.
[23], = 1 + 3y + 3y 2 + 3y 3 9I dd∆ 16π 2 − 1 , gives very good results, as is in figure 1. Perturbation theory can be recast by subtracting the MSA contributions from the perturbation terms and writing the perturbation theory as a series of corrections to the MSA. This was done by Henderson et al. [26]. The results are similar to those of the perturbation theory considered here.

A few remarks
The mean spherical approximation and perturbation theory are pleasing extensions of the classic theories of Clausius and Mossotti and Onsager for polar fluids. Not only do they make more accurate predictions for the dielectric constant of a polar fluid but they predict the other thermodynamic properties of polar fluids. Many of these ideas are applicable to polar fluids with a dispersion interaction. The simplest system of the kind is the dipolar Yukawa fluid. Szalai et al. [27] and Mate et al. [28] have considered this model to be polar fluid.
The molecules considered here are unpolarizable dipoles. Onsager considered polarizable dipoles. Valiskó et al. have generalized some of the expressions presented here for polarizable dipoles and have made simulations for a polarizable dipolar hard sphere fluid.
The author had hoped to include his lecture notes on liquid crystals as an example of a molecular fluid with a nonspherical hard core. However, despite some searching, these have not been found. If they do come to light, they can form a fourth part of this series and the "concerto" with three movements can become a "symphony" with four movements.