Liquid-gas phase behavior of polydisperse dipolar hard-sphere fluid: Extended thermodynamic perturbation theory for central force associating potential

The liquid-gas phase diagram for polydisperse dipolar hard-sphere fluid with polydispersity in the hard-sphere size and dipolar moment is calculated using extension of the recently proposed thermodynamic perturbation theory for central force (TPT-CF) associating potential. To establish the connection with the phase behavior of ferrocolloidal dispersions it is assumed that the dipole moment is proportional to the cube of the hard-sphere diameter. We present and discuss the full phase diagram, which includes cloud and shadow curves, binodals and distribution functions of the coexisting daughter phases at different degrees of the system polydispersity. In all cases studied polydispersity increases the region of the phase instability and shifts the critical point to the higher values of the temperature and density. The larger size particles always fractionate to the liquid phase and the smaller size particles tend to move to the gas phase. At relatively high values of the system polydispersity three-phase coexistence is observed.


I. INTRODUCTION
In this paper we consider the liquid-gas phase behavior of polydisperse dipolar hard-sphere mixture. Recently, liquid-gas phase equilibria in monodisperse dipolar hard-sphere fluid, Yukawa dipolar hard-sphere fluid and Shtockmayer fluid were studied using thermodynamic perturbation theory for central force (TPT-CF) associating potential [1,2,3,4]. In this study we propose extension of the TPT-CF, which enables us to investigate the phase behavior of polydisperse mixture of the dipolar hard spheres with polydispersity in both hard-sphere size and dipole moment. We call this extension as extended TPT-CF (ETPT-CF). In a certain sense ETPT-CF can be seen as a combination of Wertheim's TPT [5,6] for associating fluid with association taking place due to off-center attractive sites and TPT-CF, which allows multiple bonding of each site. In our theory we have several Wertheim's type of associating sites with the possibility for each site to be bonded to several other sites (in Wertheim's TPT each site is only singly bondable). Final expressions for thermodynamical properties of polydisperse dipolar hard-sphere fluid is written in terms of the finite number of distribution function moments, i.e. in the framework of ETPT-CF this system belongs to the family of the so-called truncatable free energy models. This property enables us to calculate the full liquid-gas phase diagram (including cloud and shadow curves and binodals) and to study effects of fractionation on the level of the distribution functions of coexisting daughter phases.

A. Analysis and classification of diagrams
We consider a multicomponent fluid mixture consisting of n species with a number density ρ = n a ρ a at a temperature T (β = 1/k B T ), where ρ a is the density of the particles of a species. The particles of the species a and b interact via the pair potential U ab (12), which can be written as a sum of the reference U ab ref (12) and associative U ab ass (12) parts U ab (12) = U ab ref (12) + U ab ass (12), (1) where 1 and 2 denote positions and orientations of the particles 1 and 2. We assume that associative part of the potential can be represented as a sum of M a ×M b terms, i.e.
U ab ass (12) = KL U ab KL (12), (2) where the lower indices K and L take the values A, B, C, . . . , respectively. These values specify splitting of the total associating potential U ab ass (12) into several particular pieces. For example in the case of the models utilized by Wertheim [6] these indices denote off-center attractive sites and in the case of the Mersedes-Benz (MB) type of the models [7] or cone models [8] they stand for the type of the hydrogen bonding arms. Hereafter we will refer to these indices as to the site indices, keeping in mind that they may have more general meaning. Here M a and M b are the number of such sites on the particles of a and b species, respectively. According to (1) and (2) Mayer function f ab (12) for the total potential (1) takes the following form: f ab (12) = f ab ref (12) + e ab ref (12) KL 1 + f ab KL (12) − 1 , where we are using the usual notation: e(12) = exp [−βU (12)], f (12) = e(12) − 1, For the sake of diagrammatic analysis we will follow Wertheim [6] and instead of circles we introduce hypercircles to represent particles in diagrammatic expansions. Herez a (i) = z a exp [−βU a (i)], i denotes position and orientation of the particle i, and U a (i) is an external field. For a uniform systemz a (1) ≡ z a . Following [1,2,5,6] we introduce the definition of the s-mer diagrams. These are the diagrams consisting of s hypercircles, which all are connected by the network of f KL bonds. The site circles, which are incident with more than m a K f ab KL bonds are called oversaturated site circles. We consider now the set of oversaturated site circles with each pair connected by at least one path formed by the circles from the same set.
The subdiagram involving this set of circles, together with the site circles adjacent to them and f KL bonds connecting all these circles , we call the oversaturated The procedure for obtaining the expression for ρ a (1) from ln Ξ remains unchanged.
The diagrams appearing in thez expansion of the singlet density ρ a (1) can be classified with respect to the number of f ab KL bonds associated with the labeledz a (1) hypercircle. We denote the sum of the diagrams with i K ≤ m a K associating bonds connected to the site K (K = A, B, C, . . .), which belongs to the particle of species a as ρ a A i A ,B i B ,C i C ,... (1). Any site K, which is connected to i K > m a K associating bons, will be denoted as K m a K In what follows we will use also a condensed version of the notation, i.e.
where {i} = i A , i B , i C , . . .. The set {i} with all indices, except one index i K , equal 0, will be denoted as i K , i.e. {i} = 0, . . . , 0, i K , 0, . . . , 0 ≡ i k , so that for any quantity Thus ρ a (1) can be written as follows B. Topological reduction The process of switching from the activity to a density expansion goes in the same fashion as in Refs. [1,2,6]. However, to proceed it is convenient to use an operator form of notation. The operators are introduced in a manner similar to that presented in references [1,6] to which we refer the reader for more details.
We associate with each labeled l hypercircle an operator ǫ a {i} (l) with the following properties: . .. The one-point quantities, which, for convenience, are denoted by x a {i} , can be presented as illustrated below: The operators ǫ a {i} are straightforward generalization of the operators introduced earlier [1, 6,9]. Thus, the rules of manipulation with the new quantitiesx a are similar to that discussed before. In particular, the usual algebraic rules apply to these quantities and analytical functions ofx a are defined by the corresponding power series. Similar, as in references [1,6,9], it is convenient to use the angular brackets to specify the operation In the case of several labeled circles the subscripts on the brackets denote the circle to which procedure (10) is to be applied.
Analyzing the connectivity of the diagrams in ρ a (1) at a labeledz a (1) hypercircle we haveρ This relation can be inverted expanding into a power series, i.e.
ρ a (1) =σ a (1) Now the diagrammatic expansions for c a {i} can be expressed in terms of irreducible diagrams. To present this result in compact and convenient form we introduce a sum of the diagrams c (0) defined as follows: C. Extended thermodynamic perturbation theory for central force associating potential Now we are in a position to rewrite the regular one-density virial expansion for the pressure P in terms of the density parametersσ a (1). Following the scheme, proposed earlier [2,4,5,6] we have expression for the pressure in operator form and explicitly where V is the volume of the system. Similarly, as in [5,6] one can verify that these expressions satisfy the regular thermodynamic relationρ a = β∂P/∂µ a , wherē ρ a = ρ a (1) d(1) and µ a is the chemical potential. This can be achieved by taking a variation of (15) (or (16)) and combining (11), (13) and (14). The corresponding expression for Helmholtz free energy is where Λ a is the thermal de Broglie wavelength. This expression is derived using the regular thermodynamic expression for Helmholtz free energy A = a N a µ a − P V together with relation which follows from (11), written for ρ a {0} . Here N a is the number of particles of species a in the system. Helmholtz free energy in excess to its reference system value A ref is obtained by subtracting corresponding expression for A ref from (17) ref is the corresponding sum of the diagrams for the reference system. Ordering the virial expansion (19) with respect to the number of associating f KL bonds and neglecting the terms with more than one associating bond we have where g ab ref (12) is the reference system distribution function and Due to single bond approximation c a {i} = 0 for all values of the set {i}, except for {i} = 0 and {i} = i K with i K = 1. This property together with (11) yield the following relations: and The set of relations (20), (21) and (23) defined all the quantities needed to calculate the Helmholtz free energy of the system (19), provided that the properties of the reference system are known.
Finally it is worth noting, that ETPT-CF theory developed here reduces to TPT1 proposed by Wertheim [6], if for all sites single-bonding condition m a K = 1 is assumed. In the other limiting case of only one site per particle ETPT-CF will coincide with TPT-CF developed earlier [1,2,4].

D. Extended TPT-CF for two sites with double-bonding condition
The theory presented in the previous section is quite general and can be applied to a number of different situations. However in the present study we are interested in the version of the theory for the model with two sites both of which can be bonded twice. More specifically we are interested in the extension and application of the theory to study the phase behavior of polydisperse dipolar hard-sphere mixture.
We assume that each of the particles in the system have two doubly bondable attractive sites, A and B, i.e. we have: M a = 2 and m a A = m a B . We also assume, that attractive interaction is acting only between the sites of the same sort. Using these suggestions, relations (11) and (12) and taking into account that the system is uniform, the density parameters σ a where K takes the values A and B and κ a K 1 = 1 + c a K 1 . These two equations give In turn, using (21), for κ a K 1 we have where Combining (20), (29) which satisfy the following set of equations: Chemical potential ∆µ a and pressure ∆P in excess to their reference system values can be obtained using standard thermodynamical relations: Finally, the average size of the clusters, which appear in the system, can be characterized by the average length of the chain L K formed by either A-bonded particles. This quantity is defined by the following where α a K,end is the fraction of singly K-bonded particles (fraction of the chain ends), α a K,mid is the fraction of doubly K-bonded particles (fraction of the chain middles) and α a 0 is the fraction of nonbonded particles. For these fractions we have Substituting these expressions into expression for L A (35) and using (28), (29) and we get the final expression for L A in terms of κ a K 1 : where if K = A thenK = B and if K = B thenK = A.

HARD-SPHERE FLUID
A. The model We consider polydisperse dipolar hard-sphere fluid mixture with a number density ρ and polydispersity in both hard-sphere diameter σ and dipolar moment d µ .
We assume, that the dipole moment is proportional to the particle volume, i.e.
Thus the type of the particle is completely defined by its hard-sphere size and hereafter we will be using σ instead of the indices a, b, . . . to denote the particle species. We also assume, that hard-sphere size of the particles is distributed Interaction between particles of species σ 1 and σ 2 in our system is described by the following pair potential: where U hs (r, σ 1 σ 2 ) is the hard-sphere potential and U dd (r, σ 1 σ 2 ) is the dipole-dipole potential, given by Here ϕ 1 and ϕ 2 denote the angles between the dipole vectors and vector, which joins the centers of the two particles, and φ 1 and φ 2 are the azimuthal angles about this vector. To proceed we have to split the total potential (40) into the reference and associative pieces. We assume, that the reference part of the potential is represented by the hard-sphere part U hs (r, σ 1 σ 2 ) and associative part by the dipoledipole potential U dd (r, σ 1 σ 2 ). At the contact distance σ 12 = (σ 1 + σ 2 )/2 the latter potential has two equal potential minima of the depth −2d µ (σ 1 )d µ (σ 2 )/σ 3 12 at "noseto-tail" configuration (ϕ 1 = ϕ 2 = 0, ϕ 1 = ϕ 2 = π. These minima are responsible for the formation of the chains of the particles in the system. In addition there are twice less deep minima (−d µ (σ 1 )d µ (σ 2 )/σ 3 12 ) at antiparallel configuration with ϕ 1 = ϕ 2 = π/2, φ 1 − φ 2 = π. The latter minima cause the formation of the network connecting the chains. According to earlier theoretical and computer simulation studies [10,11] competition between chain formation and network formation defines the existence of the liquid-gas phase transition in the dipolar hard-sphere fluid. To account for this effect we propose the following splitting of the total associative potential U ass (12, σ 1 σ 2 ) = U dd (12, σ 1 σ 2 ): where Θ(ϕ) = H (π/2 + ϕ 0 − ϕ) H (π/2 − ϕ 0 + ϕ) and H(x) is the Heaviside step function. Here ϕ 0 plays a role of the potential splitting parameter. For ϕ 0 = π/2 we have that U BB (12, σ 1 σ 2 ) = U dd (12, σ 1 σ 2 ) and U AA ( where we have dropped the lower index 1, i.e. κ K 1 (σ) ≡ κ K (σ). In order to solve this equation we propose here to interpolate the key quantity of the theory, the volume integral I KL (σ 1 σ 2 ), using a sum of N Y Yukawa terms. Since the reference system pair distribution function g ref (r, σ 1 σ 2 ) is independent of mutual orientation of the particles for the integral (31) we have wheref KK (r, σ 1 σ 2 ) is orientation averaged Mayer function for associative potential U KK (12, σσ). We assume, thatf KK (r, σ 1 σ 2 ) can be represented in the following Parameters A (n) K (σ) and z (n) K are obtained using the interpolation scheme, which is presented and discussed in the Appendix A. Using (45) and (46), we have where G (n) dr is the Laplace transform of the radial distribution function g ref (r, σ 1 σ 2 ). We will be using here Percus-Yewick approximation for hard-sphere radial distribution function, since analytical expression for its Laplace transform is known [12] e −z (n) K σ 12 G (n) where D (n) Here m l are the moments and m (n) K,l are the generalized moments of the distribution function F (σ). Symbolically expression for these moments can be represented as follows: Hereafter all the quantities denoted as m with certain set of indices will represent the generalized moments defined by (51). Corresponding expressions for m(σ) are collected in the Appendix B. Inserting (47) into equation (44) and using (48), we where C (n) K,j satisfies the following set of equations: Here Ω (n) Thus solution of the integral equation (44) for the unknown function κ K (σ) now is reduced to solution of a set of equations for 4N Y unknown constants C (n) K,j . This solution can be used to calculate κ (n) K (σ), which in turn can utilized to calculate thermodynamical properties of the system via Helmholtz free energy (32). Generalizing expression for Helmholtz free energy (32) for polydisperse system, we have Now we can use the standard relation between Helmholtz free energy and chemical potential (34), generalized to polydisperse case where δ/δ {F (σ)} denote functional differentiation with respect to the distribution F (σ). We find Here µ (n) functional derivatives δP where δC Using this expression together with expression for the chemical potential (60), we where J (n) Expression for the integral J where µ (ex) ref (σ) is the reference system chemical potential in excess to its ideal gas value.

C. Phase equilibrium conditions
One can easily see that the thermodynamical properties of the model at hand K,l . Thus the polydisperse mixture of dipolar hard spheres treated within ETPT-CF belongs to the class of truncatable free energy models [14]. This property allows us to map the phase coexistence relations onto a set on nonlinear equations for the unknown moments of the daughter distribution functions.
We assume that at a certain density ρ 0 and composition F 0 (σ) the system separates into two phases with the densities ρ 1 and ρ 2 , and compositions F 1 (σ) and F 2 (σ). Hereafter the lower index 0 refer to the parent phase and the lower indices 1 and 2 refer to the daughter phases. At equilibrium these quantities take the values, which follows from the phase equilibrium conditions, i.e.: (i) conservation of the total volume of the system, (ii) conservation of the total number of the particles of each species, (iii) equality of the chemical potentials of particles of the same species in the coexisting phases, (iv) equality of the pressure in the coexisting phases. These conditions finally lead to the following set of relations [15,16]: where µ α ) (ex) is the excess (over the ideal gas) chemical potential of the particle of species where K = A, B and {X α } represent unknowns of the problem, i.e.
The remaining 2 + 8N Y equations follows from the equality of the pressure in coexisting phases (73), from the set of equations for C (n)(α) K,1 and C (n)(α) K,2 (53) and from the normalizing condition (74) for either phase α = 1 or α = 2, Solution of the set of equations (53), (77)-(83) for a given density ρ 0 and distribution function F 0 (σ) of the parent phase gives the densities ρ α and distribution functions F α (σ) of the two coexisting daughter phases. The coexisting densities at different densities of the parent phase ρ 0 defined binodals, which are terminated when a density of one of the phases is equal to the parent phase density ρ 0 . These termination points form the cloud and shadow coexisting curves. These curves intersect at the critical point, which is characterized by the critical density ρ cr = ρ 1 = ρ 2 = ρ 0 and critical temperature T cr . The cloud-shadow curves can be obtained as a special solution of the general coexisting problem, when the properties of one phase are equal to the properties of the parent phase: assuming that the phase α = 2 is the cloud phase, i.e. ρ 2 = ρ 0 , and following the above scheme we will end up with the same set of equations (53), (77)-(83), but with ρ 2 and F 2 (σ) substituted by ρ 0 and F 0 (σ), respectively.

IV. RESULTS AND DISCUSSION
In this section we present our numerical results for the liquid-gas phase diagram of polydisperse dipolar hard-sphere fluid at different degree of polydispersity. For size distribution function F (σ) we have chosen the beta distribution given by where Here σ u and σ d defined the range of values for σ. In our calculations we have chosen σ u = 1.2947σ 0 and σ d = 0.85σ 0 . In

V. CONCLUSIONS
In this paper we propose extension of our TPT-CF approach to account for several associating sites with the possibility of each site to be multiply bonded.
The theory is applied to study the liquid-gas phase behavior of polydisperse dipolar hard-sphere fluid. We present a full phase diagram, which include cloud and shadow curves, binodals and distribution functions of the coexisting phases.

VI. APPENDIX A
Orientationally averaged Mayer functionf KK (r, σ 1 σ 2 ) were fitted empirically as a sum of Yukawa-like terms where f 1 =f AA and f 2 =f BB , i, j denotes the species of the particles, σ ij = (σ i + σ j ) /2, N = 6, and A n (σ, T ) are parameters that depend on particle size and temperature and were fitted by "polynimally-exponential" functions of different forms for first and second integrals. (1) n,i,4 y 3 , f or i = 5, 6, 7, 8, 9, (90) where y = T min T and T min = 0.13.
The x and y are the same as for the first integral.
The fitting procedure consisted of finding suitable b n,i,j and ω n,j parameters for first and second integrals by means of differential evolution optimization algorithm [17,18].