Second-order Barker-Henderson perturbation theory for the phase behavior of polydisperse Morse hard-sphere mixture

We propose an extension of the second-order Barker-Henderson perturbation theory for polydisperse hard-sphere multi-Morse mixture. To verify the accuracy of the theory, we compare its predictions for the limiting case of monodisperse system, with predictions of the very accurate reference hypernetted chain approximation. The theory is used to describe the liquid-gas phase behavior of the mixture with different type and different degree of polydispersity. In addition to the regular liquid-gas critical point, we observe the appearance of the second critical point induced by polydispersity. With polydispersity increase, the two critical points merge and finally disappear. The corresponding cloud and shadow curves are represented by the closed curves with 'liquid' and 'gas' branches of the cloud curve almost coinciding for higher values of polydispersity. With a further increase of polydispersity, the cloud and shadow curves shrink and finally disappear. Our results are in agreement with the results of the previous studies carried out on the qualitative van der Waals level of description.


Introduction
A vast majority of industrially important colloidal and polymeric materials are intrinsically polydisperse, i.e., each particle in the system is unique in size, charge, shape, chain length, etc.This feature of colloidal and polymeric systems has a profound effect on their phase behavior and may cause the appearance of new phases and new phase transitions.In addition, the phenomena associated with the phase behavior of polydisperse systems, such as fractionation, are also of technological relevance.Usually, the theoretical methods used to study polydisperse systems treat them as a mixture of an infinite number of components, each characterized by a continuous variable ξ distributed according to a certain distribution function f (ξ), e.g., for a polydisperse hard-sphere fluid, the hard-sphere size σ is usually used as such a variable, i.e., ξ = σ.The theoretical study of the phase behavior of such fluids, using the methods of the modern liquid state theory, represents a nontrivial challenge [1].The main obstacle for theoretical description arises due to the fact that now one has to deal with an infinite number of equations for coexisting phases.One of the possibilities to overcome this obstacle is to resort to the so-called truncatable free energy (TFE) models and combine them with the possibility of their analytical description.TFE models are approximate schemes, where the thermodynamic properties can be expressed by a finite number of generalized moments of the distribution function f (ξ).As a result, the formally infinite number of the phase equilibrium condition equations can be mapped onto a set of a finite number of nonlinear algebraic equations for these moments and solved using standard numerical methods.This route was recently undertaken in a number of studies, i.e., phase behavior of Yukawa and charged hard-sphere polydisperse mixtures were studied using the analytical solution of the mean spherical approximation (MSA) [2][3][4][5] and high temperature approximation (HTA) [6].More recently, HTA and dimer thermodynamic perturbation theory (DTPT) for associating fluids in combination with polymer MSA were used to investigate the phase behavior of a polydisperse mixture of Yukawa chain molecules [7,8].While the approaches based on the analytical solution of the MSA appear to be rather accurate, their application is restricted to the systems with a factorized version of Yukawa interaction, i.e., the matrix of the coefficients describing the strength of interaction is factorized into the product of two vectors.On the other hand, HTA based descriptions, being less accurate, are much more flexible and can be applied to a much larger variety of the potential models [6].
Attempting to improve the accuracy of the HTA based approaches, we present an extension and application of the second-order Barker-Henderson perturbation theory for the phase behavior description of polydisperse multi-Morse hard-sphere mixture.The paper is organized as follows: In section 2 we introduce the model and in section 3 we present a corresponding extension of the BH2 theory.Our numerical results for the phase behavior of the one-Morse version of the model are presented and discussed in section 4, while in section 5 we collect our conclusions.In addition, we include an Appendix with explicit analytical expressions for thermodynamical properties of the multi-Morse hard-sphere polydisperse mixture in question.

The model
We consider the mixture with interparticle pair potential represented by the generalized multi-Morse hard-sphere potential where ξ is the polydispersity attribute, i.e., a continuous version of the species index, σ(ξ) is the hardsphere diameter of the particle of species ξ, σ(ξ, ξ ) = [σ(ξ) + σ(ξ )]/2, z n and 0 are the the screening length and the interaction strength of the Morse potential, respectively.The form suggested for the multi-Morse potential (2.1) is similar to that used earlier for the multi-Yukawa potential [6].This form is very flexible and can be used to model a large variety of realistic potentials by an appropriate choice of the coefficients A nm (ξ) and z n , e.g., in reference [6] it is used to mimic a polydisperse Lennard-Jones mixture.
Here, N M denotes the number of the Morse potential tails and M stands for the number of terms in the sum for one Morse tail.Note that original Morse potential consists of two terms, i.e., one is attractive and the other is repulsive.In our hard-sphere Morse potential repulsive term is substituted by the hardsphere term.
The mixture is characterized by the temperature T [or β = (k B T ) −1 , where k B is the Boltzmann's constant], the total number-density ρ, and by the distribution function f (ξ

Barker-Henderson second-order perturbation theory
To describe thermodynamic properties of polydisperse Morse hard-sphere mixture, we use here the Barker-Henderson second-order perturbation theory.According to this theory, Helmholtz free energy of the system A can be written as a sum of three terms: free energy of the reference system (A ref ) and the two perturbation terms describing the contribution to the free energy due to Morse potential (A 1 , A 2 ): Here, A ref = A HS , where A HS is the free energy of the hard-sphere fluid.The first-order term is: 13605-2 where g (HS) (ξ, ξ ; r ) is the hard-sphere radial distribution function.For the second-order term, we used the macroscopic compressibility approximation (MCA): where ∂ρ/∂p HS = κ HS is the isothermal compressibility of the hard-sphere reference fluid, which is obtained from the Carnahan-Starling equation and given by where the packing fraction η is defined as η = π 6 dξρ(ξ)σ 3 (ξ).Substituting into (3.2) and (3.3) the expression for the potential (2.1), we have ) where G (HS) (ξ, ξ ; z n ) is the Laplace transform of hard-sphere radial distribution function dr r e −z n r g (HS) (ξ, ξ ; r ). (3.7) Here, we use Percus-Yevick approximation for the hard-sphere radial distribution function, since the analytical expressions for its Laplace transform is known.All the rest 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: (3.9) In the above expressions, A HS and µ (HS) (ξ) are calculated using the corresponding Mansoori et al. expressions [9].Within the framework of the BH2 approach, the model in question belongs to the class of 'truncatable free energy models', i.e., the models possessing thermodynamical properties (Helmholtz free energy, chemical potential, pressure) defined by a finite number of generalized moments.In this study, we have the following moments: (3.16) Closed form analytical expressions for thermodynamical properties (Helmholtz free energy, chemical potential, pressure) in terms of the generalized moments (3.10)- (3.16) are presented in the Appendix.

Phase equilibrium conditions
The main obstacle in theoretical studies of the phase behavior in polydisperse systems arises due to the fact that one has to deal with an infinite number of equations for coexisting phases.However, for the present 'truncatable free energy model', these equations can be written as a finite number of equations for the corresponding generalized moments of the distribution function f (ξ) [10].

Results and discussion
In this section we present numerical results for the phase behavior of polydisperse Morse hard-sphere mixture.For the species distribution function f (ξ), we have chosen a log-normal distribution, i.e., where I is the polydispersity index.Log-normal distribution frequently occurs at the colloidal and polymeric processing [1].Note that in the monodisperse limit (I = 1), this distribution is represented by the Dirac delta-function, δ(ξ − 1).On the contrary, when I becomes very large (I 1), the above distribution becomes very wide, increasing hereby the importance of the particles with a large value of ξ.All calculations were carried out for the one-Morse version of the pair potential (2.1), i.e., N M = 1, M = 1, with 13605-5 z 1 = 1.8σ 0 .We consider two types of polydispersity of the model: polydispersity only in the strength (amplitude) of the pair potential A 11 (ξ) and polydispersity in both, amplitude A 11 (ξ) and hard-sphere size σ(ξ).In the former case, we have chosen A 11 (ξ) = A 0 ξ and σ(ξ) = σ 0 , while in the latter case A 11 (ξ) = A 0 ξ and σ(ξ) = σ 0 ξ 1/4 .Here, A 0 = 1 and σ 0 is the hard-sphere size for a monodisperse version of the model at I = 1, which is used as a distance unit.In what follows the density ρ and temperature T are presented in reduced units, i.e., ρ * = ρσ 3 0 and T * = k B T / 0 .
As a first step in our numerical study, we perform the calculation of thermodynamical properties of the monodisperse version of the model, i.e., for I = 1 using the present version of the second-order BH theory and the reference hypernetted chain approximation with the bridge function due to Verlet [11,12].The latter theory is known for being very accurate in predicting the properties of simple fluids [12].The comparison of the results of both theories for the pressure and chemical potential (figure 1) at three different temperatures demonstrate that BH theory is capable of giving relatively accurate results at lower and intermediate densities.At higher values of the densities, the predictions of the BH approach are a bit less accurate.In this case, the cloud and shadow curves coincide, the critical point being located at their maximum.
For a polydisperse system (I 1), for each parent phase density ρ (0) , there is a different binodal.Each binodal is truncated at a maximum temperature, with the corresponding densities, ρ (1) and ρ (2) lying on the cloud and shadow curves, respectively.For a critical value of ρ (0) , ρ (0) = ρ cr , the corresponding binodal passes through the intersection of the cloud and shadow curves.This occurs at T = T cr and, since ρ cr = ρ (1) = ρ (2) = ρ (0) , the point (ρ cr , T cr ) is a critical point, where two coexisting phases, (1) and (2), become identical.We are interested in the phase behavior of the system at relatively large values of poly-  dispersity.For a small polydispersity, the system has only one critical point [2][3][4][5][6], which originates from the regular liquid-gas (LG) critical point of the corresponding monodisperse version of the system.With the polydispersity increase, there appeares an additional critical point induced by polydispersity.This effect on the qualitative van der Waals level of description has been observed by Bellier-Castella et al. [10,13].The second critical point, which we denote as polydisperse (P) critical point, is located at larger values of the density and at lower values of the temperature.It is present for both types of polydispersity studied (figure 2 and 3).With polydispersity increase, both LG and P critical points move towards each other and, for a certain limiting value of polydispersity, they merge.There are no critical points above this limiting value (figure 2, lower panel).With a further increase of polydispersity and at lower temperatures, we expect that the two-phase coexistence becomes unstable and there appeares a region of three-phase coexistence.For relatively high values of polydispersity studied here we observe a rather unusual shape for the cloud and shadow curves.For both types of polydispersity, they are represented by closed curves of elipsoidal-and ∆-like shapes for the shadow curves seen in figure 2 and figure 3, respectively, and closed curves of the linear shape for the cloud curves (figure 2, 3).In the latter case, the 'liquid' and 'gas' branches of the cloud curves almost coincide for the larger polydispersity (figure 2 and 3, lower panels).With a further increase of polydispersity, the cloud and shadow curves shrink and finally disappear.
Finally, in figure 4 we display distribution functions of the polydisperse Morse hard-sphere mixture with size and amplitude polydispersity at two values of the temperature, which are higher and lower than the second critical point temperature, i.e., T * = 1.9 and T * = 1.4,respectively.We present distribution functions of the coexisting phases on the critical binodal at T * = 1.4 (figure 4, upper panel) and on the coexisting cloud and shadow phases at T * = 1.4 (figure 4, intermediate panel), and T * = 1.9 (figure 4, lower panel).As usual [2][3][4][5], on the binodal the particles with larger values of ξ fractionate to the liquid phase while particles with smaller values of ξ fractionate into the gas phase.Fractionation of the particles to the shadow phases for the model at hand depends on the temperature.For the temperatures larger than the temperature of polydispersity induced critical point, we observe fractionation of the usual type, In the upper panel, we display parent f (0) (ξ) (solid line) and daughter f (1) (ξ) (corresponds to gas phase, dashed line) and f (2) (ξ) (corresponds to liquid phase, dotted line) distribution functions for critical binodal at I = 1.055 and T * = 1.4.In the intermediate panel, we show f (1) (ξ) (corresponds to gas shadow phase, dashed line), f (2) (ξ) = f (0) (ξ) (corresponds to liquid and gas cloud phases, solid line) and f (1) (ξ) (corresponds to liquid shadow phase, dotted line) distribution functions at I = 1.055 and T * = 1.4.
In the lower panel, we show the same as in intermediate panel at T * = 1.9.
i.e., the liquid shadow phase contains particles with ξ larger than those of the gas cloud phase while the gas shadow phase contains particles with ξ smaller than those of the liquid cloud phase (figure 4, lower panel).Note that distribution functions for the gas and liquid cloud phases are always the same and equal to the distribution function of the parent phase.The situation changes for the temperatures lower than the temperature of the second critical point (figure 4, intermediate panel).In this case, both liquid and gas shadow phases contain particles with lower values of ξ than those of the liquid and gas cloud phases.At the same time, the liquid shadow phase has particles of a larger value of ξ than the gas shadow phase.
This behavior of the model is related to the fact that for temperatures lower than the temperature of the second critical point, both branches of the shadow curve are located to the left of both branches of the cloud curve, i.e., the density of the shadow phases is always lower than the density of the cloud phases.

Conclusions
In this paper we present an extension of the second-order Barker-Henderson (BH2) perturbation theory for a polydisperse hard-sphere multi-Morse mixture.To verify the accuracy of the theory, we compare its predictions for the limiting case of a monodisperse system, with predictions of the very accurate reference hypernetted chain approximation.The theory is used to describe the liquid-gas phase behavior of the mixture with different type and different degree of polydispersity.In agreement with the previous study [10,13], which was carried out using the qualitative van der Waals level of description, we observe the appearance of the second critical point induced by polydispersity.With the polydispersity increase, the two critical points merge and finally disappear, i.e., for polydispersity larger than a certain threshold value, there are no critical points.The corresponding cloud and shadow curves are represented by the closed curves with 'liquid' and 'gas' branches of the cloud curves that almost coincide with the polydispersity increase.With a further increase of polydispersity, the cloud and shadow curves shrink and finally disappear.

Appendix
Here, we present expressions for thermodynamical properties in terms of the moments (3.10)-(3.16).We have: where  V with respect to the density we get the expression for the chemical potential βµ 1 (ξ):

Figure 1 .
Figure 1.The pressure as a function of the density (the upper panel) and the chemical potential as a function of the density (the lower panel) along three isotherms; the upper set of curves corresponds to T * = 2.5, the intermediate set refers to T * = 2 and the upper set belongs to T * = 1.5.Crosses are RHNC results, solid lines correspond to the Barker-Henderson second-order results.

Figure 2 .
Figure 2. Phase diagrams of the polydisperse Morse hard-sphere mixture with amplitude polydispersity only in the (ρ * , T * ) plane for three different values of the polydispersity index I , I = 1.02967 (the upper panel), I = 1.034 (the intermediate panel) and I = 1.035 (the lower panel), obtained using the Barker-Henderson second-order perturbation theory, which includes cloud (solid line) and shadow (dotted line) curves, two critical points and critical binodals (dashed lines).Filled circle denotes the position of the critical points.

Figure 3 .
Figure 3. Phase diagrams of the polydisperse Morse hard-sphere mixture with size and amplitude polydispersity in the (ρ * , T * ) plane for three different values of the polydispersity index I , I = 1.035 (the upper panel), I = 1.045 (the intermediate panel) and I = 1.055 (the lower panel), obtained using the Barker-Henderson second-order perturbation theory, which includes cloud (solid line) and shadow (dotted line) curves, two critical points and critical binodals (dashed lines).Filled circle denotes the position of the critical points.