Closed-loop liquid-liquid immiscibility in mixture of particles with spherically symmetric interaction

Thermodynamic perturbation theory for central-force (TPT-CF) type of associating potential is used to study the phase behavior of symmetric binary mixture of associating particles with spherically symmetric interaction. The model is represented by the binary Yukawa hard-sphere mixture with additional spherically symmetric square-well associative interaction located inside the hard-core region and valid only between dissimilar species. To account for the change of the system packing fraction due to association we propose an extended version of the TPT-CF approach. In addition to the already known four types of the phase diagram for binary mixtures we were able to identify the fifth type, which is characterized by the absence of intersection of the $\lambda$-line with the liquid-vapour binodals and by the appearance of the closed- loop liquid-liquid immiscibility with upper and lower critical solution temperatures.


Introduction
According to the Gibbs phase rule, the binary mixture may have up to four coexisting phases simultaneously. This fact implies that the phase behavior of the binary fluid could be very rich and complicated. Systematic study and classification of the peculiarities of the binary systems phase diagrams topologies has been undertaken more than 40 years ago by Scott and van Konynenburg [1,2]. These studies are based on the application of the van der Waals equation of state, which in most cases is capable of providing qualitatively correct description of the phase behavior. Most of the subsequent studies, carried out using quantitatively more accurate methods of the modern liquid state theory [3], have been focused on the investigation of the phase behavior of symmetric binary mixtures [4][5][6][7]. These are the mixtures with identical interaction between particles of the similar species and different interaction between particles of the dissimilar species. Phase behavior of such mixtures is defined by the competition between gasliquid and liquid-liquid coexistence. Combining Monte-Carlo computer simulation and theoretical meanfield calculations, Wilding et al. [4] identified three types of a phase diagram for square-well hard-sphere symmetrical binary fluid mixture. Similar three types of a phase diagram were detected in the symmetrical binary hard-sphere Yukawa mixture using self-consistent Ornstein-Zernike approximation (SCOZA) and Monte-Carlo computer simulation method [5,6]. At the same time, the phase diagram of the fourth type was also detected using SCOZA approach [8]. More recently, the first-order thermodynamic perturbation theory has been used to study the phase behavior of the binary Yukawa mixture with asymmetry in hard-sphere sizes [7].
In this study we are focused on the investigation of the phase behavior of symmetric Yukawa hardsphere binary mixture with additional spherically symmetric square-well associative interaction located inside the hard-core region. This additional interaction is valid only between particles of dissimilar species. The hard-sphere version of the model has been developed and studied by Cummings and Stell [9] and by Kalyuzhnyi et al. [10]. Originally, this version of the model was used as a simple hamiltonian model of the chemical reaction [9]. On the other hand, the model of this type can be regarded as a coarse grained version of the model for sterically or charge stabilized colloidal dispersions, protein solutions, star-polymer fluids, etc. [11][12][13]. Effective interaction between macroparticles of such systems has an attractive potential well at short distances and a repulsive potential mound at intermediate distances. The phase diagrams, which include two-phase gas-liquid, liquid-liquid diagrams and three-phase gas-demixed liquid diagrams have been calculated using thermodynamic perturbation theory for central force associating potential [14][15][16].
Phase behavior of the model, which is similar to the present one has been studied earlier by Jakson [17]. His model is represented by the symmetrical binary hard-sphere mixture with mean-field type of attractive interaction valid only between the same species and orientationally dependent associative interaction between particles of unlike species. Associative interaction appears due to the off-center squarewell sites. This model was used as a generic model for the phase behavior description of the binary mixtures with the possibility of hydrogen bond formation between unlike species. Combining Wertheim's TPT for associating fluids and mean-field approach, Jackson was able to show that for a certain choice of the potential model parameters, the system exhibits the closed loop liquid-liquid immiscibility with the upper and lower critical solution temperatures. It was concluded that closed loop coexistence appeares due to the presence of the highly anysotropic attraction between off-center bonding sites.
In the present work we demonstrate that the systems with spherically symmetric interaction may also have closed loop liquid-liquid immiscibility. The paper is organized as follows. In section 2 we describe the model to be considered and in section 3 we present and discuss details of the TPT-CF theory, specialized to the model at hand. Our results and discussion are included in section 4 and our conclusions are collected in section 5.

The model
We consider symmetric binary Yukawa hard-sphere mixture with additional associative interaction between dissimilar particles. The total pair potential of the model U i j (r ) is represented as a sum of hardsphere Yukawa potential U HSY i j (r ) and associative potential U ass (r ), i.e.: where the lower indices i , j denote the species of the particles and δ i j is the Kroneker delta. In our symmetric binary system, Yukawa interaction between particles of the same species is the same, i.e., U HSY 11 (r ) = U HSY 22 (r ), and between the particles of dissimilar species it is regulated by the parameter α (0 < α < 1), i.e. U HSY 12 (r ) = αU HSY 11 (r ). We have: 11 , z n and ǫ 0 are the screening length and the interaction strength of the Yukawa potential, respectively, d i j = (d i + d j )/2, d i is the hard-sphere diameter. We consider the system with hard spheres of equal size, i.e., d 1 where L is the bonding distance, ω and ε ass are the square well potential width and depth, respec- tively. In what follows we will consider the hard-sphere Yukawa potential (2.2) and (2. (r ) f ass (r ) dr. (2.5) The type of the clusters, which will be formed in the system due to association depends on the value of the bonding length L [10,14]. For values of L lying in the interval (0, d/2) only dimers can be formed in the system. When d/2 < L < 3d, the formation of the chains is possible. Further increase in L leads to an increase in the maximum number of the particles of the type 1 (or 2), which can be simultaneously associated with the particle of the type 2 (or 1), and the formation of the branched chains will be possible.

Theory
To describe thermodynamic properties of the model at hand we will utilize here thermodynamic perturbation theory for central force associative potential (TPT-CF) [14][15][16]. According to TPT-CF, Helmholtz free energy of the system A can be written as a sum of two terms: free energy of the reference system A ref and the term describing the contribution to the free energy due to association A ass : Here, A ref = A HSY , where A HSY is the free energy of the hard-sphere Yukawa fluid. To calculate A HSY , we are using the high temperature approximation. All the rest of thermodynamical quantities can be obtained using the expression for Helmholtz free energy (3.1) and standard thermodynamical relations, e.g., differentiating A with respect to the density, we get the expression for the chemical potential: and the expression for the pressure P of the system can be calculated invoking the following general relation:

High temperature approximation
Under the high temperature approximation, the expression for the free energy is: where G HS (z n ) is the Laplace transform of hard-sphere radial distribution function dr r e −z n r g HS (r ).
Differentiating the expression for Helmholtz free energy (3.4) with respect to the density, we get the following expression for the chemical potential: where µ HSY k is the hard-sphere chemical potential and Pressure P HSY of the system can be calculated invoking the following general relation: In the above expressions, A HS and µ (HS) k are calculated using the corresponding Carnahan-Starling expressions [18].

Thermodynamic perturbation theory
According to the TPT-CF for the associative part of the free energy A ass , we have: k n for l = 2, . . . , m . (3.14) Here, m is the maximum number of associative bonds per particle (the maximum number of the particles, which can be bonded to a given particle simultaneously), σ (l ) k is the density of l-times bonded particles. For the present two-component mixture, the density parameters σ (0) k and σ (1) k satisfy the following set of equations

43606-4
where K = 4π y (00) 12 (r )e (HSY) (r ) f ass (r )r 2 dr = 4πBL 2 y (00) 12 (L), (3.16) y (00) 12 (r ) represent the cavity distribution function between two Yukawa hard spheres of species 1 and 2 infinitely diluted in the original associating fluid in question. Usually, this function is approximated by the hard-sphere Yukawa cavity correlation function y (HSY) 12 (r, η) calculated for the packing fraction η. This appears to be a good approximation for the models with bonding length L ≈ d, since in this case y (00) 12 (r ) only weakly depends on the degree of the system association. However, for L < d, the actual (effective) packing fraction η eff and thus y (00) 12 (r ) are strongly dependent on the system degree of association, and the usual approximation becomes inadequate. In the present study we propose to approximate y (00) 12 (r ) by the hard-sphere Yukawa cavity correlation function y (HSY) 12 (r, η eff ) calculated for the effective packing fraction η eff , i.e., where the excluded volume V exc is: According to this expression, η eff [and thus y (HSY) 12 (r )] depends on the degree of association of the system represented by the fractions of free X (0) k and singly bonded X (1) k particles.
In the present study, the solution of this equation is obtained via numerical iteration method. On each iteration step, the new estimate for the fractions X (l ) k,new (l = 0, 1) is calculated by solving the following set of equations: which is obtained using a set of equations (3.15). Here, X old is the value of X calculated during the previous iteration step. Our iteration loop consists of two steps. In the first step, the current value of η eff is used to calculate the new values of X (l ) k using the set of equation (3.19). On the second step, we insert these values of X (l ) k into the right-hand side of the relation (3.17) to get a new estimate for η eff . This iteration loop is repeated until the following condition is satisfied. For the initial guess we have used the value of η eff = η.

The cavity correlation function for Yukawa hard sphere fluid
The cavity correlation function y (HSY) i j (r ). Here, the upper indices (HS) and (HSY) denote the hard-sphere and hardsphere Yukawa quantities, respectively, and h and c denote total and direct correlation functions, respectively. In the hard-core region δh (HSY) i j (r ) = 0 and for δc (HSY) i j (r ), we have used the expression obtained in the framework of the first-order mean spherical approximation [19]. The hard-sphere cavity correlation function y (HS) i j (r ) was calculated using Henderson-Grundke approximation [20]. Closed form analytical expressions for δc (HSY) i j (r ) and y (HS) i j (r ) are presented in the appendix.

Calculation of the phase diagram
Our calculation of the phase diagram follows closely the scheme, proposed in [5]. It is based on the solution of the set of equations that follow from the conditions of phase equilibrium, i.e., equal chemical potentials and pressures of the coexisting phases at a given temperature. The coexisting phases are characterized by (ρ, x) and (ρ ′ , x ′ ). From the Gibbs' phase rule, we expect up to four phases to be in equilibrium, i.e., the vapour (V), the mixed fluid (MF), and two (symmetric) phases of the demixed fluid (DF).
The V-MF transition is obtained by solving the set of equations:   gives the density ρ of the V or MF and the density of the DF with concentrations x(ρ ′ ) and 1 − x(ρ ′ ), in equilibrium.

Results and discussion
In this section we present our numerical results for the phase behavior of the model in question.
All the calculations are carried out at Yukawa screening parameter z n d = 1.8 and square-well width ω = 0.0000404981.
According to the previous studies [10], predictions of our theory for thermodynamical properties of the model with A i j = 0 are in a good agreement with computer simulation predictions. To test the accuracy of the theory for the model with ǫ * ass = 0, we compare theoretical and computer simulation predictions for its phase behavior. In figure 1 we show the phase diagram of the system at ǫ * ass = 0 and three values of α, i.e., α = 0.65, 0.7, 0.75. These are the system parameters for which the three types of the phase diagram were identified [4,5], depending on the position of intersection point of the λline, which represent the second-order demixing transition, with the binodals of the liquid-vapour (LV) phase transition. In addition, for comparison in the same figure, we present the corresponding computer simulations results [6]. Overall there is a reasonably good qualitative agreement between theoretical and computer simulation predictions. Predictions of the theory in the region of the LV critical point are 43606-6 about 6% higher than those of the computer simulation. As a result, while the types I and II of the phase diagrams (according to the nomenclature of references [4,5]) are theoretically reproduced for the set of the potential model parameters used to simulate the type III of the diagram, theoretical calculations still show the type II of the diagram with a small portion of stable binodals in the vicinity of the LV critical point. However, it is quite obvious that a small decrease in α will cause the theoretical phase diagram to change its type from type II to type III. This can be seen in figure 2 (panel a), where our results for α = 0.63 are shown. In the phase diagram of the type I the LV, coexistence is unstable with respect to  the three-phase coexistence between mixed fluid (MF) and demixed fluid (DF) and the λ-line ends at the tricritical point. In the type II of the diagram, λ-line ends also at the tricritical point, however, there is a portion of the LV phase diagram beeing stable in the range of the temperatures between the critical temperature and the temperature of the triple point, where one can observe LV coexistence at lower densities. Between tricritical temperature and temperature of the triple point, the MF-DF three-phase coexistence can be seen. At the triple point, there is a four-phase vapour, MF and DF coexistence. In the case of the type II diagram, the λ-line intersects the liquid binodal at the temperature slightly below the critical. In the type III of the diagram, the λ-line intersects the liquid binodal at the temperature well below the critical temperature. Here, we can see the critical end point below which there is a three phase V-DF coexistence and above which (up to the LV critical temperature), there is a LV coexistence. In the type IV of the diagram, the λ-line intersects LV binodals at the densities that are lower than the LV critical densities [8]. This occures at α = 0. We have also detected this type of the diagram, using the current approach; however, the results are not shown here.
Next, we proceed to the discussion of the phase diagrams for the nonzero value of the strength of associating interaction ǫ * = 5.2, 6.0, 6.5 at α = 0.63 (figures 3-5). Unfortunately, computer simulation results for the phase behavior of the model at hand are not available. However, taking into account a reasonable performance of the theory in the two limiting cases discussed above (A i j = 0 and ǫ 0 = 0), we expect that the accuracy of the theory for the full version of the model will be satisfactory as well. In

43606-8
Closed-loop liquid-liquid immiscibility the temperature decrease, the difference in the compositions of the coexisting liquids increases. With the increase of the strength of association ǫ * , the topology of the phase behavior in T * vs ρ * coordinate frame changes from the type III (figure 2) to type II at ǫ * = 5.2 (figure 3) and next to type I at ǫ * = 6.0 (figure 4). At the same time, one can observe the appearance of the closed loop liquid-liquid immiscibility curves with the upper stable and lower unstable critical solution points (figures 3 and 4). The stable portion of the curves increases with an increasing strength of associative interaction. Finally, for ǫ * = 6.5, the closedloop coexistence curves for demixing coexistence becomes stable ( figure 5, panel b). This corresponds to the situation when there is no intersection between LV binodals a λ-line ( figure 5, panel a). Thus, for these values of the potential model parameters, in addition to already identified four types of the phase diagram, we have identified one more, which we call the type V of the two-component mixture phase diagram topology. In a future we are planning to extend and apply our approach to study the effects of the external field [21] and porous media [22,23] on the phase behaviour of the current model.

Conclusions
In this paper we have used the TPT-CF approach to study the phase behavior of a symmetric twocomponent Yukawa mixture of associating particles with spherically symmetric interaction. Our theoretical predictions for the phase diagram of the version of the model without association appear to be in reasonably good qualitative agreement with the predictions of the corresponding Monte-Carlo computer simulation method [6]. For the model with nonzero associating potential, we were able to identify, in addition to the already known three types of the phase diagram topologies [4,5], the type V of the phase diagram. This type is characterized by the absence of intersection of the λ-line, which represents a demixing coexistence, with LV binodals. As a result, the stable closed-loop liqui-liquid immiscibility curves with upper and lower critical solution temperatures can be observed for the larger values of the temperature and density. Thus, closed-loop liquid-liquid immiscibility, which was observed earlier for the binary systems with highly directional attractive forces [17], can be also seen for the binary fluids with spherically symmetric interaction.

A. Grundke-Henderson approximation
To calculate the hard-sphere cavity correlation function y (HS) i j (r ) we use Grundke-Henderson approximation [20]. For r < d, we have: ln y (HS) i j (r ) = 3 n=0 a n r n , (A.1)