J. CHEM. SOC. FARADAY TRANS., 1994, 90(9), 1245-1250 Integral Equation Theory for Associating Liquids :Highly Asymmetric Electrolytes Jun Wang and A. D. J. Haymet" School of Chemistry, University of Sydney, NSW 2006,Australia The associating liquids formalism of Wertheim is used to study association in highly asymmetric model electro- lytes, with charge asymmetry + 1 : -10, +I : -20 and size asymmetry 1 : 13 and 1 :8. These models are designed to capture some features of ionic colloids, micelles and proteins, but fall short of a complete descrip- tion. Explicit predictions are made for the fraction of associated species in the electrolyte. The theory yields improved agreement with structural data from existing computer simulations. Moreover, it predicts correctly that a large fraction of the large, highly charged ions are associated. However, the current level of theory is not a complete description of association.1. Asymmetric Electrolytes This paper studies simple models of highly asymmetric elec- trolytes using the associating liquids formalism of Wertheim' adapted by us to treat molecules interacting uia continuous, spherically symmetric potential^.^.^ Previously we have studied exactly the same models using the hypernetted chain (HNC) integral equation4 and bridge function corrections to HNC (HNC + B).' The motivation for this study is stated simply. Solutions of polyelectrolytes such as proteins, micelles and ionic colloids may behave quite differently from simple electrolytes, due in part to the asymmetry in charge between the polyelectrolyte molecules and the surrounding counter-ions and co-ions in the solution.Our view4 is that the dynamical properties of the system, including reaction rates, can be understood only if the spatial correlations between all ions in the solution are known. There are many alternative approaches to this problem6 and we do not intend to review these here. Instead, we address just one feature of these systems, which has attracted little interest to date, namely the association of ions in solu- tion. A simple model is all that is needed to see this feature. We consider a highly charged spherical electrolyte, denoted for example XZ0-, surrounded by monovalent counterions A+. (In this calculation the concentration of co-ions B-is zero.) Our goal is to predict the equilibrium concentration of all possible associated species, such as the 'dimer' X2'-A+, the 'trimer' X20-(A+)2, and higher-order 'n-mers' X2'-(A'),-'.Most researchers seem to feel intuitively that in solution the polyions exist not as 'bare' ions but rather as highly associated groups of ions which are more nearly elec- trically neutral. We agree with this proposal and this paper reports our partial progress in addressing this goal. A weak point in our treatment is the approximation of the aqueous solvent by a structureless dielectric continuum. Else- where we report our treatment of atom-based models of water by molecular dynamics (MD) simulation' and integral equation theories.* Nevertheless, since our goal is to study ion association, and the present simplified model of the poly- electrolyte is capable of studying this effect, we present our results and leave for the future the unification of this approach with a molecular-level description of the solvent.Our results are presented in Section 4 and compared with those from Monte Carlo (MC) simulations for + 1 : -10 electrolyte by Vlachy et aL4 and with those from a MD-reference HNC (RHNC) procedure for + 1 : -20 electrolyte by Lin~e,~q" in which MD simulations for the reference system and the RHNC perturbation method have been com- bined. 2. Model Asymmetric Electrolyte The interaction between ions is given by4." where the subscripts a and b denote + for counter-ions with z, = 1 and -for polyions with valences z-= -10 and -20, gab = O, + Ob, b-' = kT, k is the Boltzmann constant and T the absolute temperature. The parameters A,, and go used in this work are those used by Vlachy et aL4 and Linse,' and are listed in Table 1.The solvent is approximated by a dielectric continuum with a Bjerrum length Be'/& x 7.15 A, to model4 aqueous solutions at T = 298 K.This model mimics salt-free aqueous solutions : the free amphiphiles present in real micellar solutions are neglected. The potential energies are plotted in Fig. 1. We note that there are two length scales in each of the potentials. The asymmetry in size is ca. 1 : 13 for the + 1 : -10 electrolyte and 1 : 8 for the + 1 : -20 electrolyte.3. Brief Summary of Associating Electrolytes Associating Liquids Theory Ion clustering in highly asymmetric model electrolytes is treated by the following the~ry.~.~ Ion association follows from the separation of the pair potential, where U$)(r) and U$!(r) are the non-associating and associ- ating parts of the potential,' and Go, is the Kronecker delta Table 1 Potential parameters for + 1 : -10 and + 1 : -20 electro-lytes 1 : -10 0.08299 0.7 9.15 18.04 18.04 18.04 1 : -20 0.09960 2.0 15.0 0.795 15.89 371.876 1246 100 90 80 70 60 I-% 50'40 30 20 10 0 &=--_________-_-------___---I I I 1 I-10 0 20 40 60 80 100 r/A Fig. 1 Potential energies of the model electrolytes studied here: +1 :-10 (-) and +1 : -20 (---) electrolytes function.Unlike the strongly associating liquids, where dimers and trimers are formed by 'chemical bond^',^*^^*^^ eqn. (2) introduces a mechanism for the association of unlike ions. Using the 'energy definition',2 the associating part of the potential energy between unlike ions U?)-(r) is written as where U, is determined by minimizing the Helmholtz free energye2v3 OrnsteieZernike Equation The total correlation functions, h$(r) =g$(r) -hOiSoj, are determined from the analogue of the Ornstein-Zernike (OZ) equation, 'ab = cab +1cad * Pd * 'db (4) where * denotes convolution. The matrices are defined to be (5) where the subscript 0 indicates a non-associated species, and the subscript 1 an associated species. The number density for species a is pa = zp:-1, where pi-(for n = 1, 2, 3, ...) denotes monomer, dimer, trimer and higher 'n-mer 'densities.Eqn. (4) is renormalized as usual to treat the long-range Coulomb intera~tion.~.~ This renormalized OZ equation in Fourier space may be written W)= {I-cz + &(k)Plm)PI -{CZ+&(k)Pl6Qk) x 1+ +P6Wl -Se(W (6) where Z is the 4 x 4 unit matrix. The matrix M(k) = 6@k) -Sqk) is obtained from the Fourier transforms of 6H(r) and SC(r).The explicit forms of matrices p, SH(r) and 6C(r) are given in our study of the +1 : -3 ele~trolyte.~ J. CHEM. SOC. FARADAY TRANS., 1994, VOL. 90 The short-ranged functions 6h$(r) and 6c$'(r), the elements of the matrices 6H(r)and 6C(r),respectively, are defined to be 6c$(r) = c$(r) +dOi6, flu:& (7) and ah;;@) = h$(r) -q$(r) (8) where U:&) is the Coulomb part of the total potential, uab(r), and q;$r) = 6oi 6oj q"DbH(r) where qgH(r)is the Debye-Huckel (DH) potential of mean f~rce.~,~ Approximate Closure We now introduce the approximation into the theory. First, we assume that the ions can form only dimers and trimers and examine the consequence of this approximation.Note that there are two distinct possible trimers in the system. We use the symbol +-+ to denote the trimer consisting of one polyion and two counterions and -+-for the other case. Secondly, we adopt the renormalized HNC-like closure rela- tions used earlier,2.3 6c$\(r) = g$)(r) -q"db(r)-Gz$\(r) -1 Sc;;b,(r) = Gz;;b,(r)[g"obo(r) -11 dcll,\(r) = 6z",(r)[g"db(r) -13 6cyb1((r) = [6z';b,(r)6z",b,(r) +Z"lbll(r) +(1 -aab)f$)-(r) + Kb(r)lg"obo(r) -6z?l(r) (9) where s"db(r)= expC-BU!$(r) + 6z"dbW +q"db(r)I fy)-(r) = exp[ -fiUy)-(r)] -1 and U!j(r) = Ub")(r)-Uf!(r) The functions account for trimer formation.For asym- metric binary mixtures, under the superposition approx-imation for trimers, we find qb(r) = (l -6ab)f$'-(r)[pt +p: I:(r)]/r + (l -6ac)p', 12(r)/r (10) C and the monomer densities pof and pi are determined self- consistently from the equations p+ = pof +4~pofpiJ+4~(pof)~pkJ++2npL(p,)'J-p-= pi +47~pofpiJ+47~(pi)~p;J-+27r~i(p;)~J+ (11) The division energy, U,, in eqn. (3) is determined by mini- mizing the Helmholtz free energy, given by fl(A -A("))/V= C [pa ln p:/pa +(pa -p",/2 a +nn(p"02 1(l -6ab)Pb,Ja1 (12)b where A'") is the Helmholtz free energy of the reference system with the interaction U$)(I-).~The fractions of dimers and trimers are obtained from the following expressions x2 = 8VofPOJ/Pt 9 x3(+ -+)= 64PL)2POJ+/pt and x3(-+-) = 64PO)2PLJ-/Pt J.CHEM. SOC. FARADAY TRANS., 1994, VOL. 90 where pt = p+ + p-. The explicit expressions for I:@), I&), J" and J in eqn. (10)-(13) are given by Wang and Ha~met.~ Eqn. (6) and (9) are solved numerically with 2048 grid points and a mesh 6r = 0.235 A by Picard iteration. The monomer densities, p;, are obtained iteratively from eqn. (11). The details of the algorithm are given by Kalyuzhnyi et aL2 and Ichiye and Haymet.14 The complete pair correlation functions are given by We note here that the monomer fractions for species a, x: = p",p" describe the proportion of 'bare' ions that exist in solu- tion.The input to the theory is the potential energy given by eqn. (1) and the total number density pa for species a. The theory predicts the structural and thermodynamic properties and the fraction of monomers, dimers and trimers. 4. Predictions for Asymmetric Electrolyte Thermodynamics, Dimer and Trimer Concentrations The reduced excess energy, E = EJp, kT, and the reduced osmotic pressure, (6, = P/ptkT, calculated via the virial route and the reduced optimal division energy, U,/kT, are collected in Table 2.The thermodynamic data, together with the HNC results, are compared with the MC results for + 1 : -10 elec- trolyte by Vlachy et aL4 and the RHNC results for +1 : -20 electrolyte by Linse.' The agreement between the present pre- dictions and the MC, RHNC results are good, although not all the thermodynamic quantities are improved over the results of HNC. A similar fact has been observed by Duh and Haymet" for a 2 : 2 electrolyte study by including 'bridge functions' into the cluster expansions for the pair correlation functions. The fractions of 'n-mer', x, (n = 1, 2, 3), are displayed in Table 3. Unlike the situation for 2 : 2 and 1 :3 model electro- lyte~,~.~where the trimer fractions are always lower than that of the dimers, the fraction of + -+ trimers is comparable to the dimer fraction for +1 : -10 electrolyte and even higher than the dimer fraction for +1 : -20 electrolyte, owing to the high asymmetries in both charge and size.Moreover, the monomer fractions for species a presented in Table 3 show that only 16% of polyions remain non-associated in the + 1 : -20 electrolyte. The -+ -trimer fractions are not displayed in Table 3 since they are four or seven orders of magnitude less than the + -+ trimer fractions for + 1 : -10. and +1 : -20 electrolytes, respectively. 1247 Complete Pair Correlation Functions and Partial Structure Factors It is well known that the HNC approximation produces unphysical descriptions for the structure and critical behav- iour of model electrolytes in the strong-coupling and low- concentration regime.'6J More specifically, the HNC approximation overestimates both counter-ion-poly-ion and counter-ion+ounter-ion correlations for highly asymmetric model electr~lytes.~~~~~~'~ This, in turn, leads to the shift of the first peak of the polyion-polyion correlation function to smaller separations than in the results from sim~lations.~~'~ Apparently, these deficiencies of HNC are because of its in- ability to account for ion association.'7 The predictions for the complete pair correlation functions and the partial structure factors from the present theory, together with the HNC results, for the + 1 : -10 and + 1: -20 model electrolytes are now stated.Fig.2-4 display the results for the + 1 : -10 electrolyte. Our associating inte- gral equation theory yields an improved counter-ion-counter-ion correlation function g+ +(r)shown by the solid line in Fig. 2. Fig. 3 shows that the peak of the counter-ion- 1.o 9++ 0.5 0.0 r/A Fig. 2 Complete pair correlation function for pairs of counter-ions for the + 1 : -10 electrolyte at c, = 0.08299 mol dm-'. (-) Present theory for the optimized division energy, U,/kT = -4.85, (---) HNC approximation; (0)MC simulation result by Vlachy et a1.4 Table 2 Thermodynamic quantities d = A,,/NkT, E = E,,/p, kT and & = P/p,kT for + 1 : -10 and + 1 : -20 model electrolytes 1: -10 0.08299 1: -20 0.0996 a Not available. Table 3 Monomer fractions x: + 1 : -20 model electrolytes 1: -10 0.082 99 1: -20 0.009 96 this work -4.85 -1.703 -2.126 0.625 aMC --.2.14 f0.05 0.625 f 0.02 HNC -1.699 -2.143 0.635 this work -6.97 -2.709 -3.152 0.596 RHNC -a -3.152 0.609 HNC -2.708 -3.161 0.606 = &/pa, for species a and the fractions of 'n-mers', x, E (pT-+ p,--l)/pt (n = 1, 2, 3), for + 1 : -10 and -4.85 0.907 0.327 0.854 0.075 0.071 -6.97 0.935 0.162 0.919 0.039 0.042 J.CHEM. SOC. FARADAY TRANS., 1994, VOL. 90 8 6 g+-4 2 0 r/A Fig. 3 Complete pair correlation function for pairs of unlike ions for the +1 : -10 electrolyte at c+ = 0.08299 mol dmV3. Lines have the same meaning as in Fig. 2. polyion correlation function, g + -(r), which is overestimated by HNC, is decreased in the present work, although a little severely as is the case for the +1 : -15 model electrolyte examined using the HNC + B approximation.' The polyion-polyion correlation function, g --(r), from the present theory, exhibited in Fig.4, is very close to the result of the HNC approximation. Bearing in mind that the unlike ions can at most form trimers in the present framework, we now turn to the +1 : -20 case. Fig. 5-7 present the complete pair correlation functions for this strongly coupled electrolyte. Fig. 5 shows that the counter-ion-counter-ion pair correlation function, g+ +(r), is improved slightly over the HNC approximation result. The primary maximum in the counter-ion-polyion pair correlation function g+ -(r), shown in Fig.6, is still underestimated as for the + 1 : -10 electrolyte. The polyion-polyion pair correlation function, g --(r), is again essentially overlapping with the HNC result, as displayed in Fig. 7. The partial structure factors, So&), for the two electrolytes under consideration are given in Fig. 8 and 9, with the HNC 1.5 1.o €I--0.5 0.0 rlA Fig. 4 Complete pair correlation function for pairs of polyions for the + 1 : -10 electrolyte at c, = 0.08299 mol dmW3. Lines have the same meaning as in Fig. 2. 1.5 1.o 9++ 0.5 0.0 __ 10 20 30 40 rlA Fig. 5 Complete pair correlation function for pairs of counter-ions for the + 1 : -20 electrolyte at c, = 0.0996 mol dm-3. (-) Present theory for the optimized division energy, U,/kT = -6.97; (---) HNC result; (0)result from an RHNC procedure by Linse." results also included.The results reported above suggest that ion association is still not being fully addressed within the present theory, in which unlike ions can only be associated into dimers and trimers. Further discussions to elucidate this point are given below. Partial Pair Correlation Functions One important concept in Wertheim's two-density integral equation theory is 'monomer depletion',' that is the process whereby the strong attraction between particles leads to association thus decreasing monomer densities. Mathemati-cally, the difficulty encountered by the HNC approximation for strongly coupled systems is due to the deep well in the attractive potential, or equivalently the large value of the Mayer function.Unlike the other routes used to unravel this problem before, such as including the 'bridge functions' in the HNC + B approximati~n~~'~ and introducing a three- body potential in HNC+3 theory," the depletion of mono- r/A Fig. 6 Complete pair correlation function for pairs of unlike ions for the + 1 : -20 electrolyte at c, = 0.0996 mol dm-'. Lines have the same meaning as in Fig. 5. J. CHEM.SOC. FARADAY TRANS., 1994, VOL. 90 0.0I,,,,,,,,,,,,B I I I 10 20 40 60 80 100 120 140 r/a Fig. 7 Complete pair correlation function for pairs of polyions for the + 1 : -20 electrolyte at c, = 0.0996 mol dm-'. Lines have the same meaning as in Fig. 5. I I I I I I , I 0.0 0.1 0.2 0.3 0.4 0.5 q/A-' Fig. 8 Partial structure factors for the +1 : -10 electrolyte at c+ = 0.08299 mol dm-'.(-) Present theory with the optimized division energy, U,/kT = -4.85; (---) HNC results. I 1 1 I I I I I-0.5I I 0.0 0.1 0.2 0.3 0.4 0.5 9lA-l Fig. 9 Partial structure factors for the + 1 : -20 electrolyte at c, = 0.0996 mol dm-'. (-) Present theory with the optimized division energy, U,/kT = -6.97; (---) HNC results. 1249 mers serves to suppress the catastrophe due to a large value of the Mayer function in Wertheim's two-density integral equation formalism.' To illustrate the importance of monomer depletion, we present the partial pair correlation functions, g$'(r), and the corresponding complete pair corre- lation functions, gob(r),for the + 1 : -10 electrolyte in Fig.10-12. Fig. 10 shows that the counter-ion-counter-ion correlation is mainly determined by the correlation between non-associated counter-ions g&+(r). Monomer depletion results in a decrease of g&,+(r) and the consequent improvement of the complete correlation function, g + +(r),over the HNC result. The partial correlation function, g&+(r), would be further reduced if more monomers are depleted. The results in Fig. 11 imply that the principal peak in the complete counter-ion-polyion pair correlation function, g+ -(r), arises predominantly from the contribution of the partial pair correlation between associated unlike ions, g:l-(r), which will be enhanced as more ion clusters are 0 10 20 30 40 50 rlA Fig. 10 Partial pair correlation functions for the + 1 : -10 electro-lyte: (-4 g++(r),(----) 90+0+('),(---) 90+l+(r),(-*-I 9:m 10 85 0 0 5 10 15 20 25 r/A Fig.11 As Fig. 10, but for the counter-ion-polyion correlation: (. . * *) 9;0-(') 0 20 40 60 80 100 120 r/A Fig. 12 As Fig. 10, but for polyions formed, i.e. more monomers are depleted. Fig. 12 shows the case for the complete and partial polyion-polyion pair corre- lation functions, where g&,-(r) plays the main role in deter- mining g --(r). 5. Conclusion Ion association in + 1 : -10 and + 1 : -20 model electro- lytes is investigated by Wertheim’s associating liquids integral equation formalism. It is predicted that only a small portion of ‘bare’ polyions exist.The effect of monomer depletion on the correlations in the + 1 : -10 electrolyte is analysed by examining the partial pair correlation functions. However, comparisons with the ‘exact’ results reveal that ion associ- ation is still not addressed sufficiently in the present formal- ism, in which ions can form at the most trimers. We believe that it is necessary to consider at a minimum the formation J. CHEM. SOC. FARADAY TRANS., 1994, VOL. 90 of tetramers, consisting of one polyion and three counterions, to improve the present formalism. Graphically, this means that four-particle diagrams will appear in the density equa- tions and closure relations. This research was supported by the Australian Research Council (ARC) (grant No.A29131271). A.D.J.H. acknow- ledges gratefully many helpful conversations on this topic with Prof. V. Vlachy, Dr. Yu.Kalyuzhnyi, Prof. M. Holovko and Mr. M. J. Booth. References 1 M. S. Wertheim, J. Stat. Phys., 1984,35, 19, 35. 2 Y. V. Kalyuzhnyi, M. F. Holovko and A. D. J. Haymet, J. Chem. Phys., 1991,95,9151. 3 J. Wang and A. Haymet, J. Chem. Phys., 1994,100,3767. 4 V. Vlachy, C. H. Marshall and A. D. J. Haymet, J. Am. Chem. Sac., 1989,111,4160. 5 J. ReSEiE, V. Vlachy and A. D. J. Haymet, J. Am. Chem. SOC., 1990,112,3398. 6 Y. N. Vorobjev, J. A. Grant and H. A. Scheraga, J. Am. Chem. SOC.,1992, 114, 3189. 7 A. Nyberg, D. E. Smith, L. Zhang and A. D. J. Haymet, J. Chem. Phys., 1994, submitted. 8 D-M. Duh, D. Perera and A. D. J. Haymet, J. Chem. Phys., 1994, in preparation. 9 P. Linse, J. Chem. Phys., 1990,93, 1376. 10 P. Linse, J. Chem. Phys., 1991,94,3817. 11 P. J. Rossky, J. B. Dudowicz, B. L. Tembe and H. L. Friedman, J. Chem. Phys., 1980,73,3372. 12 H. C. Andersen, J. Chem. Phys., 1973,59,4714. 13 H. C. Andersen, J. Chem. Phys., 1974,61,4985. 14 T. Ichiye and A. D. J. Haymet, J. Chem. Phys., 1988,89,4315. 15 D-M. Duh and A. D. J. Haymet, J. Chem. Phys., 1992,9,7716. 16 H. L. Friedman, Annu. Rev. Phys. Chem., 1981,32,179. 17 J. S. Hsye, E. Lomba and G. Stell, Mol. Phys., 1992,75, 1217. 18 B. B. Laird, J. Wang and A. D. J. Haymet, Phys. Rev. E, 1993, 47.2491. Paper 3/07172C; Received 6th December, 1993