Closed-loop phase equilibria of a symmetrical associating mixture of square-well molecules examined by Gibbs ensemble Monte Carlo simulation

Closed-loop phase equilibria of a symmetrical associating mixture of square-well molecules examined by Gibbs ensemble Monte Carlo simulation Lowri A. Davies Department of Chemistry, University of Sheffield, Sheffield S3 7HF, United Kingdom George Jackson Department of Chemical Engineering and Chemical Technology, Imperial College of Science, Technology and Medicine, Prince Consort Road, London SW7 2BY, United Kingdom Luis F. Rull Departamento de Fı ´sica Ato ´mica Molecular y Nuclear, Universidad de Sevilla, Apartado 1065, Sevilla 41080, Spain 共Received 21 June 1999兲 A closed loop of liquid-liquid immiscibility for a simple model binary symmetrical mixture of square-well monomers with a single short-ranged interaction site has been recently observed using the Gibbs ensemble Monte Carlo technique 关L. A. Davies, G. Jackson, and L. F. Rull, Phys. Rev. Lett. 82, 5285 共1999兲兴. This model system has unfavorable mean-field interactions between unlike components which leads to phase separation at intermediate temperatures; the addition of a directional bonding site leads to association and miscibilty of the system at low temperatures. In this work we present a detailed study of the effect of a variation in pressure and of the strength of the bonding interaction on the phase equilibria of such a model system by Gibbs ensemble simulation. The phase diagram is dominated by regions of liquid-liquid immiscibility which are bounded at high temperatures by an upper critical solution temperature and by a lower critical solution temperature 共LCST兲for specific values of the pressure and association strength. This closed-loop region is seen to increase in size as the pressure of the system is increased. For weak bonding interaction strengths the system does not possess a LCST and is seen to exhibit regions of two-phase vapor-liquid coexistence which are separated from the region of liquid-liquid immiscibility by a three-phase line. The phase equilibria of the same model system is also determined using the statistical associating fluid theory as adapted for potentials of variable range; the theory provides a good description of the closed-loop immiscibility and other features of the phase diagram. PACS number共s兲: 64.60.⫺i, 61.20.Qg, 05.20.Jj I. INTRODUCTION The existence of closed-loop regions of liquid-liquid immiscibility in the phase diagrams of aqueous fluid mixtures such as aliphatic and aromatic alcohols 关1兴, amines 关2兴, and nonionic surfactants 关3兴can be readily understood in terms of the nature of the intermolecular interactions. The association of unlike species via directional interactions such as hydrogen bonds leads to the low-temperature miscibility of the system. These bonds break as the temperature of the system increases resulting in phase separation above the lower critical solution temperature 共LCST兲due to the residual weak unlike interactions. The increase in the kinetic energy, and hence the entropy of the molecules at higher temperatures ensures that the system is miscible above the upper critical solution temperature 共UCST兲. A detailed description of this type of reentrant phase behavior is given in Ref. 关4兴. Early theoretical studies of this particular phenomenon of fluid phase equilibria used lattice models with orientationally dependent intermolecular interactions 共see Ref. 关5兴for a review兲. Such approaches are based on the assumption that a lattice gives an adequate description of the structure of a liquid, which is a significant oversimplification. A description of closed-loop immiscibility obtained with continuum fluid theories is more realistic, since features such as pressure and density can be directly related to those of a real fluid. The first continuum study of a system with a well-defined directional intermolecular potential which exhibits closedloop immiscibility 关5兴used the Wertheim approach 关6–9兴to describe the association interactions within the framework of the statistical associating fluid theory 共SAFT兲关10–13兴. More recently, the link between the existence of the LCST and short-ranged, directional association interactions has been established for a continuum model using computer simulation 关14兴. Here, we illustrate how continuum fluid theories, more specifically the statistical associating fluid theory as adapted to potentials of variable range 共SAFT-VR兲关15,16兴, can be used to give complementary results in a study of the effect of pressure and association interaction strength on the phase diagram of a simple model system which exhibits reentrant miscibility. It is well known that continuum fluid theories such as SAFT can be applied to obtain an accurate prediction of the phase behavior of real systems, providing suitable models are used. The seminal study of Scott and van Konynenburg 关17,18兴used the van der Waals equation of state to describe, and hence classify, the types of phase behavior exhibited by binary fluid mixtures. Type VI behavior, which incorporates regions of closed-loop immiscibility, is the only type which cannot be predicted theoretically with the van der Waals equation of state. The SAFT methodology consists of a perturbation theory about a hard-sphere reference system, where PHYSICAL REVIEW E MARCH 2000VOLUME 61, NUMBER 3 PRE 61 1063-651X/2000/61共3兲/2245共12兲/$15.00 2245 ©2000 The American Physical Society the dispersive interactions are given by a simple hightemperature expansion, and a term is included to account for the association of species. The extent of the closed-loop region in the model system studied in Ref. 关5兴with a simplified 共mean-field兲SAFT equation of state depends on both the pressure of the system and the strength of the short-ranged interaction between species in the mixture. A global study of the phase behavior of aqueous mixtures of alcohols using the SAFT approach predicts the existence of closed-loop regions of immiscibility for model associating systems of intermediate chain lengths with specific values of the site-site interaction energy 关19兴. Similar models were used to give an accurate prediction of the phase equilibria of aqueous systems of 2-butanol, butoxyethanol and long-chain nonionic alkyl polyoxyethylene surfactant molecules (CiEj) with the SAFT approach 关20,21兴. This simple approach is shown to give an adequate description of the phase behavior of systems where association is the dominant contribution to the Helmholtz free energy. A recent extension of the SAFT methodology employs the Barker and Henderson high-temperature perturbation theory 关22–24兴to give a more accurate description of the dispersive interactions, and can also be used to predict the phase behavior of systems which interact via potentials of variable range. This so-called SAFT-VR approach 关15,16兴is described in more detail in Sec. III. A major advantage of the SAFT approach is that only the monomer Helmholtz free energy and the contact value of the monomer cavity function are required to describe the equation of state of chain molecules. The success of this approach lies in its ability to provide an accurate prediction of the thermodynamics and hence the phase behavior of systems where molecular shape and associating phenomena are dominant features. Additionally, the SAFT approach is very versatile; it can be used to describe the phase behavior of a number of different systems 共see Ref. 关15兴for a brief review兲. The molecular based nature of the SAFT-VR equation allows the results to be compared with computer simulation. Here we obtain the phase behavior of a symmetrical model mixture using both the Gibbs ensemble Monte Carlo simulation method and the SAFT-VR equation of state. The system consists of a binary mixture of equal sized spheres with unfavorable mean-field interactions 共which lead to liquid-liquid phase separation at intermediate temperatures兲 and a single association site which provides the directional interactions between unlike species 共which leads to miscibility at low temperatures兲. The particles interact with a squarewell potential: uij共rij兲⫽ 再 ⫹⬁if rij⬍ ␴ ij, ⫺␧ij if ␴ ij⭐rij⬍␭ij ␴ ij, 0ifrij⭓␭ij ␴ ij, 共1兲 where rij is the distance between two particles. The contact distance is ␴ ij and the parameters ␭ij and ␧ij are the range and depth of the potential well for the i-jinteraction, respectively. For our particular system ␴ ⫽ ␴ 11⫽ ␴ 12⫽ ␴ 22 ,␧ ⫽␧11⫽␧22 , and ␧12⫽0, with ␭⫽␭11⫽␭12⫽␭22⫽1.5. A single square-well association site of fixed depth ␧a,band width ␭a,b⫽0.55 is placed at a distance rd⫽0.25 ␴ from the center of each sphere. Bonding is only permitted between unlike species, i.e., between site aon component 1 and site b on component 2. This particular model system is known from previous studies to to exhibit regions of closed-loop immiscibility 关5,14兴. A series of simulations are performed for different values of the pressure and site-site interaction energies in order to illustrate their effect on the phase diagram. The results obtained are compared with those obtained using the SAFT-VR approach for corresponding conditions. II. GIBBS ENSEMBLE MONTE CARLO SIMULATIONS The Gibbs ensemble Monte Carlo 共GEMC兲method 关25,26兴is the most common direct simulation technique used for the determination of phase equilibria in fluid systems, particularly in mixtures 关27,28兴. The isothermal-isobaric (NPT) version of the technique consists of a separate simulation in two regions aand bwhich are in thermodynamic equilibrium but not in physical contact. These regions have volumes Vaand Vband contain Naand Nbparticles, such that V⫽Va⫹Vband N⫽Na⫹Nb. The two regions are in thermodynamic equilibrium but not in physical contact. Three different Monte Carlo moves are performed in order to satisfy the conditions of thermodynamic equilibrium: particle displacements and reorientations within either subsystem, to maintain equality of temperature; volume changes of either subsystem, to maintain equality of pressure; and particle interchanges between the two subsystems, to maintain equality of chemical potential. The energy Ej, volume Vj, and composition xi j⫽Ni j/Njof particles of type iin subsystem jthus vary during the course of the simulation. The acceptance criterion for each of these moves in the NPT version of the Gibbs ensemble is governed by the pseudo-Boltzmann probability distribution, see Refs. 关26,29兴for details: PGibbs⫽exp 冋 ln 冉 N1! N1 a!N1 b! 冊 ⫹ln 冉 N2! N2 a!N2 b! 冊 ⫹NalnVa ⫹NblnVb⫺PVa kT ⫺PVb kT ⫺Ea共Na兲 kT ⫺Eb共Nb兲 kT 册 . 共2兲 We examine a symmetrical binary mixture of N⫽N1⫹N2 particles at a constant temperature Tand a constant pressure P.Simulations are performed in cubic boxes with the particles initially arranged on a face-centered-cubic 共fcc兲lattice. The usual periodic boundary conditions and minimum image convention are used 关30兴. Initial guesses for the coexisting densities and compositions at each pressure and temperature are made by using the corresponding SAFT-VR solutions 共see following section兲; it is important to ensure that the overall composition of the system lies somewhere between the compositions of the two coexisting phases. The chemical potential is determined with the Widom test particle technique 关31兴as adapted to the GEMC approach 关32兴, in order to ensure that phase equilibria is achieved. One simulation cycle consists of Ndisplacements and reorientations in each box, one volume change for either box, and a specific number of particle interchanges. The maximum displacement and volume change are adjusted to give an acceptance ratio of 2246 PRE 61LOWRI A. DAVIES, GEORGE JACKSON, AND LUIS F. RULL between 30 and 40%, and the number of insertions is controlled so that between 1 and 3% of particles are interchanged each cycle. Simulations are performed with systems of N⫽512, N⫽1000, and N⫽1728 particles at a series of coexistence pressures for different values of the site-site interaction energy. An initial simulation of 50000 cycles is performed to equilibrate the subsystems, before averaging for between 100000 and 250000 cycles. III. SAFT-VR EQUATION OF STATE FOR SQUAREWELL MIXTURES The SAFT-VR equation of state for a mixture of associating chain molecules is written in terms of four separate contributions to the Helmholtz free energy 关15,16兴 A NkT⫽AIDEAL NkT ⫹AMONO NkT ⫹ACHAIN NkT ⫹AASSOC NkT ,共3兲 where Nis the number of chain molecules in the mixture, kis Boltzmann’s constant, and Tis the temperature. In this equation AIDEAL is the ideal free energy, AMONO is the residual free energy due to the monomer segments, ACHAIN is the residual contribution due to the formation of chains of monomers, and AASSOC is the residual contribution to the free energy due to intermolecular association via sites placed on the monomer segments. We present the general expressions for each contribution to the Helmholtz free energy in the above equation together with those which are specific for the associating system examined in this work. The free energy of an ideal n-component mixture is given by 关33兴 AIDEAL NkT ⫽兺 i⫽1 n xiln ␳ i⌳i⫺1⫽x1ln ␳ 1⌳1 3⫹x2ln ␳ 2⌳2 3⫺1, 共4兲 where xi⫽Ni/Nis the mole fraction, ␳ i⫽Ni/Vis the number density, and ⌳iis the thermal de Broglie wavelength of species i. The monomer free energy is AMONO NkT ⫽ 冉 兺 i⫽1 n ximi 冊 AM NskT⫽ 冉 兺 i⫽1 n ximi 冊 aM⫽aM,共5兲 where miis the number of spherical segments in each chain i, so that m⫽1 for monomers, and Nsis the total number of segments. The monomer free energy per segment of the mixture aM⫽AM/(NskT) is obtained from the Barker and Henderson high-temperature expansion 关22–24兴: aM⫽aHS⫹ ␤ a1⫹ ␤ 2a2,共6兲 where aHS is the free energy for a mixture of hard spheres, ␤ ⫽1/kT,a1and a2are the first two perturbation terms associated with the attractive energy. The free energy of the reference hard-sphere mixture is obtained from the expression of Boublı ´k关34兴and Mansoori et al. 关35兴: aHS⫽6 ␲␳ s 冋 冉 ␨ 2 3 ␨ 3 2⫺ ␨ 0 冊 ln共1⫺ ␨ 3兲⫹3 ␨ 1 ␨ 2 1⫺ ␨ 3 ⫹ ␨ 2 3 ␨ 3共1⫺ ␨ 3兲2 册 . 共7兲 In this expression ␳ s⫽Ns/Vis the number density of the mixture in terms of the number of spherical segments. Note that in general ␳ s⫽ ␳ (兺iximi), where ␳ is the total number density of the mixture. In this case, since we consider only monatomic species ␳ s⫽ ␳ . The reduced densities ␨ lare defined as ␨ l⫽ ␲ 6 ␳ s 冋 兺 i⫽1 n xs,i共 ␴ i兲l 册 ,共8兲 where ␴ iis the diameter of spherical segments of chain i, and xs,iis the mole fraction of segments of type iin the mixture, which is given by xs,i⫽mixi 兺 k⫽1 n mkxk .共9兲 Since xs,i⫽xiand ␴ ⫽ ␴ 1⫽ ␴ 2in the system studied here, the free energy of the reference hard-sphere mixture reduces to the Carnahan and Starling expression 关36,33兴 aHS⫽4 ␩ ⫺3 ␩ 2 共1⫺ ␩ 兲2,共10兲 where ␩ ⫽ ␲␳␴ 3/6 is the packing fraction of the pure component, which is equivalent to ␨ 3in Eq. 共8兲. The mean-attractive energy a1in the perturbation expansion is given by a1⫽兺 i⫽1 n 兺 j⫽1 n xs,ixs,ja1 ij,共11兲 where a1 ij⫽⫺2 ␲␳ s␧ij 冕 ␴ ij ⬁rij 2gij HS共rij; ␨ 3兲drij,共12兲 and gij HS is the radial distribution function for a mixture of hard spheres. The integral is transformed by applying the mean-value theorem 关15兴giving an expression for a1in terms of the contact value of gij HS : a1⫽⫺ ␳ s兺 i⫽1 n 兺 j⫽1 n xs,ixs,j ␣ ij VDWgij HS关 ␴ ij; ␨ 3 eff兴,共13兲 where ␣ ij VDW⫽2 ␲ ␧ij ␴ ij 3共␭ij 3⫺1兲/3 共14兲 is the van der Waals attractive constant for the square-well interaction i-j. The contact value of the radial distribution function for the hard-sphere reference system, gij HS关 ␴ ij; ␨ 3 eff兴 is evaluated at an effective packing fraction ␨ eff, using the expression of Boublı ´k关34兴 PRE 61 2247CLOSED-LOOP PHASE EQUILIBRIA OF A . . . gij HS关 ␴ ij; ␨ 3 eff兴⫽1 1⫺ ␨ 3 eff⫹3 ␴ ii ␴ jj ␴ ii⫹ ␴ jj ␨ 2 eff 共1⫺ ␨ 3 eff兲2 ⫹2 冉 ␴ ii ␴ jj ␴ ii⫹ ␴ jj 冊 2 ␨ 2 eff2 共1⫺ ␨ 3 eff兲3.共15兲 For the binary mixture studied here ␴ ⫽ ␴ 11⫽ ␴ 22 ,␧⫽␧11 ⫽␧22 ,␭⫽␭11⫽␭12⫽␭22 , and ␧12⫽0, so that the meanattractive energy of Eq. 共11兲reduces to a1⫽x1 2a1 11⫹x2 2a1 22⫽共x1 2⫹x2 2兲a1 SW,共16兲 where a1 SW⫽⫺ ␳ s ␣ VDWgHS关 ␴ ; ␩ eff兴,共17兲 with ␣ VDW⫽2 ␲ ␧ ␴ 3共␭3⫺1兲/3. 共18兲 The Carnahan and Starling equation for the contact value of the pair correlation function is used 关36,33兴, gHS关 ␴ ; ␩ eff兴⫽1⫺ ␩ eff/2 共1⫺ ␩ eff兲3,共19兲 with ␩ eff⫽ ␨ 3 eff . The parameterization for ␩ eff obtained for the pure square-well fluid 关15兴is used, where ␩ eff共 ␩ ,␭兲⫽c1共␭兲 ␩ ⫹c2共␭兲 ␩ 2⫹c3共␭兲 ␩ 3,共20兲 and the coefficients c1,c2, and c3are given by 冉 c1 c2 c3 冊 ⫽ 冉 2.25855 ⫺1.50349 0.249434 ⫺0.669270 1.40049 ⫺0.827739 10.1576 ⫺15.0427 5.30827 冊 ⫻ 冉 1 ␭ ␭2 冊 .共21兲 This corresponds to the MX1 or MX3 mixing rules of Ref. 关16兴. The fluctuation term of the free energy is given by a2⫽兺 i⫽1 n 兺 j⫽1 n xs,ixs,ja2 ij,共22兲 and each of the terms a2 ij are obtained with the local compressibility approximation 共LCA兲关22,23兴, a2 ij⫽1 2KHS␧ij ␳ s ⳵ a1 ij ⳵␳ s,共23兲 where KHS is the isothermal compressibility for a mixture of hard spheres which is given by the Percus-Yevick expression 关37兴 KHS⫽ ␨ 0共1⫺ ␨ 3兲4 ␨ 0共1⫺ ␨ 3兲2⫹6 ␨ 1 ␨ 2共1⫺ ␨ 3兲⫹9 ␨ 2 3.共24兲 For our particular system a2⫽x1 2a2 11⫹x2 2a2 22⫽共x1 2⫹x2 2兲a2 SW,共25兲 with a2 SW⫽1 2KHS␧ ␳ s ⳵ a1 SW ⳵␳ s,共26兲 where KHS is now the pure component expression KHS⫽共1⫺ ␩ 兲4 共1⫺ ␩ 兲2⫹6 ␩ 共1⫺ ␩ 兲⫹9 ␩ 2.共27兲 Since the system examined here consists only of monatomic species the contribution to the free energy in Eq. 共3兲due to chain formation is zero (ACHAIN⫽0). Within the framework of the theory of Wertheim the contribution to the free energy in the SAFT-VR equation of state due to the association mediated by sisites on molecules of species iis described by 关11兴 AASSOC NkT ⫽兺 i⫽1 n xi 冋 兺 a⫽1 si 冉 lnXa,i⫺Xa,i 2 冊 ⫹si 2 册 .共28兲 The first sum is over the number of species iin the mixture and the second sum is over all sites aon a molecule of type i. The fractions of molecules of species inot bonded at a particular site a,Xa,i, is given by solution of the simultaneous equations: Xa,i⫽1 1⫹兺 j⫽1 n 兺 b⫽1 sj ␳ xjXb,j⌬a,b,i,j ,共29兲 where ⌬a,b,i,j⫽Ka,b,i,jfa,b,i,jgij SW共 ␴ ij兲共30兲 is specific for each a-bsite-site interaction, and incorporates the volume available for bonding Ka,b,i,j, and the strength of the association via the Mayer function fa,b,i,j⫽exp(␧ab /kT) ⫺1 of the square-well potential. The contact value of the radial distribution function for the square-well interaction gij SW( ␴ ij) is evaluated within the SAFT-VR approach for mixtures 关15,16兴using a first order perturbation expansion gij SW共 ␴ ij兲⫽gij HS共 ␴ ij兲⫹ ␤⑀ ijg1共 ␴ ij兲.共31兲 The contact value of the radial distribution function for the hard-sphere reference system gij HS关 ␴ ij; ␨ 3兴is evaluated with the actual packing fraction ␨ 3, using the expression of Eq. 共15兲. The first order perturbation term g1( ␴ ij) is obtained from a self-consistent calculation of the pressure using the Clausius virial theorem 关15兴: 2248 PRE 61LOWRI A. DAVIES, GEORGE JACKSON, AND LUIS F. RULL gij SW关 ␴ ij; ␨ 3兴⫽gij HS关 ␴ ij; ␨ 3兴⫹ ␤ ␧ij 冋 gij HS关 ␴ ij; ␨ 3 eff兴⫹共␭ij 3⫺1兲 ⫻ 冉 ␭ij 3 ⳵ gij HS关 ␴ ij; ␨ 3 eff兴 ⳵ ␭ij ⫺ ␳ s ⳵ gij HS关 ␴ ij; ␨ 3 eff兴 ⳵␳ s 冊 册 . 共32兲 For our system, the contact value of the square-well distribution function of Eq. 共32兲simplifies to gSW关 ␴ ; ␩ 兴⫽gHS关 ␴ ; ␩ 兴⫹ ␤ ␧ 冋 gHS关 ␴ ; ␩ eff兴⫹共␭3⫺1兲 ⫻ ⳵ gHS关 ␴ ; ␩ eff兴 ⳵ ␩ eff 冉 ␭ 3 ⳵ ␩ eff ⳵ ␭⫺ ␩ ⳵ ␩ eff ⳵ ␩ 冊 册 ,共33兲 where the Carnahan and Starling expressions 关36,33兴for gHS关 ␴ ; ␩ 兴and gHS关 ␴ ; ␩ eff兴are used with the total and effective packing fractions, respectively. The parametrization for ␩ eff for the pure square-well fluid is given in Eqs. 共20兲and 共21兲. The components in the mixture examined here have only have one interaction site per sphere and bonding is only permitted between unlike components in the mixture, i.e., between site aon component 1 and site bon component 2 when the two sites are within a distance rc. The expression for the contribution to the Helmholtz free energy due to association, Eq. 共28兲, thus simplifies to AASSOC NkT ⫽x1 冉 lnXa⫺Xa 2⫹1 2 冊 ⫹x2 冉 lnXb⫺Xb 2⫹1 2 冊 . 共34兲 The fractions Xaand Xbof molecules of species 1 and 2 not bonded at sites aand b, respectively, are given by Xa⫽1 1⫹ ␳ x2Xb⌬a,b共35兲 and Xb⫽1 1⫹ ␳ x1Xa⌬a,b,共36兲 where ⌬a,b⫽Ka,bfa,bgSW共 ␴ 兲.共37兲 The volume available for bonding between two sites aand b which are positioned at a distance rdfrom the center of spheres with diameter ␴ and have an interaction range rcis given by 关10兴 Ka,b⫽4 ␲ 2 ␴ 兵 ln关共rc⫹2rd兲/ ␴ 兴共6rc 3⫹18rc 2rd⫺24rd 3兲 ⫹共rc⫹2rd⫺ ␴ 兲共22rd 2⫺5rcrd⫺7rd ␴ ⫺8rc 2⫹rc ␴ ⫹ ␴ 2兲 其 /共72rd 2兲.共38兲 In order to determine the conditions of phase equilibria for the model system, the pressure and chemical potential must be obtained from the expression for the Helmholtz free energy. The chemical potential ␮ iof species iis given by ␮ i kT⫽ 冉 ⳵ A/kT ⳵ Ni 冊 T,V,Nj⫽i ,共39兲 where Niis the number of molecules of species i. The pressure Pmay be calculated through the compressibility factor Zas Z⫽PV NkT⫽兺 i⫽1 n xi ␮ i kT⫹A NkT.共40兲 These functions are used in the numerical determination of the phase behavior of the mixture, using a simplex method 关38兴. IV. CRITICAL BEHAVIOR Generally, finite size effects in the Gibbs ensemble are considered to be small for regions away from the critical point 关39兴, so that simulation data from these regions can be used as experimental data to estimate the critical temperature Tc(N) via extrapolation of a Wegner expansion 关40兴which includes scaling corrections 关41兴. This method has been used successfully to predict the critical parameters of the pure Lennard-Jones fluid 关42兴, the pure square-well fluid 关43兴, and selected square-well mixtures 关44兴. Before such a fitting procedure can be performed it is necessary to identify the universality class 关45兴of the phase transitions occurring at the UCST and the LCST for the symmetrical associating mixture. As the critical point between two phases aand bis approached the critical exponent ␤ is defined in terms of the critical temperature of the finite system as 关46兴 ⌬x共N兲⫽ 兩 xa⫺xb 兩 ⫽B0共N兲 冏 1⫺T Tc共N兲 冏 ␤ ,共41兲 where B0(N) is the leading amplitude term. We assume that the phase transitions of the binary mixture at the LCST and the UCST both belong to the Ising universality class, so that the critical exponent has the universal value of ␤ ⫽0.325. The coexistence curve thus has a cubic form in the critical region. Providing that the correlation length ␨ aremains less than the linear box length L, the value of ␤ ⫽0.325 can be assumed to be valid for the duration of the simulations. A change in the value of the critical exponent from ␤ ⫽0.325 共universal兲to ␤ ⫽0.5 共mean-field兲only occurs when the correlation length becomes greater than L关47兴. A value of ␤ ⫽0.5 results in a quadratic coexistence curve close to the critical temperature. This ‘‘crossover’’ of critical exponents has only been observed in Gibbs ensemble simulation studies of a two-dimensional lattice gas 关47兴and in a two dimensional square-well mixture 关48兴. Simulation studies of the square-well fluid in three dimensions 关48,44,49兴, show no sign of this crossover behavior. This is in agreement with the opinion that this type of crossover cannot be observed in three dimensional systems using regular simulation techniques 关50兴. Away from the critical region the power law of Eq. 共41兲 fails to accurately describe the shape of the coexistence PRE 61 2249CLOSED-LOOP PHASE EQUILIBRIA OF A . . . curve. Corrections to the renormaliZation-group 共RG兲-based scaling laws must be included in order to obtain more accurate results 关41兴. Such corrections can be written as an expansion in t⫽1⫺T/Tc(N), in a so-called Wegner 关40兴ex- pansion 兩 xa⫺xb 兩 ⫽B0 兩 t 兩 ␤ ⫹B1 兩 t 兩 ␤ ⫹⌬1⫹B2 兩 t 兩 ␤ ⫹2⌬1⫹•••,共42兲 where ⌬1⫽0.5 is a RG gap exponent, and the Biterms are the correction amplitudes. This expansion leads an estimate of the critical temperature while assuming a universal value for the critical exponent. A similar expansion can be written for the diameter of the coexistence curve 关41兴 共xa⫹xb兲 2⫽xc⫹C1 兩 t 兩 ␺ ⫹C2 兩 t 兩 ⫹C3 兩 t 兩 ␺ ⫹⌬1⫹•••,共43兲 where xcis the critical composition, the Citerms are the coefficients of the expansion, and ␺ is an exponent which characterizes the anomaly in the diameter of the curve. These two expansions can be combined to give one expression for the coexisting compositions x⫾⫽xc⫹C1 兩 t 兩 ␺ ⫹C2 兩 t 兩 ⫹C3 兩 t 兩 ␺ ⫹⌬1⫹••• ⫾1 2共B0 兩 t 兩 ␤ ⫹B1 兩 t 兩 ␤ ⫹⌬1⫹B2 兩 t 兩 ␤ ⫹2⌬1⫹•••兲,共44兲 where x⫺and x⫹represent the smaller and larger coexistence compositions, respectively. Similar expressions can be written in terms of the coexistence densities. It is important to note that the Wegner extended scaling of the GEMC simulation data for simple model systems such as the square-well fluid provides estimates for the critical parameters with an equivalent accuracy to the more thorough analysis of the critical region obtained by using finite-size scaling approaches 共see Refs. 关44,48,51兴兲, despite being significantly less computationally demanding. For the associating square-well mixture studied here this expression can be simplified due to the symmetry of the system. The critical composition is known, xc⫽0.5, and the Cicoefficients can be disregarded, since the diameter is symmetrical about xc⫽0.5. Additionally, inclusion of only the first correction to scaling term B1is known to give an adequate description of the liquid-liquid coexistence region 关44兴. The expression for the compositions of the coexisting phases becomes x⫾⫽1 2⫾1 2共B0 兩 t 兩 ␤ ⫹B1 兩 t 兩 ␤ ⫹⌬1兲.共45兲 Assuming values of ␤ ⫽0.325 and ⌬1⫽0.5 the above equation can be fitted to the simulation data using a standard least-squares procedure 关38兴to give estimates of the critical temperature. For coexistence curves which possess both an upper and a lower critical solution temperature it is possible to use two such power laws, one to estimate the UCST, and one to estimate the LCST. These two expressions can be combined to give a complete description of the closed-loop coexistence curve. A hyperbolic tangent tanh共y兲⫽exp共y兲⫺exp共⫺y兲 exp共y兲⫹exp共⫺y兲共46兲 is used as a switching function between the expansion for the UCST and that for the LCST. Hence, the expansion fitted to data in the region of the UCST is used above a certain temperature Tswitch , and the expansion fitted to data in the region of the LCST is used at temperatures below this point. The compositions of the coexisting liquid phases in the closedloop region are obtained as x⫾ loop⫽F共T兲x⫾ UCST⫹„1⫺F共T兲…x⫾ LCST ,共47兲 where F(T) is given by F共T兲⫽1 2⫹1 2tanh共y兲,共48兲 with y⫽T⫺Tswitch 2␭switch .共49兲 The parameter ␭switch controls the steepness of the switching function. V. RESULTS The reduced thermodynamic variables, temperature T*⫽kT/␧, pressure P*⫽P ␴ 3/␧, and site-site interaction energy ␧a,b *⫽␧a,b/␧are used in the following discussion. It is also convenient to reduce the pressure and temperature with respect to the critical point of one of the components FIG. 1. Temperature-composition slice of the coexistence region for the symmetrical square-well mixture at a reduced pressure of P*⫽P ␴ 3/␧⫽0.756 with a bonding interaction of ␧a,b *⫽␧a,b/␧ ⫽13. The reduced temperature is defined as T*⫽kT/␧. The triangles correspond to the GEMC data for a system of N⫽512 particles, the squares to a system of N⫽1000 particles, and the circles to a system of N⫽1728 particles. The error bars correspond to one standard deviation. The dashed, continuous, and dotted curves correspond to the SAFT-VR prediction for the same model with bonding strengths of ␧a,b *⫽12, 13, and 14, respectively. The predicted three-phase line is represented by the horizontal line. 2250 PRE 61 LOWRI A. DAVIES, GEORGE JACKSON, AND LUIS F. RULL such that Tr⫽T*/Tc *and Pr⫽P*/Pc *, where the subscript c denotes the critical value of a variable. We use the critical values of Tc *⫽1.219⫾0.008 and Pc *⫽0.108⫾0.016 for the pure-component square-well system with range ␭⫽1.5 obtained from GEMC simulation results 关43兴. Constant-pressure temperature-composition Tx slices of the coexistence regions are shown in Figs. 1–4. The symmetry of the phase diagrams about x2⫽0.5 is a result of the symmetry of the model system. The distinction between gas and liquid fluid phases is arbitrary at high pressures and temperatures; here we use the term ‘‘gas’’ to describe a phase which has its origins in the vapor phase of the pure component and which thus has the lower density. The pressure of P*⫽0.756 of Fig. 1 is approximately seven times the critical pressure of the pure square-well system with ␭⫽1.5. The Gibbs ensemble simulation data with ␧a,b *⫽13 for three system sizes are shown; there is a negligible system size effect. The region of liquid-liquid immiscibility bounded by an UCST at high temperatures is clearly shown, together with low temperature regions of gas-liquid immiscibility below a three-phase line. The curves represent the SAFT-VR prediction for the equivalent model, with three different strengths of the site-site interaction energy. For a value of ␧a,b *⫽12 the theory predicts immiscibility of the system at all temperatures below an USCT. As the site-site interaction is increased to ␧a,b *⫽13 the prediction includes regions of gas-liquid immiscibility at lower temperatures, and the three-phase line can be observed at a temperature close to the critical temperature of the pure component. Increasing the strength of the bonding interaction further to ␧a,b *⫽14, leads to a complete change in the nature of the Tx FIG. 2. Temperature-composition slice of the coexistence region for the symmetrical square-well mixture at a reduced pressure of P*⫽P ␴ 3/␧⫽1.08 with a bonding interaction of ␧a,b *⫽␧a,b/␧ ⫽13. The reduced temperature is defined as T*⫽kT/␧. The triangles correspond to the GEMC data for a system of N⫽512 particles, the squares to a system of N⫽1000 particles, and the circles to a system of N⫽1728 particles. The error bars correspond to one standard deviation. The continuous, dotted, and dashed curves correspond to the SAFT-VR prediction for the same model with bonding strengths of ␧a,b *⫽13, 14, and 15, respectively. The predicted three-phase lines are represented by the horizontal lines. FIG. 3. Temperature-composition slice of the coexistence region for the symmetrical square-well mixture at a reduced pressure of P*⫽P ␴ 3/␧⫽1.28 with a bonding interaction of ␧a,b *⫽␧a,b/␧ ⫽14.5. The reduced temperature is defined as T*⫽kT/␧. The squares correspond to the GEMC data for a system of N⫽1000 particles, and the circles to a system of N⫽1728 particles. The error bars correspond to one standard deviation. The continuous and dashed curves correspond to the SAFT-VR prediction for the same model with bonding strengths of ␧a,b *⫽14.5 and 15, respectively. The dotted curve corresponds to the fit of the GEMC simulation results obtained using the Wegner expansion of Eq. 共47兲, which includes a first correction to scaling. The predicted three-phase line is represented by the horizontal line. FIG. 4. Temperature-composition slice of the coexistence region for the symmetrical square-well mixture at a reduced pressure of P*⫽P ␴ 3/␧⫽1.48 with a bonding interaction of ␧a,b *⫽␧a,b/␧ ⫽14.5. The reduced temperature is defined as T*⫽kT/␧. The circles correspond to the GEMC data for a system of N⫽1728 particles. The error bars correspond to one standard deviation. The continuous and dashed curves correspond to the SAFT-VR prediction for the same model with bonding strengths of ␧a,b *⫽14.5 and 15, respectively. The dotted curve corresponds to the fit of the GEMC simulation results obtained using the Wegner expansion of Eq. 共47兲, which includes a first correction to scaling. The predicted three-phase lines are represented by the horizontal lines. PRE 61 2251CLOSED-LOOP PHASE EQUILIBRIA OF A . . . slice, with miscibility above the two low-temperature regions of gas-liquid coexistence. Thus the SAFT-VR equation predicts that an increase in the site-site interaction energy from ␧a,b *⫽13 to ␧a,b *⫽14 is sufficient to bring about the disappearance of the liquid-liquid immiscibility at this pressure. The low-temperature regions of gas-liquid immiscibility obtained from simulation appear at slightly higher temperatures than those predicted by the SAFT-VR equation of state for ␧a,b *⫽13. The equivalent temperature composition Tx slice for a higher pressure of P*⫽1.08 共approximately ten times the critical pressure of the pure component兲is shown in Fig. 2. The simulation results again present a region of liquid-liquid immiscibility bounded at high temperatures by an UCST, and low-temperature regions of gas-liquid immiscibility below the three-phase line. The region of liquid-liquid immiscibility is seen to be larger than in the lower pressure case of Fig. 1, but a closed-loop of immiscibility is still not observed for ␧a,b *⫽13. For this value of the site-site interaction energy the SAFT-VR equation of state predicts a region of liquidliquid immiscibility bounded by an UCST and small regions of low temperature gas-liquid coexistence below the three phase line. For a site-site interaction energy of ␧a,b *⫽14, the low temperature gas-liquid regions predicted by the theory are larger than in the ␧a,b *⫽13 case; the onset of a closedloop region can be observed. Upon increasing the association energy to ␧a,b *⫽15 the SAFT-VR equation of state predicts the existence of a small closed-loop of immiscibility bounded at high temperatures by an UCST at T*⫽1.98 and at low temperatures by a LCST at T*⫽1.84, together with low-temperature regions of gas-liquid coexistence. At this stage we can conclude that increasing the pressure from P*⫽0.756 to P*⫽1.08 for a system with an interaction energy of ␧a,b *⫽13 does not bring about a sufficient change in the phase diagram of the system for a closed-loop of immiscibility to be observed by GEMC simulation. An increase in the strength of the bonding interaction, together with a slight increase in pressure is required in order to observe a closed-loop of immiscibility 关14兴. The temperature-composition slice for a pressure of P*⫽1.28 is shown in Fig. 3. Here, Gibbs ensemble simulations are performed for a mixture with a site-site interaction energy of ␧a,b *⫽14.5; the data is reported in Tables I and II. The region of closed-loop immiscibility can be clearly observed, bounded by both an UCST and a LCST. No lowtemperature regions of gas-liquid coexistence were isolated for this system. The SAFT-VR predictions at a pressure of P*⫽1.28 with site-site interactions of ␧a,b *⫽14.5 and ␧a,b * ⫽15 are also shown. For ␧a,b *⫽14.5 the theory predicts a region of liquid-liquid immiscibility bounded by an UCST and regions of low-temperature gas-liquid coexistence at temperatures below the three-phase line. Increasing the strength of the site-site interaction to ␧a,b *⫽15 leads to the prediction of closed-loop immiscibility, with upper and lower boundaries which are similar to those obtained using Gibbs ensemble simulation. It is encouraging to observe the TABLE I. Liquid-liquid coexistence data obtained from NPT Gibbs-ensemble Monte Carlo simulations for the symmetrical mixture of square wells with a range ␭⫽1.5, and a single interaction site of reduced depth ␧a,b *⫽␧a,b/␧⫽14.5. The fixed variables during the simulation are the number of particles N, the reduced pressure P*⫽P ␴ 3/␧⫽1.28 and the reduced temperature T*⫽kT/␧. The packing fractions ␩ and mole fractions x2in the coexisting liquid phases are labeled l1and l2, respectively; the uncertainties correspond to one standard deviation. NT* ␩ l1 ␩ l2x2,l1x2,l2 1728 1.58 0.289⫾0.004 0.289⫾0.004 0.810⫾0.021 0.195⫾0.020 1728 1.59 0.289⫾0.004 0.290⫾0.003 0.814⫾0.033 0.185⫾0.035 1728 1.60 0.291⫾0.003 0.292⫾0.004 0.835⫾0.026 0.159⫾0.028 1728 1.61 0.292⫾0.003 0.292⫾0.003 0.860⫾0.017 0.136⫾0.017 1728 1.62 0.289⫾0.004 0.290⫾0.004 0.862⫾0.020 0.140⫾0.023 1728 1.63 0.290⫾0.005 0.290⫾0.005 0.864⫾0.022 0.135⫾0.027 1000 1.64 0.284⫾0.006 0.279⫾0.006 0.860⫾0.032 0.177⫾0.031 1000 1.65 0.292⫾0.007 0.285⫾0.007 0.886⫾0.026 0.176⫾0.036 1000 1.68 0.282⫾0.007 0.277⫾0.006 0.864⫾0.027 0.166⫾0.025 1000 1.70 0.275⫾0.007 0.273⫾0.007 0.859⫾0.032 0.152⫾0.031 1000 1.72 0.272⫾0.007 0.274⫾0.010 0.852⫾0.032 0.131⫾0.041 1000 1.74 0.271⫾0.005 0.267⫾0.005 0.878⫾0.024 0.150⫾0.024 1000 1.75 0.266⫾0.004 0.266⫾0.004 0.849⫾0.027 0.154⫾0.027 1000 1.76 0.264⫾0.006 0.261⫾0.005 0.850⫾0.033 0.170⫾0.029 1000 1.78 0.262⫾0.007 0.257⫾0.006 0.854⫾0.032 0.193⫾0.030 1000 1.80 0.262⫾0.005 0.260⫾0.004 0.860⫾0.028 0.154⫾0.029 1728 1.81 0.254⫾0.004 0.254⫾0.004 0.827⫾0.029 0.177⫾0.029 1728 1.82 0.256⫾0.003 0.255⫾0.003 0.833⫾0.028 0.178⫾0.030 1728 1.83 0.254⫾0.003 0.255⫾0.003 0.816⫾0.021 0.168⫾0.021 1728 1.84 0.254⫾0.003 0.254⫾0.003 0.824⫾0.020 0.169⫾0.021 1728 1.85 0.251⫾0.004 0.249⫾0.004 0.813⫾0.029 0.202⫾0.028 1728 1.86 0.247⫾0.005 0.247⫾0.005 0.796⫾0.038 0.214⫾0.037 2252 PRE 61 LOWRI A. DAVIES, GEORGE JACKSON, AND LUIS F. RULL existence of a closed-loop of immiscibility for this system with a single association site, both by GEMC simulation and with the SAFT-VR equation of state. The use of the Wegner expansion of Eq. 共47兲which includes a first correction to scaling, with ␭switch⫽0.1 and Tswitch *⫽1.63, to fit the simulation data, is seen to give an excellent description of the shape of the closed loop. The estimated values for the critical temperatures, together with the values of the coefficients B0and B1obtained are given in Table III. A temperature-composition slice obtained at a higher pressure of P*⫽1.48 by Gibbs ensemble simulation is shown in Fig. 4. As for the case of the P*⫽1.28 slice, the strength of the site-site interaction is fixed at ␧a,b *⫽14.5, and the simulation results indicate the presence of a region of closed-loop immiscibility. No low-temperature immiscible regions are found at this pressure using the GEMC method. For this higher pressure the extent of the closed-loop region has increased from that of the P*⫽1.28 state, which is in line with the theoretical prediction of the effect of pressure on the closed loop 共as described later兲. The prediction of phase behavior obtained with the SAFT-VR equation for two different interaction energies of ␧a,b *⫽14.5 and ␧a,b *⫽15 is also shown. In this case, the theoretical results do not indicate the existence of a LCST for the system at these conditions. The narrow ranges of pressure and interaction strength which result in a region of closed-loop immiscibility is thus illustrated: for a constant interaction strength of ␧a,b *⫽15 an increase of pressure from P*⫽1.28 共Fig. 3兲to P*⫽1.48 共Fig. 4兲the theory predicts the loss of the LCST for the system since the region of liquid-liquid immiscibility merges with the low-temperature gas-liquid coexistence regions. The fit of the simulation data obtained using the Wegner expansion of Eq. 共47兲with ␭switch⫽0.1 and Tswitch *⫽1.70 is also shown in Fig. 4. As for the lower pressure case, this method accurately describes the form of coexistence curve. The corresponding estimates of the critical temperatures and parameters are given in Table IV. In order to monitor the amount of association within the system, we calculate the fraction of molecules bonded Xj. Bonding only occurs when two unlike species are such that their off-center interaction sites overlap. The extent of association for the simulations performed at P*⫽1.28 and P* ⫽1.48 are shown in Fig. 5. For temperatures close to the UCST approximately 10% of the molecules are bonded and a dramatic increase of association is observed as the LCST is approached. This indicates that the low-temperature miscibiltiy of the system below the LCST is directly related to the association of unlike species. The effect of varying the strength of the site-site interaction on the phase behavior of the model system as predicted by the SAFT-VR equation of state is shown in the PT projection of Fig. 6 共also see Ref. 关5兴兲. The vapor pressure curves of the pure components lie on the same line due to the symmetry of the mixture. The dashed curve which originates at the vapor-liquid critical point of the pure component and moves to higher pressures and temperatures represents the gas-liquid critical line. For a site-site interaction energy of ␧AB *⫽15 this gas-liquid critical line goes through a minimum in temperature close to the critical point of the pure component. The critical line then moves to higher pressures, reaching a maximum before ending at the three-phase line. The continuous curve which originates at high pressures is the three-phase line for the mixture, along which two liquid TABLE II. Fraction of molecules bonded in each of the two phases Xaand Xbobtained from NPT Gibbs ensemble Monte Carlo simulations for the symmetrical mixture of square wells with a single interaction site of reduced depth ␧a,b *⫽14.5 at a reduced pressure of P*⫽1.28. See Table I for details. NT*XaXb 1728 1.55 0.423⫾0.030 0.416⫾0.033 1728 1.56 0.297⫾0.028 0.305⫾0.031 1728 1.57 0.308⫾0.030 0.321⫾0.038 1728 1.58 0.246⫾0.018 0.257⫾0.027 1728 1.59 0.242⫾0.026 0.239⫾0.030 1728 1.60 0.220⫾0.028 0.214⫾0.031 1728 1.61 0.193⫾0.020 0.187⫾0.022 1728 1.62 0.185⫾0.023 0.190⫾0.027 1728 1.63 0.181⫾0.025 0.181⫾0.031 1000 1.64 0.182⫾0.034 0.215⫾0.029 1000 1.65 0.157⫾0.031 0.211⫾0.032 1000 1.68 0.170⫾0.030 0.194⫾0.027 1000 1.70 0.164⫾0.032 0.173⫾0.032 1000 1.72 0.162⫾0.029 0.147⫾0.036 1000 1.74 0.134⫾0.026 0.155⫾0.023 1000 1.75 0.156⫾0.027 0.159⫾0.025 1000 1.76 0.150⫾0.030 0.161⫾0.026 1000 1.78 0.141⫾0.027 0.170⫾0.024 1000 1.80 0.130⫾0.025 0.139⫾0.027 1728 1.81 0.141⫾0.023 0.146⫾0.022 1728 1.82 0.140⫾0.023 0.144⫾0.022 1728 1.83 0.146⫾0.016 0.136⫾0.015 1728 1.84 0.137⫾0.017 0.134⫾0.017 1728 1.85 0.140⫾0.019 0.143⫾0.017 1728 1.86 0.140⫾0.021 0.147⫾0.018 TABLE III. Estimated critical constants for the symmetrical square-well mixture with a single bonding site of reduced energy ␧a,b *⫽14.5 at a reduced pressure of P*⫽1.28. The ‘‘error’’ denotes the maximum possible error which is estimated from the combined errors in the compositions obtained from the simulation data. UCST LCST Tc *(N)1.915⫾0.689 1.557⫾0.009 B02.32⫾1.85 3.00⫾0.21 B1⫺2.61⫾3.20 ⫺4.94⫾0.54 TABLE IV. Estimated critical constants for the symmetrical square-well mixture with a single bonding site of reduced energy ␧a,b *⫽14.5 at a reduced pressure of P*⫽1.48. See Table III for details. UCST LCST Tc *(N)2.068⫾1.038 1.466⫾0.239 B02.01⫾2.27 2.03⫾0.89 B1⫺1.45⫾3.01 ⫺1.51⫾1.12 PRE 61 2253CLOSED-LOOP PHASE EQUILIBRIA OF A . . .