Langmuir 1991, 7, 1383-1387
Ordering in Charged Colloidal Dispersions and Macroionic Solutions: A Density-Functional Approximation Paul Salgi and Raj Rajagopalan’ Department of Chemical Engineering, University of Houston, Houston, Texas 77204-4792 Received October 19,1990. I n Final Form: February 11,1991 We examine a density-functional approach for constructing the phase diagram of charged colloidal dispersionsinteracting through screened Coulombic (Yukawa)potentials. The interaction potential used accounts for the finite size of the colloids (or, “macroions“). The theoretical calculations predict the three phases (body-centered cubic (bcc), face-centered cubic (fcc), and liquidlike) that have been observed experimentally. The bcc structure, observed in the region of low electrolyteconcentrationsand relatively low volume fractions of the macroions, becomes unstable at higher electrolyte concentrations, and the theory predicts a bcc-to-fcc or liquid-to-fcctransition. The predictions are qualitatively consistent with the recently reported experimental results based on small-angle synchrotron X-ray scattering.
Introduction Order/disorder transitions in dispersions of monosized, charged colloidal spheres have attracted considerable attention in the last decade since these dispersions mimic a number of phenomena of interest in atomic systems. Since methods are available for controlling the size,charge, and uniformity (i.e., sphericity as well as uniformity of size and charge) of the colloidal particles (usually polystyrene or other latex particles or silica particles), such dispersions serve as extremely attractive model manybody systems for studying freezing and melting transitions and mechanisms, a t least for spherically symmetric interactions. In the case of neutrally charged (usually polymer-coated, sterically stabilized)spheres, one observes the well-known Kirkwood-Alder hard-sphere phase transition. The other extreme, namely, highly charged spheres dispersed in electrolyte solutions, presents a much more interesting case. The dispersions in this case can be thought of as highly asymmetric electrolytes, with the colloidalparticles serving as “macroions” in a sea of regular ions. The effective (mean-field) interaction potential between the colloidal particles under these conditions can be represented sufficiently accurately by screened coulombic potentials (Yukawa potentials), whose strength and range can be controlled experimentally. Therefore, the study of the structure and dynamics of such dispersions provides a convenient way to investigate phase transitions in systems interacting through soft,longranged potentials. In fact, because the motion of the individual particles and the phase structure of the dispersions can be examined directly by using optical methods, it has been possible to combine digital imaging techniques and optical observations to follow in detail the melting transitions and the topological defects, as has been demonstrated by Murray and co-workers and others for two-dimensional sy5tems.l The closest system that has been examined systematically (theoretically) in the literature is a collection of pointlike Yukawa particles. Robbins et aL2have presented the phase diagram for this (effectively two-parameter) system using a combination of molecular dynamics and lattice dynamics methods and application of empirical
* To whom correspondence should be addreesed.
(1) Murray, C. A.; Wenk, R. A. Phys. Rev. Lett. 1989,62,1643. Murray, C. A.; Van Winkle, D. H.; Wenk, R. A. Phose Transitions 1990,21, 93. See also Tang, Y.; Armstrong, A. J.; Mockler, R. C.; OSullivan, W. J. Phose Transitions 1990,21, 75. (2) Robbins, M. 0.;Kremer, K.; Greet, G. S . J. Chem. Phys. 1988,88, 3286.
Lindemann and Hansen-Verlet criteria for melting and freezing, respectively. The liquid-bcc-fcc transitions predicted by Robbins et ala2seem to be consistent with the experimentalresults for (polystyrenelatex) dispersions at low densities and large screening lengths reported by Monovoukas and Gast,3 after the charge of the particles is renormalized4 to a lower value. However, in the case of dispersions, the finite size of the particles begins to have a significant effect a t low charges, and it is therefore important to examine ways to include the effects of the core in the calculations. This point, in fact, has been also emphasized by Rosenberg and Thir~malai.~ The substantial influence of the core may be seen from the experimental results of Sirota et for charged polystyrene dispersions. These authors employed highresolution, small-angle, synchrotron X-ray scattering (in order to avoid potential problems due to multiple scattering encountered in light scattering) and probed the phase structure over a fairly wide range of particle concentrations. Their results demonstrate that, in addition to the formation of face-centered cubic (fcc), bodycentered cubic (bcc), and liquidlike phases (aspreviously established”),glassy phases appear even in monodispersed systems, at lower particle number densities than would be expected for pointlike particles. The purpose of the present paper is to examine if an approach based on density-functional theories could be used to study phase transitions in colloids. We consider colloids interacting through repulsive Yukawa potentials with a hard core and present a procedure based on a density-functional formulation that incorporatesthe effect of the repulsion as a perturbation from hard-sphere interactions. Density-functional theories of freezing have recently been found to be very successful in a few cases. For instance, various formulations of these theories“l2 predict the phase transition parameters of a hard-sphere (3) Monovoukas, Y.; Gaet, A. P. J. Colloid Znterface Sci. 1989,128, 533. (4) Alexander, S.; Chaikin, P. M.; Grant, P.;Morales, G. J.; Pincus, P.; Hone, D. J. Chem. Phys. 1984,80, 5776. (5) Roeenberg, R. 0.;Thirumalai, D. Phys. Rev. A: Gen. Phys. 1987, 36,5690. H. D.; Sinha, S. K.; Chaikin, P. M.; Axe, (6) Sirota, E. B.; Ou-Yang, J. D.; Fujii, Y. Phys. Reo. Lett. 1989, 62, 1524. (7) See, for a review of the literature on ordering in diepersions, Hirt-
zel,C. S.;Rajagopalan,R. ColloidalPhenomeno: Adoonced Topics;Noyea Publishers: Park Ridge, NJ, 1985. (8) Baus, M.; Colot, J. L. Mol. Phys. 1988, 55, 653. (9) Lutako, J. F.; Baus, M. Phys. Rev. Lett. 1990,64, 761.
0 1991 American Chemical Society
Salgi and Rajagopalan
1384 Langmuir, Vol. 7, No. 7, 1991
fluid very accurately. Curtin and Ashcroft13have further extended the above formalism to construct the phase diagram of Lennard-Jones fluids. In the present paper we show that the density-functional theory formalism can be used to predict all the basic features of equilibrium phase transitions observed in charged dispersions. The theoretical results are compared with the recent experimental data6 on dispersions at high volume fractions and the similarities and differences are examined. Theoretical Formulation We take the interactions between macroions to have the following Yukawa form
20.0
I-
3
m
3
.-(d
L
U
-5.0
2
1
0
3
4
5
r/r
where r is the interparticle center-to-center distance, u is the diameter of the particles, $0 is the surface potential, t is the dielectric constant of the solvent medium, and co is the permittivity of free space. Here, K is the usual DebyeHiickel inverse screening length, determined by the ionic strength of the solution. The surface potential $ais related to the electronic charge 2 on the macroion to a good approximation by14 $o = Ze/[aeota(2
+ KU)]
(2)
The above Yukawa form (eq 1) followsfrom the free-energy formulation suggested by Verwey and Overbeek,14where the electrostatic potential in the electrical double layers surrounding the particles is obtained from the solution of the linearized Poisson-Boltzmann equation (inthe DebyeHiickel approximation). In this sense, eq 1represents the effective (mean-field) interaction between the colloidal particles15 dispersed in an electrolyte solution, which is theoretically represented as an effective dielectric. Even though the Debye-Hiickel approximation breaks down in the region of strong potential gradient close to the surface of each sphere, it can still be used if the charge of the colloidal particle is suitably ren~rmalized.~ We use the above interaction potential to construct the interparticle correlation functions (in homogeneous colloidal fluids) needed in the density-functional formalism. Figure 1 presents an example of the pair-potential used here. The major difference between this and the potentials used in Robbins et a1.2and Rosenberg and Thirumalai5is that we include the hard-sphere cutoff in eq 1. This difference becomes important at high effective temperatures and at large particle concentrations. As done by Curtin and Ashcroft13 for Lennard-Jones fluids, we use a hard-sphere perturbation approximation to separate the Helmholtz free energy of the colloid crystal into primarily entropic and internal energy contributions. The entropy of the hard-sphere reference system, however, is obtained in the modified weighted-density approximation (MWDA)" instead of WDA,'OJ3 and the internal energy is determined from an expansion about the uniform liquid value. In addition to yielding results of comparable accuracy in predicting the hard-sphere free energy,MWDA (10) Curtin, W. A.; Ashcroft, N. W. Phys. Rev. A: Gen. Phys. 1985,32, 2909. (11) Denton, A. R.; Ashcroft, N. W. Phys. Rev. A: Gen. Phys. 1989, 39,4701. (12) Tarazona, P. Mol. Phys. 1984.52, 81. (13) Curtin, W. A.; Ashcroft, N. W. Phys. Rev. Lett. 1986,56, 2775. (14) Verwey, E. J. W.; Overbeek, J. Th.G. Theory of the Stability of Lyophobic Colloids; Elsevier: Amsterdam, 1948. (15) Chen, S.-H.; Sheu, E. Y. In Micellar Solutions and Microemulsions: Structure, Dynamics and Statistical Thermodynamics,Chen, S.H.,Rajagopalan, R., Ede.; Springer-Verlag: New York, 1990; Chapter 1.
Figure 1. Example of the Yukawa potential used in the theory. The diameter of the particle is u, and the contact potential is UO.
has the advantage of being computationally much simpler than WDA. In what follows, we will briefly outline the derivation of the free energy functional used here and the interested reader is referred to ref 13 for more details. The excess free energy per particle for an inhomogeneous fluid is first expanded around a reference hard-sphere system of diameter u (the actual diameter) and may thus be written as
fb1 fHSIP;al + (1/2N)
1dr dr' u(r' - r)pHs(2)[rt,r;p,u1(3)
where u is the full pair-potential. Here,fw[p;u] andpHs(2). [r',r;p,u] are the excess free energy per particle and pair distribution function for the hard-sphere (HS) system of diameter u and are unknown quantities. The second term in eq 3 is approximated in the following way. We first expand the exact f[p] about that of the liquid to all orders ",
f[~= 1 fq(po) - Z ( l / N B n ! )
1dr, ... drlc(")(r, ...rl)
n-2
AdrJ
... Ap(rl)
X
(4)
The dn)are the nth-order direct correlation functions of the liquid, which can be written as dn)= [cHs(")+ Adn)], defining Adn) as the difference between the actual of the full system and that of the HS reference system. The quantity p(r) is the single-particle density, PO is the average number density, and Ap(r) is the difference between the two. We then use the homogeneous limit of eq 3 for the liquid free energy, sum the hard-sphere contributions to all orders, approximate the HS reference energy by using the MWDA, and neglect the higher order terms Ad") (n > 2). Following this sequence of steps, one easily obtains the excess free energy per particle of the inhomogeneous fluid as
2BN
d r dr' Ad2)(r - r') Ap(r) Ap(r') ( 5 )
where PHSL(~) is the distribution function of the hardsphere liquid of diameter u. Equation 5 is used to compute the excess free energies of the crystalline and liquid phases. For the liquid phase, the first term on the right-hand side is obtained by using the Carnahan-Starling equation's
Langmuir, Vol. 7, No. 7, 1991 1386
Ordering in Charged Colloidal Dispersions
and the third term vanishes. The second term can be evaluated analytically for Yukawa potentials.17 For the crystalline phase the third term is most easily evaluated in Fourier space. In order to do this, we assume the singleparticle density distribution in the crystalline phase to be a sum of identical Gaussians centered on the given lattice sites R (either bcc or fcc), namely p(r) =
( 4 ~ ) exp[-a(r ~ ’ ~ - RI2]
(6)
R
where a determines the Gaussian width and is the only parameter that is varied to minimize the free energy of the crystals. The case a = 0 corresponds to a uniform liquid, whereas a02 >> 1represents sharp, nonoverlapping Gaussian peaks. Additional comments concerning the use of the Gaussian approximation for the crystal density are presented in the next section. As mentioned earlier, for the crystalline phase, the first term in eq 5 is obtained by using the MWDA. In this approximation,ll the excess free energy per particle of the crystalline phase is given by f€is[P;rJl
fO(h
0.28 0
100
200
300
400
500
(7)
with
-5
0.28
Liquid
0
where c~s(~)(G;;)denotes the Fourier component of the uniform hard-sphere two-partide direct correlation!unction at the reciprocal lattice vector G and density p and fo is the excess free energy of the uqiform hard-sphere fluid. This is an implicit equation in p which is solved hy using standard iteration techniques. The quantity f’o(p) denotes the firsiderivative of fo with respect to the density p evaluated at p. Even though the Carnahan-Starling16 equation is the most accurate for determining fo, we use, for consistency, the Percus-Yevick (PY)approxim-ation16 for both fo and CHS(~),since the weighted density p is low enough for the PY equation of state to be accurate.13 The Yukawa direct correlation function is obtained by using the analytical solution of the Ornstein-Zernike equation in the mean-spherical approximation (MSA)due to Hayter and Penfoldlg with the rescaling correction (RMSAP introduced when needed. To the excess free energies of both the liquid and the solid, one must add the ideal-gas contribution
where A is the deBroglie thermal wavelength. The second term in the above equation vanishes for the liquid phase.
Results The phase diagram is constructed by comparing the free energies of the three phases, namely, liquid, fcc, and bcc (the only three equilibrium phases that have been observed experimentally), for each set of the relevant parameters. Figure 2 shows the phase diagram (volume fraction of the colloid versus electrolyte concentration at constant charge) obtained from the present study for typical experimental values of the different parameters. First, we note that the computed phase diagrams lead to (16) Carnahan, N. F.; Starling, N. E. J. Chem. Phys. 1969,51, 635. (17) Foiles, S.M.; Ashcroft, N. W. J. Chem. Phys. 1981, 75, 3694. (18) See, for example,Haneen, J. P.;McDonald, I. R. Theory of Simple Liquids, 2nd ed.; Academic Press: New York, 1986. (19)Hayter,J. B.; Penfold, J. Mol. Phys. 1981, 42, 109. (20) Hansen, J. P.; Hayter, J. B. Mol. Phys. 1982, 46,651.
>
1 -I 1 ’
0.26
c 0.24
0
100
200
300
400
500
Figure 2. Calculated phase diagram as a function of the volume fraction of the colloid and the electrolyte concentration,C H ~for , two values of the effective charge 2 (a) Z = 175; (b) Z = 250. The parameters used correspond to those in ref 6, i.e., particles of diameter 0.091 pm in 0.90.1 methanol/water suspension with a relative dielectric constant of 38.
physically consistent results. The results show bcc structure in the region of low volume fractions and low electrolyte concentrations, and as the volume fraction is increased a bcc-to-fcc transition is also predicted. At higher electrolyte concentrations the bcc structure is no longer stable and a liquid-to-fcc transition is observed. Moreover, as the charge of the particles is increased, the region of stability of the bcc structure is increased, which is in agreement with the results of Shih et al.,2l who found that for particles with low surface charges only the fcc phase is observed while the bcc phase forms at low volume fractions for systems with higher charges, with the fcc phase stable at higher volume fractions. For comparison, we have reproduced the experimental phase diagram of Sirota et al.! mentioned in the introductory section of this paper, in Figure 3. An examination of this figure reveals that the experimental results are in qualitative agreement with our calculations, although the volume fractions at which the bcc phase is observed in the experiments are lower than what is predicted by our calculations. A direct comparison of the theoretical predictionswith the experimentsis, however, difficult since the charges of the particles used in these experiments are not available sufficiently accurately. In fact, determination of the effective charges of the particles in a dispersion is not straightforward and cannot always be done in an (21) Shih, W. Y.; Aksay, I. A.; Kikuchi, R. J. J.Chem. Phys.1987,86, 5127.
1386 Langmuir, Vol. 7, No. 7,1991 35
30
(
t
I
'
I
I
Salgi and Rajagopalan ,
I
,
(
,
I
Char$ed PdlystyrLne Latex in 90% methanol +IO% water (vol.)
1
f
0.03 r I
=P
0.02
0.01 I . /
1 0 1
0
0.00 0
'
i 0
'
l
'
i 200
'
l
'
l
Sirota et ul 1990 l l ' l ' l '
400
600
0 800
Concentration HCI (micromoles)
Figure 3. Experimental phase diagram of ref 6. The lines are drawn to guide the eye and have no theoretical basis.
unambiguous manner. The effective charge Zis a function of both the macroion volume fraction and the electrolyte concentration, and it is usually estimated, using suitable models, from measurements of other properties of the dispersions. Sirota et al. have used shear modulus measurements on colloid crystals to estimate the charges of the dispersions used in their experiments, and the estimated charges vary from roughly 65 e to 135 e. Our phase diagrams in Figure 2 have been, however, calculated for fixed charges, since our goal here is to present merely a parametric study of the diagram. One may note, however, that as the effective charge Z is increased from 175 to 250, the agreement between the predictions and the experiments improves. The differences between the experimental data and the predictions can be improved further if higher effective charges can be justified. In this respect, since measurements such as shear moduli give information about particle interactions, one might ask whether it is possible to use a single effective charge Z to fit both the experimental phase diagram and the measured shear moduli, as in ref 2. Since only one experimental measurement of the shear modulus (at 4 = 0.12 and CHC~ = 0) is reported in ref 6, we have not attempted any detailed calculations of the shear modulus. In addition, the RMSA procedure used here for obtaining the direct correlation function for the liquidlike phase fails to converge for higher values of Z when the inverse screening length k (=xu) is low. We return to this point in the concluding part of this section. Finally, we note that the phase diagrams shown in Figure 2 have a region where the melting line appears to be horizontal. Theselines are not exactly horizontal but have a very small slope, which is not visible on the scale of the figures. The reason for the presence of the inflection point in the melting line is, however, not clear. WhileFigure 2 is convenientfor examining thestructure of the dispersion as a function of volume fraction, the effect of the range and the %oftness" of the repulsion on the structure is best seen when the phase diagram is plotted differently. In Figure 4 we present the results of our calculations for a fixed volume fraction in terms of the reciprocal contact potential U0-I(in multiples of ~ B Tat) phase transition as a function of the dimensionless inverse screening length k. (This representation of the phase diagram is analogous to the phase diagram presented by Robbins et a1.,2but ours is a three-parameter system, with the volume fraction 4 serving as the third parameter.) One can see from Figure 4 that, upon increasing the value
2
4
6
8
10
k Figure 4. Phase diagram in the Uo-1 vemw k plane for a fixed volume fraction of 0.31.
of the contact potential and thus the strength of the electrostatic interactions, the system undergoes a firstorder freezing transition from a liquidlike state to a bcc crystal for large screening lengths (small k) and to an fcc crystal for low screening lengths (large k). For a fixed value of the contact potential (with UOless than about 30),the dispersion assumes a liquidlike structure in the absence of any added electrolyte and crystallizes into a bcc crystal as the range of repulsion is reduced; on further reduction in k one observes a structural transition to an fcc phase. In the limit of k going to infinity, the system behaves like a neutral hard-sphere system and should then return to a liquidlike structure at the volume fraction considered in Figure 4. Although not shown in the figure, this return to a "liquid" phase can be, in fact, predicted with our model, which reduces in this limit to the modified weighted-density approximation (MWDA). As it has been already established, MWDA predicts a liquidlike stable structure for the neutral hard-sphere system a t this volume fraction." Finally, we note that the extension of the bcc/ fcc line is not shown for U0-l < 0.015in Figure 4 since the contact potential becomes too high for RMSA to converge. We note here that the Gaussian approximation for the crystal density may have an influence on the prediction of the bcc/fcc line. This Gaussian parametrization of the crystal density is based on the factthat the potential energy of a particle about its crystal lattice site is approximately harmonic. A more general way to parametrize p(r) is the usual Fourier expansion. The Fourier method can address possible anharmonicities and anisotropies of the crystal density. It has been shown, in the context of a density functional study of freezing (Laird et that the Fourier expansion and Gaussian approximation yield very similar phase diagrams for hard spheres and LennardJones systems, although there exist noticeable differences in the crystal densities predicted by the theory. This means that despite the high anharmonicity of the hardsphere fluid, density-functional approaches using the Gaussian approximation yield very reasonable transition parameters. The Gaussian approximation has also been used successfully in a study of the freezing of onecomponent plasmas (BarratZ3),which freeze into a bcc structure. However, lattice dynamics and molecular dynamics studies2J have shown that anharmonic effects are important in the study of the bcc-fcc transition in pointlike Yukawa systems. In fact, a simple Einstein (22) Laird, B. B.;McCoy, J. D.;Haymet, A. D.J. J. Chem. Phys. 1987, 87, 5449. (23)Barrat, J. L. Europhys. Lett. 1987, 3, 523.
Ordering in Charged Colloidal Dispersions
approximation leads to a temperature-independent b c e (Shih fcc transition at X = 1.72,where Xis equal to et and Hone et al.%). It is therefore possible that the Gaussian approximation may represent a potential difficulty when it comes to predicting the bce-fcc transition line using our approach. Nevertheless, we choose the Gaussian approximation here, as usually done, in view of its simplicity. The differences between the theory and experiments in the region of low electrolyte concentrations (where the repulsive tail is soft) can be attributed to the fact that the method used here is based on a hard-sphere perturbation technique which in principle becomes inadequate when the potential becomes long-ranged.26 It is interesting that the theory nevertheless correctly predicts the existence of bcc phase in the region of low screening; this, we believe, is in part due to the fact that the MWDA accounts, to some extent, for the multibody effects. As pointed out by Barrat,23these multibody contributions must be included to at least some extent to predict a bcc structure using a DFT approach. Another reason for the lack of quantitative accuracy may be the use of RMSA for the second-order direct correlation function. We have attempted to use the reference hypemetted-chain (RHNC) integralequation method instead of RMSA. However, because of the relatively high volume fractions considered, the RHNC procedure also did not converge in a number of cases and thus did not prove to be very useful for a systematic scan k). An alternative procedure of the parameter space (UO-', to deal with the low-electrolyte-concentrationand lowvoiume-fraction region would be to use a perturbation method with the charged hard-sphere model as reference system and an approach similar to that of Barrata to obtain the free-energy functional of the reference system. Finally,we note here a previous attempt at using densityfunctional theory to study freezing in colloidaldispersions. (24) Hone,D.; Alexander,&; Chaikin, P. M.; Pincus,P. J. Chem. Phys. 1983, 79, 1474. (25) One may note that the uniform limit of eq 3, which we use to
evaluate the free energy of the liquidlike dispersion, may not necessarily be accurate in the region of low electrolyte concentrations and low volume fractions because of the softness and the long range of the potential under such conditions. Nevertheless, this has little effect on the location of the phase transition since the same functional is used for both the uniform and nonuniform phases and is subtraded out subsequently.
Langmuir, Vol. 7, No. 7,1991 1387
Alexander et al.26 have used the Ramakrishnan-Yussouff theory?' reformulated in a slightly different way, to study the freezingtransition of pointlike Yukawa systems. Their approach,% which is basically a second-order theory, neglects higher-order effects and uses a two-orderparameter expansion of the effectivepotential. The results of Alexander et al. seem to predict that pointlike Yukawa systems freeze into a bcc structure for A values up to about 6. It is interesting to note that this is the only secondorder DFT approach (that we are aware of) which favors the bcc structure over the fcc structure. Several DFT studiesz3@have shown the necessity of including higherorder correlations to stabilize the bcc phase.
Concluding Remarks We have shown that an approach based on perturbationbased density-functionaltheory predicts the crystallization of a highly charged dispersion into body-centered cubic and face-centered cubic phases in the appropriate ranges of the screening length. More accurate predictions in the region of large charges and low electrolyte concentrations may be possible if sufficiently accurate liquid-state theories are available for calculating the free energies and direct correlationfunctions under similar conditions (ashas been noted in another context elsewhere).30
Acknowledgment. We thank the National Science Foundation and the Petroleum Research Fund, administered by the American Chemical Society, for partial support of the work presented here. We also thank the NSF Pittsburgh Supercomputing Center for computational time. (26) Alexander, S.; Chaikin, P. M.; Hone, D.; Pincue, P. A.; Schaefer, D. W.Phys. Chem. Liq. 1988,18, 207. (27) Ramakrishnan,T.V.;Yusaouff,M. Phy8.Reu.B: CondeM.Matter 1979,19,2776. (28) A number of differences exists between this work and o w . In contrast to the pointlike particles used in Alexander et al.,a we consider
the hard core in the potential explicitly in the present work. The freeenergy functional ueed is,therefore, a eecond-order expansion about the hard-sphere reference syetem. However, higher-order correlatiom in the reference system are incorporated in an approximate way through the use of MWDA, asmentioned earlier. Other differencesin the formulation of Alexander et al.salso exist, for which we refer the reader to the original paper. (29) Curtin, W. A. J. Chem. Phys. 1988,88, 7050. (30) Seeref 14 andHirtze1, C. S.;Rajagopalan,R. 1nMicekarSohtion.a a n d Microemulsions: Structure, Dynamics, a n d Statistical Springer-Verlag: ; Thermodynamics; Chen, S.-H., Rajagopalan, R., as. New York, 1990; Chapter 7.
