Langmuir 1992,8, 2898-2906
2898
Probing Interaction Forces in Colloidal Fluids through Static Structural Data: The Inverse Problem Raj Rajagopalan Department of Chemical Engineering, University of Houston, Houston, Texas 77204-4792 Received June 23, 1992. In Final Form: August 17, 1992
Determination of the details and the nature of the interaction forces in dense or strongly interacting dispersions has been approached in recent years using radiationscatteringexperiments and modem theories of liquid-state physics. Fitting the experimental structure factors to calculated ones based on assumed potentials, however, leads to ambiguous results. We present a perturbation-theory-based predictorcorrector method that recovers the interaction potentials from static structure factors with remarkable accuracy. The versatility and accuracy of the method are illustrated using three classes of potentials (one with a harshly repulsive core, one with a soft core, and a DLVO potential with primary and secondary minima) with theoretically generated structure factors. The method leads to analytical, closed-form expressions for the potential in terms of structure factors, is easy to apply, and enhances the usefulness of radiation scatteringdata. The practical issues relevant to the applicationof the method to experimental data are also outlined. Introduction The importance of understanding the interaction forces in colloidal dispersions has been recognized for a number of years, as evident from the attention given to this subject in numerous textbooks and monographs dating back to the early work of Verwey and O~erbeek.l-~ This subject has assumed a new significance in recent years with the emergence of the realization5 that the structure and properties of colloids (particularly at large concentrations or when they are strongly interacting) can be predicted and used profitably with the help of current advances in liquid-state p h y s i ~ s . ~Central g to this realization is the treatment of colloids as macroscopic analogs of atomic liquids, often emphasized by the use of the terms “supramolecular fluids”and “macrofluids”. This development is especially important since it also implies that colloids can be used as excellentmodel many-body systems to study phenomena in condensed matter physicsthat are otherwise inaccessibleexperimentallyat convenient length and time scales.lOJ1 Partly due to the above developments and partly motivated by the important role played by colloidal and “surface” forces in numerous industrial processes and (1)Verwey, E. J. W.; Overbeek, J. Th. G. Theory of the Stability of Lyophobic Colloids; Elsevier: Amsterdam, 1948. (2)Hunter, R. J. Foundations of Colloid Science; Oxford University Press: Oxford, UK, 1987,1989; Vols. I and 11. (3) Hiemenz, P. C.; Rajagopalan, R. Principles of Colloid and Surface Chemistry, 3rd ed.; Marcel Dekker: New York, 1993. (4)Israelachvili, J. N. Intermolecular & Surface Forces, 2nd ed.; Academic: London, UK, 1992. (5)Viewing colloids as ‘macrofluids” is already implicit in the work of van’t Hoff, Einstein, and Perrin [see the English translation of Perrin’s Les Atomes (Perrin, J. Atoms, Ox Bow Press: Woodbridge, CT, 1990)l. Although the concept of treating colloids as macrofluids has become popular and very useful since about 1975,the credit for applying a rigorous statisticalthermodynamic formulation to colloids for the first time appears to belong to Kirkwood and Mazur (Kirkwood, J. G.; Mazur, J. J. Polym. Sci. 1952, 9, 519). See, also, ref 6. (6)Hirtzel, C. S.; Rajagopalan, R. Aduanced Topics in Colloidal Phenomena, Noyes Publishers: Park Ridge, NJ, 1985. (7)Safran, S.A,, Clark, N. A., Eds. Physics of Complex and Supermolecular Fluids; Wiley: New York, 1987. (8) Ise, N., Sogami, I., Eds. Ordering and Organisation in Ionic Solutions, World Scientific: Singapore, 1988. Rajagopalan, R., Eds. Micellar Solutions and (9)Chen, S.-H., Microemulsions: Structure, Dynamics and Statistical Thermodynamics; Springer-Verlag: New York, 1990. (10) Murray, C. A.; Wenk, R. A. Phys. Rev. Lett. 1989, 62, 1643. (11)See the examples in ref 6 through 9 and the Preface to ref 9.
interfacial phenomena, numerous experimental innovations have also been advanced over the years for determining colloidal forces directly or indirectly. Notable among these is the so-called“surfaceforce apparatus” (also known as ‘Jacob’s Box”), which can measure directly the force of interaction between two coated or uncoated mica surfaces immersed in liquid^.^ Although this technique corresponds to measuring the forces between two flat surfaces (i.e., not between two particles), it has contributed to an enhanced understanding of electrostatic and steric forces in colloids as well as to the interpretation of structural forces4and the potentials of mean force between surfaces in the presence of added polymer or micellar entities.12J3 The measured forces have also been shown to be consistent with second viral coefficients obtained from light scattering measurements, at least in some cases.14 It appears that a direct measurement of forces between a colloidal particle and a surface may also be possible through the use of total internal reflection microscopy and evanescent wave scattering measurements.15 An important point to note about the above techniques, however, is that they measure forces between two objects (in a medium); i.e., in the context of colloids, the measurements correspond to forces between particles at extreme dilution. Therefore, it is important to explore methods that can provide pair-interaction forces (at least effectivepair-interaction forces) from measurements that are obtained from dense (i.e., concentrated) or strongly interacting systems. An important case of the latter is offered by micellar systems (either ionic or nonionic); the micelles, being thermodynamic entities, cease to exist upon dilution. The measurement of effective (or, ‘dressed”) (12)Richetti, P.; KBkicheff, P. Phys. Reu. Lett. 1992, 68, 1951. (13)Parker, J. L.; Richetti, P.; KBkicheff, P.; Sarman, S. Phys. Reo. Lett. 1992, 68, 1955. (14)Gee, M. L.; Tong, P.; Israelachvili, J. N.; Witten, T. A. J. Chem. Phys. 1990, 93,6057. (15)Prieve, D.C.; Luo, F.; Lanni, F. Faraday Discuss. Chem. SOC. 1987, 83, 22. Instead of using total internal reflection microscopy, evanescent wave Scattering from a large number of particles near a surface canalsobe meaauredusingstaticlightscatteringtechniques(Schumacher, G. A.; van de Ven, T. G. M. Langmuir 1991, 7, 2028) in order to get statistically better estimates for the distribution of distances of the particles from the surface. This information can also be used to deduce the particle/surface interaction potentials, but this procedure requires the assumption of a potential of interaction in advance of the calculations and, therefore, cannot be considered a “direct” method.
0 1992 American Chemical Society
Interaction Forces in Colloidal Fluids
chargesand interaction forces in such systemsthus requires methodsthat can accountfor (and "unscramble"the effects of) collective interactions implicit in the measured properties which are used as the basis for inferring the interaction forces. The objective of this paper is to propose such a method and to demonstrate its performance. For this, we use the so-called static structure factors obtainable from radiation scattering techniques (such as X-ray, neutron, or light scattering) as the source of information about the local structure of a given dispersion. (The definition of the structure factor and its relation to related structural quantities and interaction forces are given in a subsequent section.) We show that it is possible to obtain highly accurate information about the interaction forces in dense colloidalsystems from static structure factors. No a priori information about the nature of the functional form of the interaction potentials is needed, and the method is analytically and numerically straightforward. Also important is the fact that the method described here enhances the usefulness of scattering techniques, which have been shown in recent years to be very powerfuland indispensable in the study of colloidal systems. The above problem falls under a broad category known as inverse problems and carries with it a number of related issues such as whether the problem is well-posed or illposed, the uniqueness and existence of the solution, and the stability of the solution to errors in experimental data. Althoughmany of these issues can be addressedadequately and in a reasonably rigorous manner, the scope of the present paper will be restricted to an outline of the inverse problem and a method of solution, since our primary goal is to introduce a new method and to illustrate its performance. In particular, we present (i) a brief overview of the inverse problem relevant to interactions in colloids, (ii) the development of the details of an inversion method which does not require any preconceived functional form for the potential nor any a priori assumptions concerning the nature of the interaction forces, and (iii) a demonstration, using the method developed here, that shows not only that the details of the microscopic interactions are preserved in the structure but that these can, in fact, be recovered to a remarkable extent. The paper is organized as follows. We begin, in the next section, with a brief background on past attempts on the use of scattering data for inferring information about the details of interaction forcesin atomicliquidsand liquidlike collections of particles. We also draw attention in this section to the broader significance of the inverse problem as well as to the failure of the past attempts. Following this, the essential details of a perturbationtheory-based, predictor-corrector inversion method are presented. This method is then applied to three classes of "fluids" whose interaction potentials cover a broad range of repulsive as well as attractive forces. Finally, we discuss the implications of these results and other related issues that are relevant to the solution of the inverse problem for colloids. Analysis of Static Structural Data: The Inverse Problem In the past, X-ray, neutron, and light scattering techniques have been used primarily to follow (often qualitatively) the development of structure and ordering in the materials under study and to infer the global features of the interactions. In contrast, the inverse problem of extracting effective interaction potentials from observed structure has received very limited attention, for reasons
Langmuir, Vol. 8, No. 12, 1932 2899
that will become apparent from the discussion below. Nevertheless,the inverse problem is a very important one since the effective potentials contain a wealth of fundamental information on density dependence and possible long-range oscillatory behavior of the interaction forces (particularly in the case of Coulombic systems) and on cooperativeeffects. It is worth noting here that the above inverse problem has a much broader significancethan may be gathered from the specific focus of this paper. A detailed discussion of this falls outside the scope of the present paper. However, the broader implications of the inverse problem, overviews of past attempts to address the inversion of static structural data, and the pros and cons of the available options are already presented It is, therefore, sufficient for our purpose el~ewhere.'~J~ here to recognize the following two common approaches that have been popular in the literature; these are (i) inversions based on integral-equation theories, such as the Percus-Yevick (PY) and the hypernetted-chain(HNC) equations,16and (ii) fitting of static structure factors to those calculated from assumed potentials. Method i is a formal inversion technique, albeit highly approximate and inadequate, and dates back to the pioneering work of March and c o - w ~ r k e r s in ~ ~the J ~ context of interactions in rare-gasliquids and liquid metals. In contrast to method i, method ii is, in essence, a trial-and-error procedure. It is, nevertheless, the more popular of the two, both in the area of classical liquids as well as in the modern literature on colloid science. It is essential to note that both approaches have serious shortcomings. As has been already illustrated,16t20 the use of PY and HNC integral equations to extract effective potentials from known static structure factors is quite simple and straightforward, although the quality of the results is far from satisfactory. Numerous examples of the failure of classical integral-equation closures for inversion are given by Ailawadi.16 Conclusions similar to the above can also be drawn from the results of Nieuwenhuis and Vrij,2l who used the HNC closure to extract pair-potentials for poly(methy1 methacrylate) latex dispersions in benzene. Although not a formal inversion method, method ii is the more common one and has been used extensively in colloid science in the study of dense dispersions,ionic and nonionic micellar solutions, and microemulsions. For instance, this approach has been found very useful for determining parameters such as the aggregation numbers and the "dressed" charges of micelles and some biomolecules from small-angle neutron scattering data.22923 However, we emphasize that this method requires a priori information on the nature and the functional form of the interaction potentials and, therefore, could lead to serious errors since such information is not always available in advance. Any assumption concerning the details of the potential thus biases the results one obtains from the (16) Ailawadi, N. K. Phys. Rep. 1980,57, 241. (17) Sen, A. K.; Rajagopalan, R. Phys. Chem. Liq. 1991,22, 195. (18) Johnson, M. D.; March, N. H. Phys. Lett. 1963,3,313. Johneon, M. D.; Hutchinson, P.; March, N. H. Proc. R. SOC. London, A 1964,282, 283. (19) March, N. H. Liquid Metals: Concepts and Theory; Cambridge University Press: Cambridge, UK, 1990. (20) Rajagopalan,R.; Hirtzel, C. 5.J . Colloid lnterfoce Sci. 1987,118, 422. (21) Nieuwenhuis,E. A.; Vrij, A.J. ColloidlnterfaceSci. 1979,72,321. (22) Chen, S.-H.; Sheu, E. Y. In Micellar Solutions and Microemulsions: Structure, Dynamics, and Statistical Thermodynamics; Chen, S.-H.,Rajagopalan, R., Eds.; Springer-Verlag: New York, 1990, Chapter 1, pp 3-27. (23) Chen, S.-H.; Sheu, E. Y.; Klaus, J.; Hoffmann, H. J. Appl. Crystollogr. 1988,21, 751.
Rajagopalan
2900 Langmuir, Vol. 8, No.12, 1992 I
I
@. I
0
0
I
I 2
1
qu
Figure 1. hperimentalstructurefactorfor a charged dispersion “fitted” with structure factors calculated using two different potentials shown in Figure 2. The calcuhMS(q)for the Yukawa potentialis based on the rescaled mean-sphericalapproximation. The S(q)for the Sogamipotentialis based on Browniandynamics computer experimenta. The particle diametar is denoted by u. The data and the theoretical calculations are from Tata et al.% The f i e illustrates that unless proper attention is paid to the details of the experimental structure factor (includingthe low-q range), fitting methods could lead to ambiguous or incorrect resulte. See also Salgi et aL90for additional details.
of the failure is very instructive. However, these are not directly relevant to the subject of this paper, and therefore we close this section by noting merely that perhaps the single most important flaw in these methods is that, due the approximations implicit in them, they overlook or miss the basic physics underlying the relation between the interaction forces and the microstructure of the material. In the following, we outline a method that attempts to rectify the above problems.
A Perturbation-Theory-Based Inversion: Theory One of the keys to unscrambling the information in the static structure of the diaperisone is to devise a solution procedure that correctly sorts out the individual features of the potential which determine the come-grain and the fiie-grain structures of the fluid. For this, we appeal to the so-called van der Waals (vdW) picture of ‘liquids”, which posits that the core of the potential and the effective density of the fluid determine the coarser features of the structure while the rest can be attributed to the tail of the potential. This thesis forms the backbone of perturbation theories in statistical thermodynamicsand is satisfied by a large class of materials including the ones of interest here.% This vdW picture allows one to expand the configurational integral in a series in terms of a perturbation relative to a convenient reference potential. Since excellent descriptions of various forms of perturbation formalisms (for the forward problem), including the socalled optimized random-phase approximation (ORPA), which will be used here, are readily available in the literat~re,2~ the present section will be restricted to the central details of the inversion method used in this paper. Additional details and relevant equations (and extensions to cases not presented here) are given in GuBrin.% We begin with the experimental static structure factor, S(q),at known intervals of the scattering vector, q. The radial distribution function,g(r), is related to S(q)through g(r) = [p(2d3I-’
dq [S(q.)- 11 exp(iqr)
(1)
The magnitude q of the scattering vector, q, is defined as
= (47rn/X) sin ( W 2 ) where B is the angle at which the intensity of the scattering is measured, X is the wavelength of the radiation used in the experiment, and n is the refractive index of the carrier fluid in which the particles are suspended. (In the case of X-ray and neutron scattering, n is replaced by unity in the above equation.) The vdW picture implies that the (as yet unknown) potential, u(r), can be written as q
-1
I I I I I 10 20
0
rla
Figure 2. Yukawa and the Sogami potentials for the dispersion whose structure factor is shown in Figure 1. The experimental parameters are given in Tata et al.% method. An example of such a situation may be seen from the results presented in Tata et alea for a charged colloid. Tata et al. find that their experimentaldata can be ‘fitted”, with about the samelevel of accuracy, usingeither a purely repulsive, monotonically decreasing Yukawa potential or a potential displaying a minimum and a long-range attractive tail (see Figures 1 and 2)!25 A number of reasons can be advanced for the failure of the above methods and, indeed, an analysis of the sources (24) Tata, B. V. €2.; Sood, A. K.; Keaavamoorthy, R. pram^^ 1990,
34,23. (26) The attractive potential shown in Figure 2 (and used in Figure 1) is the -ami potential, which hati been a source of some controversy(see Sogami, I. Phys. Lett. A 1983,96,199. Sogami,I.; Ise, N. J. Chem. Phys. 1984,81,6320. Overbeek, J. Th.G. J. Chem. Phys. 1987,87, U )but , it ia not our intention here to be drawn into this debate. The purpcae of Figure 1 is to merely point out that “fitting” procedureti can lead to drastically different results and that they are not necamrily reliable since they focuo primarily on the global features of the structure factor. See the following paper for a further discussion of thir point: Salgi, P.; GuBrin, J.-F.; Rajagopalan, R. Colloid Polym. Sei. 1992,270, 786.
u(r) = uo(r)+ up(r) (2) where uo(r) is a suitable reference potential (usually an appropriate representation of the hard core) and up(r)is the perturbation (usually the soft, slowly varying tail). The potential uo(r) can be replaced, provisionally, by an effective hard-sphere potential U d ( r ) of diameter d, determined using an appropriate prescription (see eq 41, to define a trial potential, UT(?-),i.e. u&) = u&) + up@) (3) Numerous options exist for the selection of d (for a given (26) Castillo, C. A.; Rajagopalan, R.; Hirtzel, C. S. Aeu. Chem. Eng. 1984,2, 237. (27) Hansen, J.-P.; McDonald, I. R. Theory of Simple Liquids, 2nd ed.;Academic Press: New York, 1986. (28) GuBrin, J.-F. PhD Dissertation, University of Houston, 1992.
Langmuir, Vol. 8, No.12, 1992 2901
Interaction Forces in Colloidal Fluids
u ~ ( r ) )we ; recommend the following criterion, which is known to lead to better thermodynamic consistencym $[ayd(r)/adl(exp[-Buo(r)1- exp[-Bud(r)l) d r = 0 (4) where yd(r)is the so-called cavity function27corresponding to U d ( r ) . The structure factor corresponding to UT@) can be derived from ORPA and is given by
+ Sd(q)p&p(q)l
(5) where p is the number density, B1is the product of the Boltzmann constant, kg, and the temperature, T,Sdq) is the structure factor of a hard-sphere fluid of diameter d, and up(q)is the Fourier transform of up(r). Associated with this up@)is what is often referred to as a “renormalized potential” +(r)that makes the radial distribution function inside the core become zero, as demandedby the excludedvolume criterion. The Fourier transform of the function +(r) is related to Sd(Q) and up(@by ST(q) = Sd(q)/[l
-[Sd(q)12flup(q)/[l + &k‘)PBUP(q)l (6) The renormalized potential +(r)can be used to derive a higher level of approximation for g(r) than implicit in eq 5 +(q)
g(r) = Yd(r) exp[+(r) - Buo(r)l (7) which takes into account the softness of the core of the potential. The inversion of a given structure factor proceeds as follows. One first needs an estimate of the effective hardsphere diameter, d, that represents the first approximation to the excluded-volume effect of the core of the potential on the structure. This can be obtained from the wellknown relation27between the isothermal compressibility and the known structure factor at q = 0, i.e. [l/S(0)1 = [a(BP)/dplT (8) where the pressure P (the osmotic pressure in the case of colloids) can be substitutedfrom an equation-of-statesuch as the Carnahan-Starling equation for hard-sphere fluids. Once d is estimated, &(q) can be obtained using the solution of the Percus-Yevick equation.27 Then, eq 5 and the given structure factor data [substituted in place of ST@)] lead to a f i t estimate of the perturbative part of the effective potential and subsequently to the renormalized potential through eq 6. Since +(r)and Yd(r) are now known,the “softened” core w(r)can then be determined from eq 7 with the experimental radial distribution function,g,,,(r), in place of g(r). The first approximation to the effective potential is thus given by BuW = In 1 3 ~ ~ ~ ) / ~ e+x+(r) p ( ~+)Bup(r) l (9) The f i t approximation may need further refinements since the effective hard-sphere diameter, d, is estimated originally rather crudely using eq 8. One can refine the result from eq 9 through successiveiterations by returning to the above scheme with a new hard-sphere diameter obtained (using eq 4) for the extracted potential. We call eq 9 the ‘predictor”, and the refinement of the hard-sphere diameter the “corrector”. The above predictor-corrector method converges very rapidly (often within a few iterations) and, in addition, incorporates the essential physics more rigorously than the integral-equation methods. Moreover, the given structure factor is matched exactly by that of the inverted potential (in contrast to the fitting method (ii) discussed
earlier; see, for example, Figure 1);this implies that the extracted potential is a truly structurally consistent effective potential which contains all the structural information contained in the measurement.
Results We now illustrate the above method for a few classes of fluids representative of a broad spectrum of interaction potentials found in colloidal systems and simple liquids. A general discussion of the method and of some related issues follows subsequently in the next section. As mentioned earlier, so far radiation scattering techniques have been used in colloid science (asin other areas of condensed matter physics) primarily to investigate the general features of the microstructure of the dispersions and to study the accompanying structural of phase transitions. As a consequence,high-quality data, covering a broad range of scattering vectors (this may require combining light scattering and small-angle X-ray or neutron scattering), are not available for well-defined model colloids. In view of this and in order to test the inversion method for well-defined potentials, we take the usual approach followed in statistical mechanics. That is, we consider some known potentials and use a reliable forward theoretical procedure or computersimulations to generate static structure factors, S(q),and radial diatribution functions,g(r). These are then u d in the inversion scheme to investigate how faithfully the features of the known potentials are recovered by the inversion method. Although a number of options exist for generatingS(q)’s from given potentials, it is essential that the ones chosen provide sufficiently accurate results for the structure factors since, as discussed in the subsequent sections of the paper, the accuracy of the inverted potential is sensitive to the accuracy of S(q) in the low-q region. Since sufficiently accurate theoretical techniques are not available for high densities currently, the need for careful computer simulationsis inevitable. We have chosen Monte Carlo (MC) simulations (with a large number of particles, ranging from about 3000 to 11000) for generatingg(t) and S(q). The need for highly accurate computer simulations [with appropriate estimates for statistical errors in the simulations (e.g., confidence intervals)] even for testing the accuracy of the various closures for the integral equations is illustrated in Salgi et al.,3O and some of the observations and details relevant in the present context may be found there. In addition to MC simulations, we have also wed a number of standard integral-equation theories for generating S(q)’s. Most of the integral-equation theories are seldom sufficiently accurate for testing inversion procedures of any kind. Our purpose in considering them, however, is to use them merely as a convenient means to generate insufficiently accurate S(q)’s so that the sensitivity of the inversion procedure to inaccuracies in S(q)’s may be examined. For the purpose of this paper we shall restrict our attention to the so-called reference hypernetted-chain (RHNC) method here. The RHNC method is a solution procedure devised for the Ornstein-Zernike integral equation for the total correlation function, h(r) [=g(r)- 13, with a closurerelation that attempts to account for the higher-order contributions not included in simpler closures such as the hypernetted-chain (HNC) closure. Details of this approach may be found in standard monographs on liquid-state physics.27 There are two reasons for considering the RHNC theory in the present (30)Sa&, P.;GuBrin, J.-F.; Rajagopalan, R. Colloid Polym. Sci. 1992,
~
(29)Lado, F. Mol. Phys. 1984,52,871.
270, 785.
Rajagopalan
2902 Langmuir, Vol. 8, No. 12, 1992
ITwo Model Potentials 1 12-6 LJ Potential
1.5
1
Exp-6 Potential
--
1.0
w
E
'1.
1
\
1 . L
m
0.5
I
*
i i
Id*
Monte Carlo
-
0.0
2
3
4
5
6
7
8
Figure 3. Model Lennard-Jonesand Exponential-6potentials used for calculating the structure factors shown in Figures 4-6. See the text for details.
context: (i) This theory is one of the best forward methods available currently,and we have found it to be very accurate over a wide range of densities for two of the potentials considered in the present paper. Thus it provides a convenient way for generating S(q)'swith arbitrarily small intervals in q-space. (ii) The RHNC theory has been used along with molecular dynamics in a predictor-corrector inversion scheme by Reatto et recently, and the inclusion of the RHNC method for generatingS(q)'s allows us to draw some parallels and emphasize some differences between our inversion method and that of Reatto et al. Three types of potentials have been selected here for testing the inversion procedure: 1. The first is a Lennard-Jones (LJ)potentid, which is chosen in such a way that, although it is not a colloid potential, it mimics the essentialfeatures of a typical colloid potential in an electrostatically stabilized dispersion with a moderate "secondary" minimum. This implies that the electrostatic "barrier" is sufficiently large that the colloidal particles do not sample interparticle separations where the potential may display a primary minimum. As wellknown, the LJ potential has the form BuLJ(r)= 4~e[(u/r)~~ - (u/~)'I (10) where e is the minimum in the potential and u is the collision diameter (defined as the value of r at which the potential is zero (see Figure 3)). 2. The next is an exponential-6 potential, which has a much softer core than the LJ potential. This potential has the form
BuExpG (r) = A[exp(-Br/u) - (u/r)'I (11) where A and B are two parameters that are chosen according to a prescription described below. 3. Finally, we choose a colloidal potential of the Dejaguin-Landau-Verwey-Overbeek (DLVO)form3that includesmore complexfeaturestypical in colloidalsystems. The functional form and the parameters are described below. There is a number of reasons for the choice of these potentials. In addition to the fact that they cover the types of features expected in colloid potentials, they also have additional attributesthat allow us to test the inversion method rigorously, at the limits of the assumptions that have gone into the invereion formalism. For instance, the LJ potential has a steep core and a moderate range of (31) Reatto, L.;Levesque, D.; Weis, J. J. Phys. Rev. A 1986,33, 3451.
0
1
2
3
4
5
6
q (l/A, Figure 4. Structure factors calculated using Monte Carlo simulationsand the referencehypernetted-chain(RHNC)theory for the U potential shown in Figure 3. attraction that is fairly slowly varying. The assumptions implicit in the inversion method (e.g., decoupling the potential in terms of a harsh reference potential and a slowly varying perturbation) are expected to be satisfied relatively easily by this class of pair-potentials. On the other hand, the exp-6 potential used here has a relatively softer core (while the attractive tail remains the same), namely, a condition under which the replacement of the core with an effective hard-sphere potential is stretched to a limit. Similarly, the colloid potential chosen also contains features that allow one to test the inversion in other ways. As illustrated later, this potential has a steep core (as would be expected in most dispersions), but the "tail" of the potential has sharply varying attraction and repulsion, in addition to being rather long-ranged. More complex cases are considered by Gu6rin.28 We begin with the resulta for the LJ and exp-6 potentials. The two pair-potentials are shown in Figure 3. We have arbitrarily chosen u to be 3.4 A and Be to be 0.51 in eq 10. The parameters for the exp-6 potentials have been chosen in such a way that the minimum in the potentials and its location remain the same for both potentials (see Figure 3). As also evident from Figure 3, the tails of both potentials are practically identical, and only the cores differ. This allows us to test the relative accuracy of the inversion method with respect to two qualitatively and quantitatively different cores while maintaining the rest of the potentials the same. The structure factors for these two potentials are calculated for a range of densities using MC simulations and the RHNC theory as input to the inversion. The MC results for a reduced density of pa3 = 0.8 (i.e., p = 0.0204 A-1) are shown in Figures 4 and 5 for the U andexp-6potentials, respectively. These are based on simulations in which 10 976 particles are used, with a maximum cutoff distance of half the length of the simulation The averages are taken over at least 20 X 1 0' configurations, after confirming that the effects of the initial configuration are eliminated. Also shown in Figures 4 and 5 are the results based on the RHNC theory. It is evident from the figures that the RHNC theory leads to very accurate values of S(q) for the cases shown. The effect of the softness of the core on the long-range configurational structure (i.e., the low-q range of the structure fador) is strikingly brought out in Figure 6, which (32) The radial distributionfunction can be extended beyond the c u t off distance using the Verlet algorithm (Verlet,L. Phya. Rev. 1968,166, 201) or other methods, if necessary, although extreme care muet be exercised to make sure that the g ( r ) and the corresponding S(q) are mutually consistent.
Interaction Forces in Colloidal Fluids 1.4
Langmuir, Vol. 8, No. 12, 1992 2903 3
m
* oe,
1.2
Exponential-6 Potential
:
LJ (12-6) Potential 2
1.0 n
LI
A
3 cn
W
0.8
7 a
1
0.6
0 0.4
0.2
-1 0
1
2
3
4
5
6
q (l/A) Figure 5. Structure factors calculated using Monte Carlo simulationsand the referencehypernetted-chain(RHNC)theory for the Exponential-6potential shown in Figure 3. The triangle represents a point obtained by the Verlet
2
3
5
4
r
6
7
8
(A>
Figure 7. Comparison of the known W potential with the potential invertedfromthe correspondingstructurefactorshown in Figure 4. Exp-6 Potential
2 1.0 n
w
W
m 0.5
0
1
2
3
q
4
5
2
6
(l/A)
Figure 6. Structure factors calculated using the reference hypernetted-chain(RHNC)theory for the two potentials shown in Figure 3. shows the structure factors from Figures 4 and 5 on a single plot. One sees a sharp increase in S(q) for low-q values for the exp-6 potential indicating the compressibility of the fluid interacting through the softer potential (see, for example, eq 8 for the relation between S(0) and the isothermal compressibility). In contrast, the dominance of the core of the potential in determining the long-range statistical ordering in the LJ fluid may be seen by the larger value of the first peak in the corresponding S(q) and the low value of the isothermal compressibility. Thus the structure of the LJ liquid under the given conditions of temperature and density is consistent with the assumptions implicit in the inversion scheme, and an excellent agreement is seen between the known potential and the inverted potential in this case (see Figure 7). What is noteworthy is that almost a similar level of success is seen in the case of the exp-6potential as well, as illustrated in Figure 8. The inversion method recovers both the core and the tail of the potential extremely well. It is also worth noting that the calculations converge within about 5-10 iterations in each case, the convergence is uniform, and the intermediate values of the potential at any r "bracket" the final value. Next we examine a more complex interaction potential, namely, a DLVO potential, which is described below in some detail. The model colloid potential used here has a hard core (of unit diameter), a strong London-van der Waals attraction, UA,a weak electrostatic repulsion, U E , and a steep Born-type repulsion near the core, U B
3
4
5
r
6
7
8
Figure 8. Comparison of the known Exponential-6 potential with the potential inverted from the corresponding structure factor shown in Figure 5.
&(r) = bu,(r) + &E(r) + @B(r) (12) where r is the center-to-center distance of separation between the particles in units of particle diameter and @ is again ( k ~ T ) - l The . expression for the London-van der Waals interaction u~between two identical spheres is given by3i26
= - ~ & / ( 2 r ~ )+ 1/[2(r2 - 1)1+ In [(r2- 1)/r2]j (13) where
Nh = @H/6 The Hamaker constant, H,in the above equation is a material property. The repulsive contribution, UE,in eq 12 arises from the electrostatic interaction between the diffuse portions of the electrical doublelayers surrounding each particle. The general expressions for electrostatic interaction energies are rather complicated and depend on whether the interacting surfaces keep their potentials constant or their charge densities constant. We have wed an expression that reproduces the essential features well for particles with constant surface-charge densities and thick electricaldouble layers,lJ Le., small inverse-screening length, KU < 2.5 Bu&) = NE[exp(-Ka(r - l))l/r
(14)
where a is the radius of the particles (equal to 0.5 in our dimensionless unite) and the premultiplying factor NEis
Rajagopalan
2904 Langmuir, Vol. 8, No. 12, 1992
5t
0
4
Inverted from RHNC S(q) at Volume Fraction 0.15
11 0
1.00
1.05
1.10
r/d
1.15
1.20
1.25
Figure 9. Comparison of the known DLVO potential with the potential invertedfrom the corresponding RHNC structurefactor (for a volume fraction of 0.15).
I I
1.00
I
1.05
Known u(r) Inverted from RHNC S(q) at Volume Fraction 0.4
r/ d
1.10
1.15
Figure 10. Comparison of the known DLVO potential with the potential invertedfrom the corresponding RHNC structurefactor (for a volume fraction of 0.4).
related to the contact potential. The other repulsive contribution, UB,despite its very short range, needs to be taken into consideration for dispersions with moderate or low electrostatic interactions in order to adequately represent the steep repulsion at near-contact distances. An expression of the following form has been used here for the Born repulsion U B so that the appropriate shape of the pair-potential could be generated &(r) = NB(r (15) with NB being a scale factor that we call the Born parameter. (Instead, a hard-sphere interaction could have been used. Other options are also possible, but the one chosen serves the purpose well without any loss of generality.) The numerical values that have been used for the dimensionless parameters are NLO= 1.025, NE= 10.0, KO = 3.75, and NB = 2.275 X The resulting pair-potential haa both a primary minimum and a secondary minimum, with a moderate energy barrier in between. The attractive force immediately beyond the primary minimum is significant in magnitude. The potential resulting from eqs 12-15 is shown in Figures 9 and 10, along with the potentials inverted from the structure factors for two different volume fractions. In this case the structure factors and the radial distribution functions (not shown) used for the results shown in Figures 9 and 10 have been generated using the RHNC theory. We have found that the RHNC results compare very well with the results of Monte Carlo simulations for low and moderate densities (e.g., volume fractions up to
about 0.2). At larger volume fractions (about 0.3or above), even though the radial distribution functionsappear quite accurate, the structure factors in the low-q region are in error by about 20% or more. Since the S(q)values at low q's w i l l have significant contributions to the extracted potential near the core of the potential, we have considered two volume fractions, namely, 0.15 and 0.4,to contrast the results with each other. One may note from Figure 10 that the inversion method predicts the near-field interactions (i.e., close to the primary minimum) quite accurately even with the structurefactor results for the highest volume fraction (0.4)tried here. It can also be seen from Figures 9 and lothat the inverted pair-potentials are in excellent agreementwith the known potentid. All the major features of the potential and the shape of the potential are reproduced very well. As noted earlier, the tail of the potential (i.e., the potential beyond the primary minimum) has a steeply varying attraction followed by a slowly varying repulsion as r increases. (The tail shows a very small, roughly 0.05k~T,secondary minimum in this case.) Nevertheless, the inversion methods perform quite well. Overall, the results demonstrate that the inversion method presented here is capable of reproducing both the near-field and the farfield behavior of the potential very well, and in a physically consistent manner, as long as sufficiently accurate structural data are available. The comments and observations concerning the stability, robustness, and convergence of the method made earlier apply in this case as well. We now proceed to a general discussion of the method and its application to experimental data and a brief discussion of why the present method succeeds where others fail.
General Discussion One essential feature of almost all (formal) inversion methods that have been used in the past is that they depend, directly or indirectly, on successful retrieval of the direct correlation function, c(r), from the known structure factor. This function is related to S(q)by the relation c(q) = [S(q)- lI/S(q) (16) where c(q) is the Fourier transform of c(r). However, c(r) is sensitive to errors in S(q),and this may be compounded further because of numerical errors introduced by Fourier transforms (when the q-space information is transformed to real space). This is a major reason for the failure of inversion approaches based on simple integral-equation closures such as the PY and HNC approximations, which we discussed briefly in an earlier section.= This statement also applies to a number of other approaches such as the ones based on the mean-spherical approximation (MSA) and generalized mean-spherical approximation(GMSA), which are not discussed here but are outlined elsewhere." The same is true in the case of the single-step RHNC method (i.e., without a ycorrector" step as done in our method), although the latter is generally superior to the others mentioned above because of the appearance of the reference hard-sphere bridge-function correction (to the HNC closure) in the predictor. The single-step RHNC inversion, however, fails at high densities as shown by (33) Of course, in the caae of the "fitting"methods diecueeed in the context of Figure 1, the 'extracted" potentiale and their parametere are alreadyhighly b i d by the a priori ~ ~ m pthatt goiinto ~ the~assumed potentiale. Even if there are strong reasons for selecting a particular form of the potential in advance,the fitting procedures are handicapped by the well-knowninsensitivityof the forward theories (Le.,the structure factomcalculated using the theories) to someof the detailsof the potentials.
Langmuir, Vol. 8, No. 12, 1992 2906
Interaction Forces in Colloidal Fluids
Reatto et al.,3l whereas our inversion leads to correct potentials. Even though c(r) does not enter our method directly, it is implicitly present in the calculation of up(r), but ita effect is diminished to some extent since up@)is proportional to the difference between c(r) of the entire potential and C&), the direct correlation function of the hard-sphere reference fluid. A detailed discussion of the nuances involved in this and related aspectsof the inversion methods requires a considerable amount of additional details about the physical basis and the mathematical structure of various inversion schemes and will not be attempted here as it would require a digression from the primary focus of the present paper. We, however, note that the same issue concerning inaccuracies in c(r)is what makes single-step inversions inadequate. And, for the same reasons, the corrector also needs to be designed very carefully. In particular, what information is “fed back” to the corrector at each subsequent iteration and what is left out are crucial. This is the reason we chose to correct only the influence of the effective excluded volume at eachstep;this is, of course, also consistent with the general conceptual underpinning of the present perturbationtheory-based approach. We shall discuss these and other related issues in detail in a future publication. Here, we return to some of the other issues important in the context of colloidal systems. Three issuesthat are important in applying the inversion methods to practical systems are (1)the sensitivity of the method to “errors” in S(q) and the related issue of the extension of experimentally measured S(q)’s beyond the finite range over which they are usually accessible, (2) the effects of polydispersity of the dispersions, and (3) the effects of multibody interactions. The extension of experimental S(q)’s is very important from a practical perspective for two reasons: (i) the zero-q limit (i.e., q values smaller than qpea~,the location of the first peak in S(q))is not always accessibleexperimentally to a sufficient extent, and one may need to extend the data in this region by suitable approximations; (ii) a similar situation may also arise for large q values (Le., one may find that sufficiently large values of q for which S(q)tends to unity are not reachable experimentally). The issue of experimental errors in S(q) has been addressed by GuBrin%to a limited extent and our calculations show that a 10% error in the most important region of S(q),i.e., the low-q region (especially for q I 0.8qpea~)causes roughly a 5% error in the region of the minimum in the potential. This is in sharp contrast to previous methods (including the PY and HNC methods), where one sees a doubling of the error in the region of the minimum in the potential.l6sM A similaranalysis for errors in the region beyond the second peak in S(q) reveals that the inversion is insensitive to errors for large q’s. Thus the proposed method is very stable relative to numerical and experimentalerrors. These and related issues will be discussed elsewhere, but a word of caution is warranted here. The fact that the inversion is stable with respect to errors in the large-q region does not imply that the truncation of S(q)in the large-q region is of no consequence. The truncation of S(q)before S(q) reaches the asymptotic limit of unity sufficiently closely (or the nonavailability of experimental data in this region) has a severe effect on the quality of the radial distribution function derived from the data. This is an important problem that remains essentially unaddressedadequately at present, although maximum-entropy formalisms and other techniques have been explored to some extent28in this context. (34) Ballentine, L. E.; J. C. Jones, Can. J. Phys. 1973,63, 1831.
It is also well-known that the effects of polydispersity are equally important. The structure of polydisperse colloids has not been studied in great detail so far, but a reasonable amount of analysis based on PY equations (for hard spheres and adhesive hard-spheres) and other integral-equation formalisms (for charged spheres with and without attraction) is currently available (see a forthcoming review on this topic3s). These resulta show that the average structure factor deviates only slightly from the structure factor based on the average particles as long as the standard deviation of size-polydispersity (or polydispersity in other parameters) is within about 10% of the average size (or the average value of the appropriate parameter). In such a case, the inversion will lead to an average effective potential. If decoupling approximation^^^ can be used, then one can again obtain an average effective potential. On the other hand, for binary or multicomponentsystems, the inversion method presented here can be extended appropriately (as shown by GuBrin%),if the partial structure factors are known. The effects of multibodyinteractions, if significant, would also lead to an effective potential when the present inversion method is used. In the present paper, we have relied on computer-generatedor theoretically determined S(q)’s based exclusively on pair-additive interactions in order to avoid the complications due to multibody interactions and to focus specifically on the conceptual foundations of the inverse problem and on the detaila of the proposed inversion method. The examination of multibody interactions requires more elaborate analyses of the structure and ita relation to the potentials, but is feasible to some extent. Some additional discusaion on this is presented elsewhere17using the inversion of neutronscattering data36 on argon near ita triple point. The inversion schemepresented here may also be supplemented further by imposing other “consistency” criteria (such as forcing osmotic pressure or some other chosen property based on the inverted potential to be the same as that of the actual system) so as to guarantee a useful effective potential in more complex cases. Finally, it is also important to note that the present work, and inversionmethods,in general, are of significance to the forward problem as well. In fact, as has been pointed out by March,lginversion methods, if accurate, provide a much better and a much more rigorous test of the accuracy of structure factors calculated by forward theories. For example, it may be possible, using method such as the one presented in this paper, to check the accuracy of the theoretical techniques available for predicting the structure of dense liquids and colloids. We have shown,28in one suchexample involvinga long-range,oscillatory,liquidmetal potential, that the RHNC theory and another integral-equation formalism due to Zerah and Hansena7 (known as the HMSA method) lead to different structure factors for the same input potential under identical conditions of temperature and density. Our computer simulations show that, under the conditions used in the example, the RHNC predictions for the structure factor are accurate, whereas the HMSA theory underestimates the low-q structure factor by about 1695, and thia is confirmed by the results of the inversion. Specifically, the present inversion method, when applied to the RHNC (orthe MC) structurefador, recoversthe original potential; on the other hand, when the HMSA S(q) is used in the ~
(35)Salgi,P.; Rajagopalen,R.Submittedfor publication in Adu. Colloid Interface Scr. (36)Yamell, J. L.;Katz, M. J.; Wenzel, R. G.; Koenig, S. H. Phys. Rev. A 1973, 7 , 2130. (37)Zerah,G.; Hansen, J.-P. J. Chem. Phya. 1986,84, 2336.
2906 Langmuir, Vol. 8, No. 12, 1992
inversion, a slightly different potential (with about a 10% differencenear the core from the actual potential) results. Since subsequent MC simulations using the potential inuerted from the HMSA S(q)data reproduce the HMSA S(q),it is evident (i) that the present inversion correctly identifies the potential appropriate for the HMSA S(q) and (ii) that the HMSA closure is not accurate at the temperature and density used. Notice that these observations also imply that the performance of the inversion procedure presented here is not tied to, or dependent on, any specific closure relation or any theoretical technique used for producing the test data. The last statement is an important one since we have used the RHNC method to examine the results of our inversion in some of the examples discussed in this paper. The spirit of the RHNC theory is akin to that of the perturbation theory used in the present paper for inversion, since the former blends, in a certain sense, the integral-equation formalism with the perturbation approach.% However, the inversion formalismused here and the RHNC theory are sufficiently different, and, as evident from the above discussion and the results presented earlier, the use of the RHNC theory for generating S(q)does not favorablyinfluencethe success of the present inversion. This is also evident from the results presented for the DLVO fluid (especially at the higher density), for which the RHNC structure factor is inaccurate in the low-q region; the inversion in this case leads to an appropriately different potential, which reproduces the input S(q) when subjected to a rigorous test using computer simulations. One may also emphasize at this point the importance of the need for carefully executed computer simulations for testing any inversion. We chose to use the RHNC theory (and others not reported here) to generate variations in S(q)’s so that the sensitivity of the inversion to inaccuracies in S(q) could be examined relatively easily, but the RHNC theory is not used as a replacement for computer simulations. As noted by Reatto et al.,31even the RHNC method,%the best liquid-state theory currently available, does not have the required accuracy at very high densities (close to the triple point) in the forward route for providing data of sufficient quality for inversion (see, also, the excellent review by Reattom for a discussion of related issues). This has, in fact, prompted Reatto et al. to propose the use of computersimulationsas the corrector at each step. In combinationwith the RHNC method (in the reverse route) as the predictor, their predictorcorrector scheme offers a powerful inversion method, although the use of simulations as the corrector carries with it the attendand need for smoothing and extending the simulated data at each step in the inversion. Our objective in the present paper, however, has been to develop an inversion scheme that does not require computer simulations within the inversion scheme. The results we have presented are encouraging. We shall present a comparison the present method with that of Reatto et al. in a future publication. Closing Remarks When confronted with the task of extractinginformation about interaction forcesfrom experimental structure factor data (assuming that they are sufficiently accurate, unambiguous, and available over a broad enough range of q ) one has, aa mentioned in the beginning of the paper, at least two options from which to choose, namely, (i) a (38) Rosenfeld, Y.; Ashcroft, N. W. Phys. Reu. A 1979,20, 1208. (39) The RHNC method is also known 88 the modified hypernettedchain (MHNC) method; the latter term is used by Reatto et aLS1 (40) Reatto, L. Philos. Mag. A 1988,58,37.
Rajagopalan
“formal” inversion method of the type suggested here or (ii) fitting the observed structure factor using assumed forms of interaction potentials based on physical grounds. As noted the second option may lead to misleading results since any bias implicit in the selection of the trial potential may in turn bias the outcome. On the other hand, when there are sufficientlystrong physical grounds for choosing the form of the potential, the fitting method is relatively straightforward and haa been found to be quite useful in analyzing data on association ~ o l l o i d s . Such ~~*~~ an approach has also been used relatively successfully by Dharma-WardanaandAersq1for atomicliquids. However, the above method is prone to seriouserrorswhen the nature or the details of the interaction potentials are not clear in advance; a case in point is steric or polymer-induced interactions in colloids. (In the case of colloids as well as in simple liquids, the fitting method is often found appealing since it appears to circumvent the need for extending the experimental data if the latter are available only over a limited range of q. Nonetheless, it is important to note that, regardless of the method used, the problem remains ill-posed if the data are not available (or cannot be extended) over a sufficiently large range of q.) The first option, namely, resorting to rigorously-based formal inversion schemes (such as the schemesuggested by Reatto et al.) avoids some of the difficulties mentioned above, but requires a substantial amount of effort. One of the reasons why we have explored the perturbation-theorybased method we describe in the present paper is that it provides very good results, at least for the range of conditions explored here, with very little investment in effort. Further, it is relatively easy to implement and requires very little computationalresources (nothingmore than a desk-top computer or an equivalent is needed). In closing, we note that past attempts to gather microscopic interaction potentials from the macroscopic structure by fitting assumed forms of potentials have led to the general belief that such microscopic detaile are not preserved in the structure. This belief has been largely due to the often poor resolution of the forward formalisms. We have demonstrated,through the results presented here, that such details not only are preserved but can, in fact, be recovered to a large extent if the inversion method places proper emphasis on the dominatingfeatures of the interaction potentials. Acknowledgment. I thank the Editors of Langmuir for inviting me to submit a manuscript for the special issue of the journal in honor of Jean Perrin, whose work has influenced me greatly in directing my o m interests toward molecular physics and supramolecular fluids. I also thank the National Science Foundationfor partial support of the work reported here (through Grant No. CTS91000139). Acknowledgment is also made to the donors of The Petroleum Research Fund, administered by the American Chemical Society. The computing time for the simulations, provided by the Pittsburgh Supercomputing Center through Grant No. PHY900036P, is also appreciated. Dr. L. Reatto (Dipartimento di Fisica dell’Universita, Milano, Italy) and Dr. D. Levesque (Laborahire de Physique Thhoriqueet Hautes Energies, Universite de Paris XI, Orsay, France) were kind enough to share with us the original data used in their paper,3l for which I am grateful. I have benefited greatly from discussions with Dr. A. K. Sen (SahaInstitutefor Nuclear Physics, Culcutta, India), a former associate of mine, and with Drs.J.-F. Guhrin and P. Salgi, two of my former students. (41) Dharma-Wardana,M. W. C.; Aers, G. C. Phys. Rev. E 1988,28, 1701.
