Liquid-gas phase behavior of Stockmayer fluid with high dipolar moment

A multi-density version of the thermodynamic perturbation theory for associative fluids with central force associative potential is developed. The theory is used to describe the liquid-gas phase behavior of Stochmayer fluid with different values of the dipole moment ranging from low to high. Results of the theory are in a reasonable agreement with corresponding computer simulation results in a whole range of the dipole moment values studied. Contrary to previous computer simulation findings our calculations predict the occurrence of the liquid-gas phase coexistence for the soft-sphere dipole fluid model.


Introduction
Stockmayer fluid is a simple and convenient model, which is often used to study the properties of the polar fluids [1] and ferrofluids [2][3][4], in particular their phase behavior.Perhaps first theoretical predictions for the liquid-gas phase diagram of Stockmayer fluid was given by Stell, Rasaiah and Narang (SRN) [5,6].Their calculations were carried out using the Padé approximation based on the perturbation theory due to Pople [7] and Stell et al. [5,6].Comparison of the theoretical results with corresponding computer simulation results, carried out later [8,9], shows that predictions of the theory are accurate in the region of relatively low values of reduced dipole moment µ * d 2 2, where µ * d = µ d / LJ σ 3 LJ , µ d is the dipole moment and σ LJ and LJ are the Lennard-Jones (LJ) potential parameters.In the region of higher values of µ * d the theory becomes substantially less reliable.According to computer simulation studies [10] in this region, strong dipole-dipole interaction causes the particles to form chains in "nose-to-tail" configuration.Similar behavior has been also detected for the dipolar hard-sphere fluid [11].Thus, for the theory to be successful, adequate account for the chain formation is crucial.Recently we have developed thermodynamic perturbation theory (TPT), which explicitly takes into account formation of the chains due to strong dipole-dipole interaction [12].The theory is based on the multi-density approach for associating fluids with central force type of associating potential, which was developed earlier by one of the authors [13].In turn the latter approach represents an appropriate reformulation of Wertheim's multi-density theory for associating fluids [14,15].We refer to our version of the TPT as TPT central force (TPTCF).
TPTCF proposed in [12] utilizes three-density version of the multidensity formalism developed earlier [13].In this paper we propose extension of the TPTCF with the possibility of using an arbitrary number of the densities.This enables the theory to be used in describing the network forming fluids.We apply TPTCF to the study of the phase behavior of Stockmayer fluid with the dipole moment µ * d from a wide range of values (1 µ * d 2

36
) and compare our theoretical results against currently available computer simulation results.We also consider the model of the dipole fluid with the pair potential represented by Stockmayer potential with attractive part of the LJ potential multiplied by the variable parameter λ.For λ = 1 this potential is identical to original Stockmayer potential and for λ = 0 it reduces to soft-sphere dipole potential.This model was used to study the effects of the dispersive interaction on the phase behavior of the dipolar fluids [10].The paper is organized as follows.In section 2 we present details of the theory and in section 3 specify the potential model to be studied.In section 4 we discuss our choice of the reference system and its description.Results and discussion are presented in section 5 and in section 6 we collect our conclusions.

Thermodynamic perturbation theory for central force associative potential
As in our earlier studies [13,12] we represent the pairwise total potential energy U (12) and corresponding Mayer function f (12) as a sum of the reference U ref (12), f ref (12) and associative U ass (12), F ass (12) pieces, i.e.: and where F ass (12) = e ref (12) Here z(i) = z exp [−βU (i)], i denotes position and orientation of the particle i, and U (i) is an external field.For a uniform system z(1) ≡ z.
Following [12][13][14] we introduce the definition of the S-mer diagrams.These are the diagrams consisting of S circles, which are connected by F ass bonds.The circles, which are incident with more than m F ass bonds are called oversaturated circles.The set of all possible S-mer diagrams can be constructed in three steps: (i) generating the subset of all possible connected diagrams with only F ass bonds, (ii) inserting combined bond e ref = f ref + 1 between all pairs of circles which are adjacent to the same oversaturated circle and are not connected by a direct F ass bond and (iii) taking all ways of inserting a f ref bond between the pairs of circles which were not connected during the previous two steps.As a result the diagrams, which appear in ln Ξ and ρ(1), can be expressed in terms of the S-mer diagrams: ln Ξ is the sum of all topologically distinct simple connected diagrams consisting of S-mer diagrams with S = 1, . . ., ∞ and f ref bonds between pairs of circles in distinct S-mer diagrams.
The procedure for obtaining the expression for ρ(1) from ln Ξ remains unchanged.The diagrams appearing in the z expansion of the singlet density ρ(1) can be classified with respect to the number of F ass bonds associated with the labeled white z(1) circle.We denote the sum of the diagrams with k (0 k < m) F ass bonds connected to the white z(1) circle as ρ k (1) and the sum of the diagrams with m or more F ass bonds connected to the white z(1) circle as ρ m (1).Thus for ρ(1) we have: The process of switching from the activity to a density expansion proceeds in the same fashion as in [13,12].Analyzing the connectivity of the diagrams in ρ(1) at a white z(1) circle we have where and Here c (0) is the sum of all topologically distinct simple irreducible diagrams consisting of S-mer diagrams with S = 1, . . ., ∞ and f ref bonds between circles in distinct S-mers.All circles in c (0) are black circles carrying the factor σ m−l if the corresponding circle is incident with l m F ass bonds.The circles, which are incident with more than m F ass bonds carry the factor σ 0 .
Using the above introduced density parameters σ k (1) the virial expansion for the pressure P can be written in the following form: where V is the volume of the system.This expression satisfies the regular thermodynamic relation β (∂P/∂µ) = ρ, where ρ = ρ(1)d( 1) and µ is the chemical potential.Corresponding expression for Helmholtz free energy A is where Λ is the thermal wavelength.This expression is derived using the regular thermodynamic expression for Helmholtz free energy A = N µ − P V together with the relation which follows from (4).For the excess Helmholtz free energy ∆A = A − A ref we have where N is the number of the particles in the system, ∆c ref are the corresponding quantities for the reference system.Ordering the virial expansion (11) with respect to the number of associating bonds F ass and neglecting the terms with more than one associating bond we have ref (12) ref (12) and where g (2) ref (12) is the reference system pair distribution function.Equality (14) together with (5) yields the following relation between the densities Set of the relations ( 12), ( 13) and ( 15) define all the quantities needed to calculate Helmholtz free energy of the system (11), provided that the properties of the reference system are known.For the uniform system, free energy expression (11) will be simplified where σ 0 and σ 1 satisfy the following set of equations and ref (12)f ass (12)d( 2).

The model
We consider the fluid with Stockmayer pair potential, which is represented by the sum of the Lennard-Jones (LJ) potential U LJ (r) and dipole-dipole potential U dd (12) where and Here σ LJ and LJ are LJ potential size and energy parameters, respectively, µ d is the dipole moment, θ 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.

Definition of the reference system and its properties
Following [12] we have Here q plays a role of the potential splitting parameter, since at different values of q different portions of the total potential will be included into its reference and associative parts.In the limiting case of q = 0 we have: U ass (12) = 0 and U ref (12) = U (12).Due to our previous study [12] the optimal choice for this parameter is q = 3 and this value is used here.According to the above choice of the reference system potential, the potential energy minima at "nose-to-tail" configuration (θ 1 = θ 2 = 0, θ 1 = θ 2 = π) are included into the associative potential U ass (12).These minima are responsible for the formation of the chains in the system.According to SRN scheme [5,6], Helmholtz free energy of the reference system A ref can be written using the Padé approximation where ref is the LJ contribution to Helmholtz free energy, which was calculated using the equation of state for the LJ fluid proposed by Kolafa and Nezbeda [16], LJ (r)dr [U ref (12)] LJ (r 12 , r 13 , r 23 )dr where g LJ (r) is the LJ fluid radial distribution function (RDF) and dΩ denotes integration over the orientations.After carrying out the angular integration in ( 27) and (28), we have LJ (r)dr , LJ (r 12 , r 13 , r 23 )dr 2 dr 3 , where expressions for w 2 (r) and w 3 (r 12 , r 13 , r 23 ) are presented in the Appendix and for the triplet distribution function g LJ (r 12 , r 13 , r 23 ) we have used superposition approximation, i.e. g LJ (r 12 , r 13 , r 23 ) = g LJ (r 12 )g LJ (r 13 )g LJ (r 23 ) .
Finally for the pair distribution function of the reference system g ref (12), the following approximation was utilized g ref (12) = y LJ (r) exp (−βU ref (12)) , where y LJ (r) is the LJ fluid cavity distribution function.This function was calculated using solution of the Ornstein-Zernike (OZ) equation supplemented by the closure conditions due to Duh et al. [17][18][19].

Results and discussion
Thermodynamical properties of Stockmayer fluid was calculated using three-density version (m = 2) of the TPTCF expression for excess Helmholtz free energy (16) and Padé expression for the reference system Helmholtz free energy (26).RDF g (2) LJ (r), which is needed as an input in ( 16) and (26), was calculated using OZ equation together with the closure relation suggested by Duh et al. [17][18][19].For the pressure P and chemical potential µ the standard thermodynamical relations  Predictions of the present theory for the liquid-gas phase behavior of Stockmayer fluid at different values of the dipolar moment µ d are shown in figure 1.Our results are compared with computer simulation results [8,9,20] and the results of the SRN type of the theory, to which the present approach will be reduced if in (24) and (25) one assumes the zero value for the splitting parameter q, i.e. q = 0.For the low values of the dipole moment (µ * d 2 = 1 ÷ 2) one can see a good agreement between theoretical and computer simulation predictions.With the dipole moment increase, the accuracy of the SRN approach rapidly decreases; at high values of the dipole moment (µ * d 2 = 16 ÷ 36) its predictions become unsatisfactory.In the whole range of the dipole moment values studied our theory gives predictions, which are in reasonable agreement with computer simulation predictions.In general, our results for the gas branch of the coexistence curve are more accurate than those for the liquid branch.Predictions of the present theory for the critical temperature T * cr are in a good agreement with computer simulation predictions.At the same time SRN approach gives slightly better predictions for T * cr at the lowest values of the dipole moment studied (µ * are placed nearby corresponding curves.Solid curves denote results of the present theory, dashed lines denote results of SRN approach and empty circles denote computer simulation results [8,9].Critical points are denoted by filled symbols.
Finally, using the theory developed above for Stockmayer fluid, we study the phase behavior of the model dipolar fluid with the following pair potential where 0 λ 1.For λ = 1 this potential coincides with Stockmayer potential and for λ = 0, expression (34) represents a potential model for soft-sphere dipolar fluid.This model was used to study the role of the dispersive interaction in the phase behavior of the dipolar fluids [10].According to computer simulations studies, carried out in this work, there was no evidence of the liquid-gas phase transition for λ < 0.3.However, recent computer simulation studies, carried out for the dipolar hard-sphere fluid, present the evidence for the liquid-gas phase transition [21] with the critical point located at low temperature and density.Since there is a strong similarity between hard-sphere and present (λ = 0) soft-sphere dipolar models one may expect that liquid-gas phase transition for the latter model exists and occurs at a temperatures lower than those studied in [10].Thermodynamical properties of the model at hand were calculated using the scheme outlined above.Since the introduction of λ modifies LJ part of the interaction in (34) we were not able to use the equation of state due to Kolafa and Nezbeda [16].The term A (0) ref , which appears in the expression for Helmholtz free energy (26), was calculated using perturbation theory of Barker and Henderson [22].Our results for the phase behavior of the present model together with the corresponding results of computer simulation study [10]   the critical density (figures 3 and 4b) and for the liquid branch of the coexistence curve (figure 3).Unlike computer simulation study [10] our study predicts the existence of the liquid-gas phase transition for λ < 0.3.In the limiting case of soft-sphere dipolar model (λ = 0) theoretical predictions for the critical temperature and critical density are T * cr = 0.404 and ρ * = 0.043, respectively.

Conclusions
Three-density version of the thermodynamic perturbation theory for associating fluids with central force associating potential developed earlier [12] is extended and a corresponding version of the multi-density thermodynamic perturbation theory is proposed.This extension enables the theory to be used in describing the network forming fluids.We apply the theory to the study of the phase behavior of Stockmayer fluid with the dipole moment µ * d in a wide range of its values (1 ). Results of the theory are in a reasonable agreement with currently available computer simulation results in the whole range of the dipole moment values studied.The theory is also applied to the study of the phase behavior of Stockmayer fluid with additional parameter, which makes a dispersive part of interaction variable [10].According to our calculations the liquidgas phase equilibrium exists in the whole range of the values of this parameter, including the limiting case of soft-sphere dipolar fluid model.

Appendix
Angular integration in ( 27) and (28) was carried out analytically using Maple computer algebra package.The final result reads: 33) have been utilized.Here Z = βP V /N .In what follows the dipole moment µ d , the pressure P , the temperature T and the density ρ of the system will be represented by the dimensionless quantities µ * d = µ d / LJ σ 3 LJ , P * = P σ LJ / LJ , T * = k B T / LJ and ρ * = ρσ 3 LJ , respectively.

Figure 1 .
Figure 1.Liquid-gas phase diagram of Stockmayer fluid at different values of the dipole moment µ * d in T * vs ρ * coordinate frame.The values of the squared dipole moment µ * d 2are placed nearby the corresponding curves.Solid curves denote the results of the present theory, dashed lines denote the results of SRN approach and empty circles denote computer simulation results[8,9,20].Critical points are denoted by filled symbols.

d 2 = 1 ÷ 2 )Figure 2 .
Figure 2. Vapor pressure along the coexistence curves at different values of the dipole moment µ * d .The values of the squared dipole moment µ * d 2 are placed nearby corresponding curves.Solid curves denote results of the present theory, dashed lines denote results of SRN approach and empty circles denote computer simulation results [8,9].Critical points are denoted by filled symbols.

at µ * d 2 = 4 Figure 3 .
Figure 3. Liquid-gas phase diagram of Stockmayer fluid with variable dispersive interaction (34) for µ * d 2 = 4 and different values of the switching parameter λ.The values of λ are placed nearby the corresponding curves.Solid curves denote the results of the present theory and empty circles denote computer simulation results [10].Critical points are denoted by filled symbols.

Figure 4 .
Figure 4. Critical temperature (a) and critical density (b) vs λ.Solid curves denote results of the present theory and empty circles denote computer simulation results [10].

w 2
(r) = A r 6 , w 3 (r 12 , r 13 , r 23 ) = 3 r 12 r 13 r 23 3 B(θ 1 , θ 2 ), (35)where A = 17825792/225450225 andB (θ 1 , θ 2 ) =85 (12)s(12), e ref(12)= exp [−βU ref(12)], β = 1/k B T , f (12) = e(12) − 1, T is the absolute temperature and k B is the Boltzmann's constant.Here 1 and 2 denote positions and orientations of the particles 1 and 2. Due to the decomposition of the Mayer function f (12) (2) we will have the following diagrammatic expressions for the logarithm of the grand partition function Ξ and for the one-point density ρ(1) in terms of the activity z: ln Ξ is the sum of all topologically distinct simple connected diagrams consisting of black z circles, f ref and F ass bonds.Each bonded pair of z circles has either a f ref or a F ass bond.