# Equation of state of sticky-hard-sphere fluids in the chemical-potential route

###### Abstract

The coupling-parameter method, whereby an extra particle is progressively coupled to the rest of the particles, is applied to the sticky-hard-sphere fluid to obtain its equation of state in the so-called chemical-potential route ( route). As a consistency test, the results for one-dimensional sticky particles are shown to be exact. Results corresponding to the three-dimensional case (Baxter’s model) are derived within the Percus–Yevick approximation by using different prescriptions for the dependence of the interaction potential of the extra particle on the coupling parameter. The critical point and the coexistence curve of the gas-liquid phase transition are obtained in the route and compared with predictions from other thermodynamics routes and from computer simulations. The results show that the route yields a general better description than the virial, energy, compressibility, and zero-separation routes.

###### pacs:

05.70.Ce, 61.20.Gy, 61.20.Ne, 65.20.Jk## I Introduction

The chemical potential of a fluid can be evaluated as the change in the Helmholtz free energy when a new particle is added to the system through a coupling parameter Onsager (1933); Kirkwood (1935); Hill (1956); Reichl (1980). The coupling parameter determines the strength of the interaction of the added particle to the rest of the system and usually varies between zero (no interaction) and unity (full interaction). This method provides the equation of state (EOS) of the fluid in the so-called chemical-potential route (or route). This can be considered as the fourth route in addition to the better known routes based on the pressure (or virial), energy, and compressibility equations Hansen and McDonald (2006). It must be noted that all these ways to obtain the EOS are formally equivalent.

In practice, the various thermodynamic routes have been mostly developed, under the assumption of additive pair interactions, using the so-called radial distribution function (RDF) . Within this class of interactions, the evaluation of thermodynamic properties of a classical fluid reduces to finding the corresponding RDF. Since all well-known theoretical methods to obtain give approximate solutions, with the exception of a few, simple fluid models (for example, one-dimensional systems whose particles interact only with their nearest neighbors Salsburg et al. (1953)), the EOS obtained from different routes differ in general from one another.

The route has been largely unexplored, except in the scaled-particle theory Reiss et al. (1959); Lebowitz et al. (1965); Mandell and Reiss (1975); Heying and Corti (2004); Stillinger et al. (2006). Recently, one of us used this method to obtain a hitherto unknown EOS for the hard-sphere (HS) model in the Percus–Yevick (PY) approximation Santos (2012). This method was then extended to multicomponent fluids for arbitrary dimensionality, interaction potential, and coupling protocol Santos and Rohrmann (2013). Its application to HS mixtures allowed us to provide a new EOS of this classical model in the PY approximation Santos and Rohrmann (2013) and to derive the associated fourth virial coefficient in the hypernetted-chain approximation Beltrán-Heredia and Santos (2014). Evidently, the route represents a helpful tool for the construction of new EOS and the analysis of thermodynamic properties of fluids. It is therefore of great interest to consider its application to non-HS models.

In this paper we use the route to evaluate the EOS of the three-dimensional sticky-hard-sphere (SHS) model introduced by Baxter Baxter (1968). In this fluid, impenetrable particles of diameter interact through a square-well potential of infinite depth and vanishing width. The study of this pair potential model has two advantages. First, it admits an exact analytical solution to the Ornstein–Zernike equation with the PY closure Baxter (1968), its thermodynamic properties being described in terms of two simple parameters, the packing fraction and a stickiness parameter Barboy (1974). Second, the SHS model has proved to provide an excellent starting point for the study of colloidal systems with short-rang attraction Chen et al. (1994); Verduin and Dhont (1995); Rosenbaum et al. (1996); Pontoni et al. (2003); Buzzaccaro et al. (2007), interactions between protein molecules in solution Piazza et al. (1998), and other interesting applications Noro and Frenkel (2000); Maestre et al. (2013).

We will exploit the known exact solution of the PY integral equation for both single and multicomponent SHS fluids Baxter (1968); Perram and Smith (1975) to obtain the EOS through the route and compare the outcome with the three standard routes (virial, energy, and compressibility), with the less known zero-separation (ZS) route Barboy and Tenne (1976, 1979), and with Monte Carlo (MC) simulations Miller and Frenkel (2004). As we will see, the route EOS in the PY approximation changes with the choice of the prescription followed to switch on the extra particle to the rest of the system. Despite this, the spread is typically much smaller than the one existing among the different routes. Interestingly enough, the route generally provides the best results, including the gas-liquid transition properties of the fluid.

The paper is organized as follows. In Sec. II we give the mathematical formulation of the route for SHS fluids. In Secs. III and IV we use the known exact and PY solutions of the SHS system in one and three dimensions, respectively to derive the route EOS of those systems. Section V is reserved for discussions of the results. The relevant calculations are presented in a series of appendixes.

## Ii Chemical-potential route

We consider a -dimensional system of volumen containing spherical particles of diameter with surface adhesion. The SHS interaction potential between two particles with centers separated a distance is defined by

(1) |

where , and being the Boltzmann constant and the absolute temperature, respectively. The dimensionless parameter measures the strength of surface adhesion (stickiness). The equality in (1) must be interpreted as the limit of an increasingly deeper () and narrower () square-well potential of depth and (relative) width with

(2) |

kept constant. The stickiness parameter is related to the Baxter temperature Baxter (1968) by . The pure HS model is recovered from Eq. (1) in the limit .

We now include into the system an additional particle (the solute). Its interaction with any other particle in the fluid (the solvent) is given by

(3) |

where plays the role of a coupling parameter and is a continuous function of encoding the strength of the solute-solvent attractive force. It runs from at to at .

For large , the excess chemical potential , being the contribution of the corresponding ideal fluid, can be written as follows Kirkwood (1935); Hill (1956); Reichl (1980); Santos (2012); Santos and Rohrmann (2013)

(4) |

where

(5) |

is the configurational integral of solvent particles plus one solute particle with a coupling parameter . Here, , where refers to all the translational coordinates of the solvent particles and refers to the coordinates of the solute particle. Furthermore,

(6) |

is the total potential energy, being the distance between particles and . Hence,

(7) |

where denotes the solvent potential energy.

For convenience, we decompose the right-hand side of Eq. (4) into two separate contributions,

(8) |

Henceforth, for the sake of simplicity, we adopt . The two contributions in Eq. (8) are worked out as follows. First, with Eq. (3), we write

(9) |

For , the surfaces defined by () do not overlap, so that the condition implies . As a consequence, the integrals in (7) of order two or higher in vanish. On the other hand, integration of over gives the free volume of the solute particle,

(10) |

where is the volume of a sphere of radius . Furthermore,

(11) |

being the surface of a sphere of radius . Therefore,

(12) |

With this result, Eq. (7) yields

(13) |

where

(14) |

is the configurational integral of the solvent. As shown in Appendix A, the case is singular if . This difficulty can be overcome by the choice . Therefore, taking into account that , and taking the limit in Eq. (13), the first term on the right-hand side of Eq. (8) becomes

(15) |

where

(16) |

is the packing fraction.

For one must follow another strategy because integration over in Eq. (7) depends upon the coordinates of all the solvent particles. In this case, we consider

(17) |

The solute-solvent RDF is expressed as Hill (1956)

(18) |

It follows from Eqs. (7), (II), and (18) that

(19) |

where

(20) |

is the solute-solvent cavity function,

(21) |

and in the last step of Eq. (II) we have used spherical coordinates.

Finally, inserting Eqs. (15) and (II) into Eq. (8), we obtain

(22) |

This gives the excess chemical potential of -dimensional SHS fluids as obtained from the coupling parameter procedure. To have an expression for more explicit than Eq. (21), we note that, according to Eq. (3),

(23) |

Thus,

(24) |

We may now derive the compressibility factor ( being the pressure) in the route. The familiar thermodynamic relation

(25) |

can be expressed as

(26) |

so that

(27) |

Thus, making use of Eq. (22), we obtain

(28) |

This constitutes the EOS of -dimensional SHS obtained from the route (hence the superscript in ). The better known virial, energy, and compressibility routes are worked out in Appendix B.

## Iii Sticky hard rods: Exact results

As a test of the correctness of Eq. (22), we prove in this section that it leads to the exact EOS for the one-dimensional system (). In that case, Eq. (22) reduces to

(29) |

with

(30) |

As shown in Appendix C,

(31) |

Thus, Eq. (30) may be written in the form

(32) |

Then, Eq. (29) becomes

(33) |

This result is exact and does not depend on the explicit form of in the interval . Making use of Eq. (69), it is straightforward to check that Eq. (26) is indeed satisfied.

## Iv Sticky hard spheres: Percus–Yevick theory

The excess chemical potential for three-dimensional SHS fluids () is obtained from Eqs. (22) and (24) as

(34) |

(35) |

The associated route compressibility factor is given by [see Eq. (II)]

(36) |

The evaluation of Eq. (35) requires the contact values of the solute-solvent cavity function and its derivative . These may be obtained using the PY approximation for an SHS binary mixture (see Appendix D). In particular, and are given by Eqs. (90) and (D.2), respectively.

For an explicit evaluation of Eqs. (34)–(IV) we need to specify the -dependence of (within the constraints and ). In this paper, we shall consider results based on three representative prescriptions:

(37) |

These three protocols are depicted in Fig. 1. In all of them, the solute-solvent stickiness monotonically grows from zero to the solvent-solvent value as the solute diameter () grows from zero to the solvent diameter (). At a given solute diameter, the strength of the solute-solvent attraction increases when going from A to C.

Since the PY integral equation is an approximate theory, a common RDF is expected to yield different EOSs depending on the route followed. In the case of the route, as will be seen below, an extra source of thermodynamic inconsistency arises: the EOS depends on the choice for the protocol .

### iv.1 Virial expansion

Exact | |||||||
---|---|---|---|---|---|---|---|

Undetermined | |||||||

ZS | |||||||

Undetermined | ||||||||

ZS | ||||||||

A standard method of examining different approximations in statistical mechanics is to compare the successive terms in the virial expansion of the compressibility factor. For SHS fluids,

(38) |

The virial coefficients in the virial, energy, compressibility, and chemical-potential routes can be respectively evaluated from Eqs. (59), (65), (63), and (IV), complemented by the PY results summarized in Appendix D. The virial coefficients corresponding to an additional route, the so-called ZS route Barboy and Tenne (1976, 1979), can be derived within the PY approximation from Eq. (85).

All the routes in the PY approximation yield the exact second virial coefficient:

(39) |

As for the third virial coefficient, its exact expression,

(40) |

is recovered from the virial, energy, compressibility, and chemical-potential routes, but not from the ZS route. The ZS result is

(41) |

which is especially wrong in the HS limit ().

The exact fourth virial coefficient is a sixth-degree polynomial in Post and Glandt (1986), i.e.,

(42) |

where the numerical coefficients are given by the first row of Table 1. The PY predictions for depend on the thermodynamic route. They have the structure of Eq. (42), except that . The corresponding numerical coefficients are displayed in Table 1. Note that the energy route is unable to fix the HS EOS, so that the coefficients remain undetermined in that route. All the coefficients derived from the route are rational numbers although a limited number of digits is shown in Table 1. Note that the three protocols of the route agree in the values of and . However, the coefficients – depend on the choice of . For an arbitrary function , they are

(43) |

(44) |

(45) |

(46) |

The three protocols in Eq. (37) have the common form with , , and for C, B, and A, respectively. Taking as a free parameter and using Eqs. (43)–(46), it is possible to find the optimal value of that makes for a given value . For instance, the optimal values are , , , , and for , , , , and , respectively. For the mathematical solutions of are , but these are nonphysical values violating the condition .

Figure 2 compares the exact with various PY routes, where the Carnahan–Starling (CS) Hansen and McDonald (2006) value has been taken in the case of the energy route. A very poor behavior of the ZS route is observed. In what concerns the other four routes, small deviations occur among them for low and moderate stickiness (), a good agreement with the exact values being found in that range. For , however, larger discrepancies occur, with the energy and virial routes showing the most extreme deviations with respect to the exact solution. The route predictions lie between the virial and the compressibility ones, becoming closer to the exact values as a softer stickiness prescription is used (protocol A).

Although, to the best of our knowledge, the fifth virial coefficient is not exactly known, it is worthwhile comparing the different PY predictions for it. They have the polynomial structure

(47) |

the coefficients being presented in Table 2. Again, all the coefficients are rational numbers. The dependence of on the stickiness parameter is shown in Fig. 3 (with the CS choice for the energy route). As expected, the influence of the thermodynamic route on is stronger than in the case of . The general shapes of in the and compressibility routes are intermediate between those in the virial and energy routes.

### iv.2 Weakly sticky limit

Undetermined | |||

ZS | |||

As a complement of the virial expansion (38), it is of interest to examine the leading terms in the series expansion

(48) |

of the compressibility factor in powers of the stickiness parameter. Obviously, the zeroth-order coefficient in the expansion is just the compressibility factor of the pure HS system. Equation (48) can be interpreted as a high-temperature expansion.

Making use of the results of Appendix D in Eqs. (59), (65), (63), (85), and (IV), the first-order and second-order coefficients from the different routes can be derived. The results are displayed in Table 3. Interestingly, one has , thus generalizing the results and observed in Tables 1 and 2. This reinforces that the natural extension of the energy route to HS fluids is the virial EOS Santos (2005, 2006).

The HS EOS was already derived in Ref. Santos (2012). We observe that is protocol-independent. This generalizes to any order in density the behavior observed for and in Tables 1 and 2, respectively. The expression of for arbitrary is

(49) |

The coefficients and , relative to the virial-route predictions, obtained from various PY routes (except the ZS one, in order to avoid distortion of the scales) are shown in Fig. 4. The discrepancies grow as density increases. In particular, the largest inconsistencies occur between the energy and compressibility routes. On the other hand, the route deviates only slightly from the virial route, especially in the case of .

### iv.3 Finite density and stickiness

After having examined the low-density and low-stickiness (or high-temperature) regimes, we now consider the full non-perturbative regime. The density dependence of the reduced pressure is plotted in Fig. 5 for different values of and for the three protocols (37). Of course, in the limit of zero stickiness (), the choice of the protocol becomes irrelevant and one recovers the EOS (see Table 3) corresponding to the PY theory in the route Santos (2012). As increases, the pressure decreases with respect to the HS value and the influence of the protocol is practically negligible up to . For higher stickiness, however, the values of are increasingly sensitive to the protocol chosen. We observe that, as expected on physical grounds, the stronger the relative stickiness , the smaller the pressure.

The three higher values of in Fig. 5 correspond to the gas-liquid critical values , , and predicted by the PY approximation in the virial, energy, and compressibility routes, respectively (see below). In fact the kink in at and reflects the fact that is the critical point for the existence of real solutions of the PY equation (see Appendix D).

In Fig. 6 we compare MC simulations at Miller and Frenkel (2004) with PY predictions from the different routes. Since, as discussed before, the energy route leaves the integration constant undetermined, henceforth the CS EOS Hansen and McDonald (2006)

(50) |

will be taken to complete the determination of via Eq. (63), despite the fact that the choice would be more consistent Santos (2005, 2006). The virial, compressibility, and ZS data have been obtained from Eqs. (59), (65), and (85), respectively. In all the cases, use has been made of the PY solution detailed in Appendix D.

We observe that in the low-density range () all PY routes and simulation data agree very well. For higher densities, the ZS pressure grows too rapidly and the curves corresponding to the three different protocols of the route remain rather close in comparison with those from the virial, energy, and compressibility routes, which show a larger spread. In the range , the route gives the best fits to the simulation data. In the same region, and overestimate the simulation values, while underestimates them. Up to , one has . Finally, there is a rather strong disagreement of all the PY routes at high densities, , where the simulation data exhibit lower pressure values than the theoretical ones. Aside from the ZS curve, the compressibility route shows the largest deviations from MC results on the whole range of studied densities.

### iv.4 Gas-liquid transition

MC | ZS | |||||||
---|---|---|---|---|---|---|---|---|

The two highest values of in Fig. 5 show a domain of mechanical instability, (i.e., ), which indicates a phase transition according to the route in the various protocols analyzed. As is well known, the SHS fluid has a (metastable) fluid-fluid transition which was early predicted in the compressibility Baxter (1968) and energy Watts et al. (1971) routes of the PY approximation. Critical values of the parameters and (or and ) can be determined by the conditions and . Results concerning to various PY routes are summarized in Table 4, where they can be compared to the ones obtained by Miller and Frenkel Miller and Frenkel (2004) using MC simulations.

As is known, the compressibility route produces a gross underestimation of the critical density. An even higher underestimation of is obtained from the ZS method. The critical density is much better approximated by the virial route (deviation of ) and, especially, the route (deviations of , , and for protocols A, B, and C, respectively). On the other hand, the critical value of the stickiness parameter evaluated from the virial and the compressibility routes differ significantly from those predicted by numerical experiments. For this parameter, the ZS route (deviation of ), the energy route (deviation of ) and the route (deviations of , , and for A, B, and C, respectively) give the best results. In view of the general poor performance of the ZS route, its good prediction of the critical stickiness can be viewed as accidental. In addition, it must be remarked that the critical point obtained from is quite sensitive to the choice of . If, instead of the CS EOS (50), the more consistent choice Santos (2005, 2006) is used, then no SHS critical point is predicted by the energy route. In conclusion, the parameters of the critical points obtained from the route show the best global agreement with simulations.

We have also computed coexistence curves with the familiar equal-area Maxwell construction that is applicable when . Figure 7 displays the coexistence curve and the location of the critical point derived by various PY routes and from computer simulations Miller and Frenkel (2004). For simplicity, in the case of the route only results from protocol A are shown. As may be seen in Fig. 7, the curves obtained from the virial, compressibility, and ZS routes differ substantially from computer evaluations. On the contrary, the agreement is reasonably good for the energy and routes. As already seen from Table 4, the critical Baxter temperature predicted by the energy route is more accurate than the one obtained from the route. However, the general shape of the coexistence curve (both the gas and the liquid branches) is better described by the route.

## V Conclusions

In this paper we have studied the SHS fluid using the concept of a partially coupled particle, whereby the interaction potential connecting this particle (the solute) to all other particles in the fluid (the solvent) is regulated by a charging parameter which varies from (no interaction) to (full interaction). With this method, first introduced by Onsager Onsager (1933) and subsequently developed by Kirkwood Kirkwood (1935) and other authors Hill (1956); Reiss et al. (1959), we have derived the chemical potential of -dimensional SHS fluids in terms of the contact values of the solute-solvent cavity function and its first derivative [see Eqs. (24) and (22)].

The procedure requires and in the range , where the effective diameter () of the coupled particle varies between zero to the solute diameter . Thus, the explicit evaluation of the EOS of the fluid in the route requires the structural functions of the corresponding binary system in the infinite dilution limit. While this may represent a practical disadvantage with respect to the other standard routes (which only require the RDF of the one-component fluid), it nicely complements them. As is well known, the natural variables of the free energy are the temperature (), the volume (), and the number of particles (). The internal energy, the pressure and the isothermal compressibility are directly related to , , and , respectively. Therefore, the chemical potential, being related to , completes the picture.

The route also requires a prescription for the solute-solvent interaction potential, i.e., the dependence of the solute-solvent stickiness on the coupling parameter must be specified. Here, in order to avoid the effects of trimer particle configurations, we have selected prescriptions with