Acta Chini. Slov. 1999, 46(A), pp. 531-541 Structure and thermodynamics of a two-dimensional model of electrolyte SOLUTIONS x T. Urbič and V. Vlachy Faculty of Chemistry and Chemical Technology, University of Ljubljana, Aškerčeva 5, 1000 Ljubljana, Slovenia. Abstract Monte Carlo and integral equation results in the hypernetted-chain approximation are presented for the classical two-dimensional Coulomb fluid above the critical temperature for the Kosterliz-Thouless transition. The equilibrium structure and thermodynamics were studied as a function of the concentration of the fluid. The Monte Carlo results were used to assess the accuracy of the corresponding predictions of the hypernetted-chain theory. Comparison of the pair distribution functions and thermodynamics shows that the theoretical results are in semiquantitative or better agreement with the simulations. Introduction The classical two-dimensional (2D) Coulomb fluid with logarithmically interacting charges [1] is relevant to an understanding of many different physical systems. The 2D Coulomb fluid consisting of an equal number of positive and negative charges serves as a model for interactions in thin superconducting and superfluid He films [2]. At higher temperatures the individual ions move freely and the fluid conducts electrical current. As the temperature is decreased, a critical temperature occurs at which neutral pairs of particles are formed. This is the well known Kosterliz-Thouless (KT) [3] transition between the conducting and dielectric phase. Dedicated to professor Drago Leskovšek on occasion of his 80*" birthday 531 A fluid of charged disks interacting via a logarithmic pair potential also has some relevance for solutions of rigid cylindrical polyelectrolytes. In recent papers Kholo-denko and Beyerlin [4] and Levin [5] discuss the relation between the Manning counterion condensation [6] and the Kosterlitz-Thouless phase transition. For us this model is of interest in connection with recent calculations on 2D model of water solutions [7, 8]. An extension of these calculations to study a solvation of ionic solutes would require a solution of the integral equation for a mixture of water-like molecules and ions. In this paper we present integral equation theory results for a system of charged hard disks. The system as a whole is electroneutral: one half of the particles have positive charge and the other half carry negative charge. The solvent is treated as a dielectric continuum. This model has been examined using the Monte Carlo simulation in several papers [9, 10, 11] but has not yet, to our knowledge, been studied by an integral equation technique. In this contribution the Ornstein-Zernike equation for the model in question is solved in the hypernetted-chain (HNC) approximation. We studied the model fluid in a broad concentration range above the KT critical temperature where the HNC integral equation can be solved. As a result the pair correlation functions between various species (not studied in [9, 10, 11]) and the excess internal energy and pressure were obtained. These results are complemented by canonical Monte Carlo calculations for the same model to assess the validity of the HNC approximation. The model and theory The system considered here is composed of N+ hard disks carrying a positive charge q+ and AL disks having a charge g_. The diameter of the particles is a and they are embedded in a continuous dielectric with permittivity s. The pair interaction potential for two particles separated by distance r^ is defined by r _^iLLmrü r.. > 0 «y(ry)= 2wLL° * X3~/ (1) J [ oo T'ij < a Note that the zero of the potential is chosen at r = a. If we introduce reduced units as r* = L z = sttk anc[ u* = Bi where Ah = » e" T and L is the characteristic length, we can write equation (1) in a simpler form <(r*,) = ( -ZiZiinriJ ^-J (2) The HNC approximation consists of two equations. One is the Ornstein-Zernike equation, which correlates the total correlation function, h(r) = g(r) — 1, and the 532 direct correlation function, c(r) [12] M12) = Q,(12) + J2 f Cik(13)pkhkj(32)d(3), (3) k J where pk is the number density of the component k. The second equation is given by [12] ^•(12) = htj(12) - In [1 + M12)] - ßuij(12) + 5*, (12). (4) The HNC approximation sets the unknown bridge function -By(r) to zero for all distances r. A necessary first step toward numerical solution of the Ornstein-Zernike equation for ionic systems is the renormalization procedure. For three-dimensional ionic fluids, the renormalization scheme is well known (see for example, [13]). In the present work a similar technique was applied. The total potential is separated into uij=uij+ulip (5) where u*j and «*J are the short and long-range parts defined as ( 0 r* > 1 u'ij(rij) = -zizjlnr*J. (7) We define the short-range total and direct correlation functions to be Kj = K - i>a (8) cij = cij ~ Viji \") respectively, where çV, is the long-range limit of the direct correlation function ij = -ßufj (10) and ipij is the screened 2D Coulomb potential % = -f^oK). (ii) K0(r) is the modified Bessel function, T* = ^- is the reduced temperature, k = fL X)i Pizi is the screening parameter and p* = pa2 is the reduced density. We make the further definitions V'1 = p'1 - (j). (12) 533 and T = hs - cs. (13) By substituting Eqs. (8)-(13) into Eq. (3) we obtain the final expression for the renormalized Ornstein-Zernike equation which in the matrix form written is as f = -cs + p"1 (l - Fcj"1 VeVp-1. (14) The renormalized closure has the following form elf = exp (-pXj + Tij + ipij) - 1 - Tij - tki- (15) Solution of the OZ equation (14) together with the HNC closure conditions (15) has been obtained using the direct iteration method. Forward and inverse Bessel-Fourier transforms, needed to couple the correlation functions in real and Fourier space, have been carried out on a logarithmic scale using the method developed by Taimen [14]. The computations were performed on a grid with 512 to 2048 integration points, depending on concentration. The numerical accuracy was monitored by checking the electroneutrality condition: the zero moment defect was less than 0.001 in all calculations reported below. The Monte Carlo simulations were performed at constant volume and temperature with 256 (or in some cases 512) particles in the system (square), using the standard Metropolis algorithm [12]. From two to four million configurations, after an equilibration run of several hundred thousand, were enough to obtain reliable statistics. To avoid effects due to the finite size of the system the Ewald summation method [15] was used. Results and discussion The model can be characterized by two nondimensional parameters, the reduced density p* and temperature T*. The calculations presented here apply to four temperatures T*= 0.5, 1.0, 2.0, and 4.0 and to the concentration range from p* = 0.00001 to 0.5. Computer simulations [9, 10, 11] disagree on the exact location of the critical point. According to [10] the critical temperature is 0.056, while the most recent simulations of Lidmar and Wallin [11] suggest the value of 0.032. In any case, our calculations apply to conditions well above the critical temperature for the conducting phase-dielectric phase transition. We consider first the results for the pair distribution function (pdf), g(r*), shown for various temperatures and densities in Figures 1-6. It is interesting that no structural 534 1.6 1.4 1.2 1 E 0.8 0.6 ft IQ 0" a b 0.4 0.2 0 O O O oo1 0 X3J a ^-©^©-G-G-e-e-e-G^e-e-o-e^©-^' Figure 1: Pair distribution function for like-charge ions at T* = 1.0 and p* = 0.5. Monte Carlo simulation data are presented by symbols and the continuous curve denotes the HNC results. "5? Figure 2: Figure 1. Pair distribution function for unlike-charge ions. The parameters as for 535 Figure 3: Pair distribution function for like-charge ions; T* = 2.0 and p* = 0.02. "5? Figure 4: Pair distribution function for unlike-charge ions; T* = 2.0 and p* = 0.02. 536 1.2 1 0.8 ^ 0.6 0.4 0.2 O çyQ& &Qxy&® ^°^ €>^^&^ ^©^ &<5 <9 é> &' oya& i 0 © ©' Ooe1- _L 10 15 20 Figure 5: Pair distribution function for like-charge ions; T* = 0.5 and p* = 0.02. Figure 6: Pair distribution function for unlike-charge ions; T* = 0.5 and p* = 0.02. 537 information was presented in previous papers concerned with computer simulations [9, 10, 11]. The HNC results in these figures are shown by continuous lines and the Monte Carlo results are denoted by symbols. In Figures 1 and 2 we present the like-ion, g(r*), and the unlike-ion pair distribution functions, g(r*) at p* = 0.5 and T*= 1.0. The two methods give pdfs which are in reasonably good agreement for these values of the parameters. The same holds true for the next pair of figures. In Figures 3 and 4 we show the two pair distribution functions at the lower concentration p* = 0.02 and for T*= 2.0. Figures 5 and 6 present the structure for p* = 0.02 and T*= 0.5. We note in passing that for these conditions we are close to the limits of the solvability of the HNC theory; the lowest reduced temperature where we obtained a convergent solution at this concentration was T*= 0.3. In the next two figures the thermodynamic parameters are presented as a function of the fluid density. The excess internal energy and the osmotic coefficient reported below were evaluated using standard thermodynamic equations [16]. The expression for the excess internal energy per particle is uex = ^ Y, PiP*j / ruij9ijdr, (16) Pt ij ¦> where p\ is total density. Because of the shape of the pair potential the excess energy approaches infinity for zero concentration p* = 0 [9]. The osmotic pressure P was obtained using the virial theorem [16]. An important quantity is the osmotic coefficient defined as (f) = P/pkBT. It is easy to show that the general expression for the osmotic coefficient (f) equal to (/) = 1~ 2pJF ? PIP^S f2OA - 0yU)) (17) reduces in the case of a 2D Coulomb fluid to the following formula ^=1-7^ + ^EpI^(1)- (18) P t i, j For the point ions (a = 0) we correctly obtain 0 = 1 — 1/4T*; this limiting result was first derived Salzberg and Prager [16]. The results for thermodynamic quantities are shown in Figures 7 and 8. First we present the concentration dependence of the excess internal energy at three different temperatures. Again, the results presented by lines apply to the HNC calculation and symbols represent the Monte Carlo data. We compared our calculations with the results of previous simulations. For example, our result for the reduced excess energy at T* = 0.5, p* = 0.0636 is 0.299, while the result of Caillol and Levesque 538 2.5 2 1.5 b lili i >'. \ A.\ x \ \ "—. i i i — —t------------------ --------- 0.05 0.1 0.15 0.2 0.25 P* Figure 7: The excess energy per particle as a function of p* for different temperatures T* = 4.0, 2.0, 1.0 and 0.5 going from the upper curve to the bottom curve. Monte Carlo results lie on the corresponding HNC curves and are shown by symbols (crosses). o / \ ' / , / 2.5 y x / • / 2 y y * y „ ' 1.5 ^ ,,--' + " 1 --'"'+----"'''" - .+- >-- 0.5 n i i i i i 0.1 0.2 0.3 P* 0.4 0.5 Figure 8: The osmotic coefficient as a function of p*. The upper curve is calculated for T* = 0.5 and the lower curve for T* = 1.0. Monte Carlo results for the reduced temperature 1.0 are presented by symbols. 539 [9] is 0.303. The HNC result for this quantity is 0.287. The agreement between the two simulations is reasonable; note that these authors [9] used a substantially different simulation procedure. The numerical error in our simulation is estimated to be about 5%. Next, in Figure 8, the results for osmotic coefficient are presented. Because of the logarithmic form of the potential the osmotic coefficient is a monotonically increasing function of the electrolyte concentration. The predictions of the HNC theory are in agreement with the machine calculations for the values of parameters examined in this work. The HNC approximation seems to be less accurate at higher volume fractions of the fluid model. This is consistent with our earlier results for 3D models of electrolyte solutions [19]. Conclusions The properties of the two-dimensional Coulomb fluid are of interest for several problems in physics and chemistry. In the present paper the integral equation theory in the hypernetted-chain approximation was applied to the model fluid. Parallel simulations, performed utilizing the Monte Carlo method, reveal the usefulness of the hypernetted-chain theory. This approximation is able to predict the pair distributions and thermodynamic parameters to a very reasonable degree of accuracy. This observation is in agreement with previous applications of the hypernetted-chain theory to Coulomb systems in three dimensions (see for example [17]). Note that the HNC approach requires at least an order of magnitude less computer time than the Monte Carlo method. A major disadvantage of the HNC approximation is that the theory does not yield convergent solutions for low temperatures. One way to extend the region of solvability of the HNC approximation is to construct a theory (see [18] and references therein) based on Wertheim's two-density formalism. This improvement is currently under investigation. Acknowledgements This work was supported by the Slovene Ministry of Science and Technology and in part by the US-Slovene Science and Technology Joint Fund. References [1] P. Minhagen, Rev.Mod.Phys. 1987, 50, 1001-1066. [2] D. R. Nelson and J. M. Kosterliz, Phys. Rev. Lett. 1977, 39, 1201-1205. [3] J. M. Kosterliz and D. J. Thouless, J. Phys. C 1973, 6, 1181-1203. 540 A. L. Kholodenko and A. L. Beyerlin, Phys. Rev. Lett. 1995, 74, 4679-4681. Y. Levin, Physica A 1998, 257, 408-412. G. S. Manning, J. Chem. Phys. 1969, 51, 924-933. K. A. T. Silverstein, A. D. J. Haymet, K. A. Dill, J. Am. Chem. Soc. 1998, 120, 3166-3175. T. Urbič, V. Vlachy, Yu. V. Kalyuzhnyi, N. T. Southall and K. A. Dill, to be submitted. J. M. Caillol and D. Levesque, Phys. Rev. B 1986, 33, 499-509. G. Orkoulas and A. Z. Panagiotopoulous, J. Chem. Phys. 1996, 104, 7205-7209. J. Lidmar and M. Wallin, Phys. Rev. B 1997, 55, 522-530. H. L. Friedman, A Course in Statistical Mechanics, Prentice-Hall, Inc., Englewood Cliffs, New Jersey 07632, (1985). T. Ichiye and A. D. J. Haymet, J. Chem. Phys. 1988, 89, 4315-4324. J. D. Taimen, J. Comp. Phys. 1978, 29, 35-48. J. W. Perram and S. W. de Leeuw, Physica A 1981, 109, 237-250. A. M. Salzberg and S. Prager, J. Chem. Phys. 1963, 38, 2587. V. Vlachy, Annu. Rev. Phys. Chem. 1999, 50, 145-165. Yu. V. Kalyuzhnyi and V. Vlachy, Chem. Phys. Lett. 1993, 215, 518-522. J. Reščič, V. Vlachy, L. B. Bhuiyan and C. W. Outhwaite, Mol. Phys. 1998, 95, 233-242. Povzetek V članku so predstavljeni rezultati, ki so bili dobljeni z Monte Carlo simulacijo in z rešitvijo integralskih enačb v HNC približku za klasičen dvodimenzionalni Coulombski plin nad kritično temperaturo Kosterlitz-Thoulessovega faznega prehoda. Termodinamiku in strukturo sva preučevala kot funkcijo koncentracije plina. Rezultate Monte Carlo simulacije sva uporabila za ugotovitev točnosti HNC teorije. Primerjava parskih porazdelitvenih funkcij in termodinamskih količin kaže, da se teoretični (HNC) rezultati semikvantitativno ali pa celo bolje ujemajo z rezultati dobljenimi s simulacijo. 541