Simulation study of the link between molecular association and reentrant miscibility for a mixture of molecules with directional interactions

VOLUME 82, NUMBER 26 PHYSICAL REVIEW LETTERS 28JUNE 1999 Simulation Study of the Link between Molecular Association and Reentrant Miscibility for a Mixture of Molecules with Directional Interactions Lowri A. Davies,1George Jackson,2and Luis F. Rull3 1Department of Chemistry, University of Sheffield, Sheffield S3 7HF, United Kingdom 2Department of Chemical Engineering and Chemical Technology, Imperial College of Science, Technology and Medicine, Prince Consort Road, London SW7 2BY, United Kingdom 3Departamento de Fı ´sica Atómica, Molecular y Nuclear, Universidad de Sevilla, Apartado 1065, Sevilla 41080, Spain (Received 20 January 1999) The reentrant liquid-liquid miscibility of a symmetrical mixture with highly directional bonding interactions is studied by Gibbs ensemble Monte Carlo simulation. The resulting closed loop of immiscibility and the corresponding lower critical solution temperature are shown to be a direct consequence of the dramatic increase in association between unlike components as the temperature is lowered. Our exact calculations for an off-lattice system with a well-defined anisotropic potential confirm the findings of previous theoretical studies. [S0031-9007(99)09439-9] PACS numbers: 64.60.–i, 61.20.Qg, 05.20.Jj Reentrant phase behavior is one of the most characteristic and beguiling signatures of hydrogen bonding in simple and complex fluids alike. The so-called closedloop liquid-liquid immiscibility found in a wide variety of fluid mixtures is a common example of the relationship between hydrogen-bonding association and reentrant behavior: here the region of liquid-liquid immiscibility is bounded above and below by upper and lower critical solution temperatures (UCST and LCST) (see Refs. [1,2] for general reviews). Closed-loop immiscibility was first reported for the nicotine-water mixture at the turn of the 20th Century [3]. Other well-known examples of this type of reentrant behavior include aqueous solutions of butan-2-ol [4], alkylpyridines [5], and alkylpolyoxyethylene nonionic surfactants [6], as well as a number of polymer solutions [7,8]. The common feature of these systems is the possibility of hydrogen bonding between the two components of the mixture. An increase in temperature overcomes the unfavorable enthalpic contribution (incompatibility of the two components) leading to miscibility above the UCST. The low-temperature miscibility below the LCST can be explained in terms of hydrogen bonding between unlike components which becomes the dominant feature of the free energy of the system, and results in a favorable heat of solution. One should note that in the case of polymer solutions the principal feature responsible for the LCST is the large difference in size between the polymer molecule and the solvent. Hirshfelder et al. [9] were the first to suggest directional attractive forces as a direct cause of the LCST. This view is supported by the theoretical predictions of Barker and Fock [10], who analyzed a simple lattice model with interactions that depend on the relative orientations of the molecules using the Bethe-Guggenheim quasichemical approximation; an appropriate choice of parameters for the unlike directional interactions gives rise to closedloop behavior. A more rigorous treatment of similar lattice models, using appropriate decoration transformations to the equivalent Ising models of ferromagnets [11] or alternatively using the Migdal-Kadanoff approximations of renormalization group (RG) theory [12], has also been applied to study closed-loop immiscibility in systems with highly asymmetrical interactions (e.g., see Refs. [13–18]). The results of the more recent calculations are in reasonable quantitative agreement with the asymmetric coexistence curves found in real systems, at least for pressures close to atmospheric. In addition to these analytical approaches, reentrant miscibility in similar lattice models has been studied with Monte Carlo techniques for both bulk [19] and confined systems [20]; as for the exact analytical studies, e.g., Refs. [13,15], symmetrical closed loops are seen, with the additional discovery of low-temperature self-assembled supermolecular structures. The lattice models are usually decorated lattices in which the different sides of the molecule are treated differently (directional) so that the contact energies of interaction with its neighbors depends on the relative orientation of the molecules. In all of these studies the reentrant behavior is a direct consequence of the directional interactions used to mimic hydrogen bonding. It is clear from the ensuing discussion that most of the theoretical studies of closed-loop immiscibility have been within the framework of lattice models. Although the lattice approach provides the correct qualitative description of fluid phase behavior, the representation of the structure is a vast oversimplification. In most instances only nearest neighbor interactions are taken into account. Furthermore, the density of a lattice model is difficult to relate to that of a real fluid, and hence the pressure dependence of the phase equilibria is unrealistic. Continuum models on the other hand provide a more adequate description of a fluid. Various phenomenological continuum approaches have been used to study closed-loop immiscibility. Here, the highly anisotropic interactions are treated in terms of angle-averaged intermolecular potentials which have the appropriate temperature dependence. Reatto and Tau [21] 0031-9007兾99兾82(26)兾5285(4)$15.00 © 1999 The American Physical Society 5285 VOLUME 82, NUMBER 26 PHYSICAL REVIEW LETTERS 28JUNE 1999 studied systems of spherical molecules with a temperaturedependent potential of mean force which includes solvent effects by accounting for the decrease in the size of the hydration shell as the temperature is increased; the mean spherical approximation was applied to describe the LCSTs of model micellar solutions. Goldstein [22] used a Flory-Huggins free energy expression together with a simple phenomenological term for the hydrogen-bonding interaction to describe closed-loop behavior in aqueous surfactant solutions. In fact an effective intermolecular interaction with the correct temperature dependence can always be found to reproduce reentrant behavior [23]. The problem with this type of phenomenological approach is that the exact form of the temperature dependence is unknown, and hence the description is merely qualitative and difficult to relate to a specific system. An obvious alternative is to examine continuum models with well-defined directional attractive interactions. The Wertheim [24] approach offers an appropriate statistical mechanical framework to examine these associating systems. Such an approach has been used within an augmented van der Waals perturbation theory to describe closed-loop immiscibility in mixtures of molecules with directional attractive forces [2]: the model involves a symmetrical mixture of equal sized spherical molecules with unfavorable unlike mean-field interactions (which lead to liquid-liquid phase separation at moderate temperatures), and with additional directional bonding interactions between the unlike species (which leads to association and miscibility at low temperatures). This is essentially the first continuum study of closed-loop immiscibility for a system with a well-defined directional intermolecular potential. Although the preliminary theoretical study was for a rather unrealistic prototype system, the approach has been extended to provide a quantitative description of mixtures such as aqueous solutions of butanol, butoxyethanol, and alkylpolyoxyethylene nonionic surfactants [25,26]; it is particularly gratifying to find that the correct highpressure behavior is predicted for these systems with these relatively simple continuum models. Despite the large body of work in this area (the principal studies have already been mentioned), the link between association and the closed-loop behavior has not been explicitly shown for a specific model. In this Letter we present the results of Gibbs ensemble Monte Carlo simulation for the liquid-liquid phase equilibria of a symmetrical binary mixture with directional unlike attractive interactions, which is similar to the system studied theoretically in Ref. [2]. We aim to provide the first exact calculations of closed-loop immiscibility for an off-lattice system with a well-defined Hamiltonian. As well as confirming the results of the earlier theoretical findings, the focus of this paper is to show the direct unequivocal link between the association of the unlike species and the occurrence of the LCST. The simulation results of this prototype system should also provide a basis for future studies of the LCST, such as the nature of the critical exponents and wetting behavior. We study a symmetrical square-well binary mixture in which the unlike isotropic interaction is purely repulsive: the spherical core diameters sare equal 共s苷s11 苷 s12 苷s22兲, as are the ranges of the interactions 共l苷 l11 苷l12 苷l22 苷1.5兲and the well depths of the like interactions 共´苷´11 苷´22兲, while the unlike well depth is zero 共´12 苷0兲. This system is seen to exhibit large regions of liquid-liquid immiscibility [27] due to the unfavorable interactions between unlike components. The inclusion of directional unlike attractions is necessary for the system to exhibit low-temperature miscibility. This is achieved by adding a single off-center and short-ranged square-well association site on each component to mimic the hydrogen-bonding interaction: the sites are placed halfway between the center and surface of the spherical core, have a range of lhb 12 苷0.55, a well depth of ´hb 12 苷 14.5´, and only unlike site-site association is allowed. The Gibbs ensemble Monte Carlo (GEMC) simulation technique [28,29] has become one of the most versatile methods to study phase equilibria in pure fluids and mixtures [30]. We use the standard GEMC technique in the isothermal-isobaric ensemble 共NPT兲as described in Ref. [27] although additional reorientational moves are now required. Because of the small volume of the sitesite overlap region, large system sizes (N苷1000 and 1728, with equal numbers of both components) and long simulations (⬃2.5 3105GEMC cycles) are required to obtain good statistics. A typical phase diagram for the model associating system is shown in Fig. 1(a) as a temperature-composition 共Tx兲isobaric slice. The closed-loop region of liquidliquid immiscibility is clearly discernible; the system is totally miscible above a UCST 共Tⴱ c⬃1.92兲and below an LCST 共Tⴱ c⬃1.56兲. It is clear from Fig. 1(a) that the system mixes below the LCST; the two boxes in the GEMC simulations become identical for temperatures slightly above the LCST. For larger values of the sitesite energy ´hb, the unlike interactions are strong enough for the system to mix completely; in the case of smaller values of ´hb the liquid-liquid immiscibility persists at lower temperatures [2]. An increase in the site-site interaction also increases the UCST, but the effect is less dramatic. Because of the symmetry of the intermolecular interactions, the phase diagram is symmetrical about a mole fraction of x苷1 2: the densities and the degree of association of the two coexisting phases are equal while the mole fractions in the two phases are related by x and 12x.Wedefine the degree of association in terms of Xhb, the fraction of the total number of molecules which are bonded; two molecules are considered to be associated when the off-center bonding sites overlap. The fraction of molecules bonded in the coexisting phases is depicted in Fig. 1(b). Close to the UCST about 10% of the molecules are bonded. As expected the fraction of 5286 VOLUME 82, NUMBER 26 PHYSICAL REVIEW LETTERS 28JUNE 1999 0.0 0.2 0.5 0.8 1.0 x 1.5 1.6 1.7 1.8 1.9 2.0 T* a) 0.0 0.2 0.4 0.6 Xhb 1.5 1.6 1.7 1.8 1.9 2.0 T* b) FIG. 1. (a) The temperature-composition 共Tx兲phase diagram for the symmetrical square-well mixture with unlike directional bonding interactions for a reduced pressure of Pⴱ 苷1.28 (where Tⴱ 苷kT兾´and Pⴱ 苷ps3兾´). The data points correspond to Gibbs ensemble Monte Carlo simulations of the coexisting liquid phases for N苷1000 (squares) and 1728 (circles) particles, the estimated critical point (diamonds), and the solid curve is the Wegner fit to the simulation data. (b) The fraction of the total number of molecules bonded Xhb as a function of temperature. The degree of association is shown for a single phase which, due to the symmetry of the interactions, is the same for both coexisting phases. The error bars correspond to one standard deviation. molecules bonded increases with decreasing temperature. A dramatic increase in association is seen close to the LCST, where about 50% of the molecules are seen to be bonded. This clearly indicates that the LCST and subsequent miscibility occurs when a large proportion of the molecules is associated. The fraction of bonded molecules increases more slowly as the temperature is lowered below the LCST; the large change in association appears to be closely linked to the lower critical behavior. One of the advantages of GEMC simulation is that it allows a description of phase equilibria moderately close to the critical point, at least until the critical fluctuations become larger than the box dimensions [31]. Simulations of finite systems in the canonical ensemble results in an overestimate of the critical temperature (or an underestimate for an LCST). The results of GEMC simulations can be treated as “experimental data”and analyzed using the ideas of universal critical exponents and corrections to scaling [32]. The liquid-liquid phase transition exhibited by the square-well mixture belongs to the three-dimensional Ising-type universality class. Renormalization group theory leads to scaling laws for critical phenomena which can be used to correct for the behavior away from the asymptotic critical region [33]. The resulting expansion for the mole fraction can be written as an expansion in the departure from the critical temperature t苷12T兾Tcas [27,32]: x6苷xc1C1jtjc1C2jtj1C3jtjc1D11··· 61 2共B0jtjb1B1jtjb1D11B2jtjb12D11···兲, (1) where x2and x1represent the two coexisting compositions, respectively, xcis the composition at the critical point, Ciare coefficients in the expansion of the diameter of the coexistence curve, ccharacterizes the anomaly in the diameter, D1is a renormalization group gap exponent, Biare the correction amplitudes, and bis the critical exponent. In the case of our square-well mixture we take the universal critical exponent of b苷0.325, a RG gap exponent of D1苷0.5, and include only the two leading terms B0and B1in the Wegner series to describe the liquid-liquid phase equilibria; the symmetry of the coexistence curve also means that xc苷1 2and C2苷0 (since the diameters are parallel to the temperature axis). The series for the mole fraction of the system can thus be written as x6苷 1 261 2共B0jtjb1B1jtjb1D1兲.(2) The UCST and LCST can be analyzed with two separate series, and a full description of the closed loop can then be obtained by combining the two using a hyperbolic tangent switching function. The results of such a fit of the extended scaling series to the simulation data is shown as the solid curve in Fig. 1(a); the optimum values of the critical parameters are Tⴱ c苷1.92,B0苷 2.32, and B1苷22.61 for the UCST and Tⴱ c苷1.56, B0苷3.00, and B1苷24.15 for the LCST. Clearly, our simple truncated Wegner series provides an adequate analytical description of the closed-loop envelope. Finitesize scaling (FSS) simulation techniques have recently been developed to properly describe the critical region; a full description can be found in the review by Allen [34]. We are currently undertaking a full FSS study of the LCST region to confirm universal critical behavior in this case. In conclusion, we confirm that the low-temperature miscibility in systems with directional attractive interactions is a direct result of the sharp increase in association as the temperature approaches the LCST. It is important 5287 VOLUME 82, NUMBER 26 PHYSICAL REVIEW LETTERS 28JUNE 1999 to mention that, against the generally accepted view that closed-loop immiscibility and highly directional interactions are intricately linked, this type of phase behavior has been reported for Lennard-Jones molecules [35] and attractive spheres [36], where the attractive interactions are spherically symmetrical. These studies use approximate theories to describe the thermodynamics of the systems, and the corresponding LCSTs are found at extremely low temperatures for densities which are deep in the solid region. Although it is possible that isotropic interactions could give rise to reentrance under such unphysical conditions, anisotropic interactions such as hydrogen bonding are the central feature of systems which exhibit closedloop behavior, a fact which is borne out by experiment. Here we have focused on liquid-liquid immiscibility of two phases with similar effective density which disappear at low temperatures. “Closed-loop”coexistence between phases of very different nature (e.g., vapor and liquid phases of different density) can be seen in nonassociating systems at high pressures and temperatures for systems with symmetrical interactions (see global phase behavior of Scott and van Konynenburg [37]). Reentrance can also be seen in charged-stabilized colloidal suspensions [38], which again is not caused by hydrogen bonding or molecular association; the complex interplay between Coulombic and entropic effects gives rise to the phase behavior. Reentrant phase behavior is of course not restricted to simple fluids. Reentrant liquid crystalline nematic phases can be seen as a consequence of the delicate temperature dependence of hydrogen bonding (e.g., see Ref. [39]). In fact, molecular association and dipole pairing are often used to explain reentrant behavior [40]. Hydrogen bonding offers a number of new possibilities to control the phase transitions in liquid crystalline and polymeric systems. L.A.D. thanks the EPSRC and BP for financial support. We also acknowledge support from the Royal Society, the EPSRC, the European Union, and the DGICYT (PB94-1442) for computer hardware. [1] J.S. Walker and C.A. Vause, Sci. Am. 256, 90 (1987). [2] G. Jackson, Mol. Phys. 72, 1365 (1991). [3] C.S. Hudson, Z. Phys. Chem. 47, 113 (1904). [4] W. Dolgolenko, Z. Phys. Chem. 62, 499 (1908). [5] R.J.L. Andon and J.D. Cox, J. Chem. Soc. 1952, 4601 (1952). [6] J.C. Lang and R.D. Morgan, J. Chem. Phys. 73, 5849 (1980). [7] F.F. Nord, M. Bier, and S.N. Timasheff, J. Am. Chem. Soc. 73, 289 (1951). [8] G.N. Malcolm and J.S. Rowlinson, Trans. Faraday Soc. 53, 921 (1957). [9] J. Hirshfelder, D. Stevenson, and H. Eyring, J. Chem. Phys. 5, 896 (1937). [10] J.A. Barker and W. Fock, Discuss. Faraday Soc. 15, 188 (1953). [11] M.E. Fisher, Phys. Rev. 113, 969 (1959). [12] A.A. Migdal, Sov. Phys. JETP 42, 743 (1976); L.P. Kadanoff, Ann. Phys. (N.Y.) 100, 359 (1976). [13] J. Wheeler, J. Chem. Phys. 62, 433 (1975); J. Wheeler and G.R. Andersen, J. Chem. Phys. 69, 2082 (1978); ibid. 69, 3403 (1978); ibid. 73, 5778 (1980). [14] A. Robledo, Europhys. Lett. 1, 303 (1986). [15] J.S. Walker and C.A. Vause, Phys. Lett. 79A, 421 (1980); ibid. 90A, 419 (1982). [16] R.E. Goldstein and J.S. Walker, J. Chem. Phys. 78, 1492 (1983). [17] J.-C. Lin and P.L. Taylor, Phys. Rev. Lett. 73, 2863 (1994); Phys. Rev. E 49, 2058 (1994). [18] J.-M. Romero-Enrique, I. Rodriguez-Ponce, L.F. Rull, and U. Marini Bettolo Marconi, Mol. Phys. 93, 501 (1998); Mol. Phys. 95, 571 (1998); Phys. Rev. Lett. 79, 3543 (1997). [19] M. Kotelyanskii, B. Veytsman, and S.K. Kumar, Phys. Rev. E 58, R12 (1998). [20] M. Kotelyanskii and S.K. Kumar, Phys. Rev. Lett. 80, 1252 (1998). [21] L. Reatto and M. Tau, Chem. Phys. Lett. 108, 292 (1984). [22] R.E. Goldstein, J. Chem. Phys. 84, 3367 (1986). [23] M. Gencaslan, P.H.E. Meijer, M. Keskin, and A.H.L. Levelt, J. Supercrit. Fluids 7, 107 (1994). [24] M.S. Wertheim, J. Stat. Phys. 35, 19 (1984); ibid. 35,35 (1984); ibid. 42, 459 (1986); ibid. 42, 477 (1986). [25] M.N. Garcia-Lisbona, A. Galindo, G. Jackson, and A. Burgess, Mol. Phys. 93, 57 (1998). [26] M.N. Garcia-Lisbona, A. Galindo, G. Jackson, and A. Burgess, J. Am. Chem. Soc. 120, 4191 (1998). [27] D.G. Green, G. Jackson, E. de Miguel, and L.F. Rull, J. Chem. Phys. 101, 3190 (1994). [28] A.Z. Panagiotopoulos, Mol. Phys. 61, 813 (1987). [29] A.Z. Panagiotopoulos, N. Quirke, M.R. Stapelton, and D.J. Tildesley, Mol. Phys. 63, 527 (1988). [30] K.E. Gubbins, Mol. Simul. 2, 223 (1989); A.Z. Panagiotopoulos, Mol. Simul. 9, 1 (1992). [31] B. Smit, Ph. de Smedt, and D. Frenkel, Mol. Phys. 68, 931 (1989). [32] J.V. Sengers and J.M.H. Levelt Sengers, Annu. Rev. Phys. Chem. 37, 189 (1986). [33] F.J. Wegner, Phys. Rev. B 5, 4529 (1972). [34] M.P. Allen, Computer Simulation in Chemical Physics, edited by M.P. Allen and D.J. Tildesley (Kluwer Academic, New York, 1993), pp. 49–92. [35] V.A. Mazur, L.Z. Boshkov, and V.B. Fedorov, Dokl. Akad. Nauk SSSR 282, 137 (1985); L.Z. Boshkov and V.A. Mazur, Russ. J. Phys. Chem. 60, 16 (1986); L.Z. Boshkov, Dokl. Akad. Nauk SSSR 294, 901 (1987). [36] L.V. Yelash and T. Kraska, Ber. Bunsen-Ges. Phys. Chem. 102, 213 (1998). [37] R.L. Scott and P.H. van Konynenburg, Discuss. Faraday Soc. 49, 87 (1970). [38] R. van Roij, M. Dijkstra, and J.-P. Hansen, Phys. Rev. E 59, 2010 (1999). [39] R.P. Sear and G. Jackson, Mol. Phys. 83, 961 (1994); Phys. Rev. Lett. 74, 4261 (1995). [40] P.E. Cladis, Mol. Cryst. Liq. Cryst. 85, 165 (1988). 5288