Equation of state for five-dimensional hyperspheres from the chemical-potential route

Equation of state for five-dimensional hyperspheres from the chemical-potential route

René D. Rohrmann rohr@icate-conicet.gob.ar http://icate-conicet.gob.ar/rohrmann Instituto de Ciencias Astronómicas, de la Tierra y del Espacio (ICATE-CONICET), Avenida España 1512 Sur, 5400 San Juan, Argentina    Andrés Santos andres@unex.es http://www.unex.es/eweb/fisteor/andres/ Departamento de Física and Instituto de Computación Científica Avanzada (ICCAEx), Universidad de Extremadura, Badajoz, E-06071, Spain
September 22, 2019

We use the Percus–Yevick approach in the chemical-potential route to evaluate the equation of state of hard hyperspheres in five dimensions. The evaluation requires the derivation of an analytical expression for the contact value of the pair distribution function between particles of the bulk fluid and a solute particle with arbitrary size. The equation of state is compared with those obtained from the conventional virial and compressibility thermodynamic routes and the associated virial coefficients are computed. The pressure calculated from all routes is exact up to third density order, but it deviates with respect to simulation data as density increases, the compressibility and the chemical-potential routes exhibiting smaller deviations than the virial route. Accurate linear interpolations between the compressibility route and either the chemical-potential route or the virial one are constructed.

61.20.Gy, 05.70.Ce, 61.20.Ne, 65.20.Jk

I Introduction

The thermodynamic and equilibrium statistical-mechanical behavior of liquids and dense gases can be conveniently described by means of distribution functions Hill (1956); Barker and Henderson (1976); Reichl (2009); Hansen and McDonald (2006). The most widely used one is the pair correlation function or radial distribution function (RDF) , which is particularly appropriate for the study of systems of particles interacting with a pairwise additive potential. For such fluids, there exists a set of well-established, rigorous relationships connecting thermodynamic quantities with configuration integrals over RDF. In particular, the compressibility, energy, and pressure (or virial) equations provide well-known routes to the thermodynamic properties of the fluid Hill (1956); Barker and Henderson (1976); Reichl (2009); Hansen and McDonald (2006); Santos (2014).

Thermodynamic properties can also be obtained from another route that connects with the chemical potential through the charging process of a test particle in the fluid Onsager (1933); Kirkwood (1935); Hill (1956). This represents the chemical-potential route ( route), which has remained almost unexplored until recent years, when it was used to obtain a new Percus–Yevick (PY) equation of state (EOS) for the hard-sphere system Santos (2012). Subsequently, this has been formally generalized to arbitrary multicomponent systems and applied to additive hard-sphere (AHS) mixtures Santos and Rohrmann (2013). In addition, the EOS of sticky hard spheres (or Baxter model Baxter (1968)) was derived and the critical point associated with a liquid-gas transition was captured by PY results in this way Rohrmann and Santos (2014).

All these routes (virial, energy, compressibility, and ) are formally exact and thermodynamically equivalent, provided that the exact RDF is employed Santos (2014). Actually, however, only approximate evaluations of are available for systems in dimensions with nontrivial interactions. In general, the RDF is obtained with some accuracy from numerical simulations Hansen and McDonald (2006), by solving integral equation approximations (e.g., the PY Percus and Yevick (1958) and hypernetted chain approaches van Leeuwen et al. (1959); Morita (1960)) or from density functional theories Evans (1979); Rosenfeld (1989). In this context, the route has demonstrated to provide an alternative and useful path for the study of the thermodynamic properties of fluids.

With the aim of extending the use of the route and to gain insight into its properties, in this paper we apply the method to find the EOS of the hard hypersphere fluid at spatial dimension . Besides the intrinsically interesting properties of hard particle systems, they are an important basis for constructing more complicate models, so that there are active theoretical efforts to study them in dimensions (see, for instance, Refs. Colot and Baus (1986); Luban and Michels (1990); Santos et al. (2001); Finken et al. (2002); González-Melchor et al. (2001); Clisby and McCoy (2004); Lyberg (2005); Clisby and McCoy (2006); Torquato and Stillinger (2006a, b); Rohrmann and Santos (2007); Adda-Bedia et al. (2008); van Meel et al. (2009); Lue et al. (2010); Leithall and Schmidt (2011); Piasecki et al. (2011); Estrada and Robles (2011); Zhang and Pettitt (2014) and references therein).

The route is based on a charging process where a test particle (solute) with tunable interaction is inserted into the bulk fluid (solvent particles). As a consequence, the application of the method requires the solute-solvent RDF of the corresponding binary system in the infinite dilution limit. To our knowledge, the only systematic theoretical method for the evaluation of the RDF of AHS fluid mixtures at dimensions higher than is the so-called rational-function approximation (RFA) Rohrmann and Santos (2011a). Its simplest implementation provides the solution of the Ornstein–Zernike relation coupled with the PY closure for AHS mixtures in odd dimensionalities Rohrmann and Santos (2011b). Therefore, we adopt the RFA technique at the PY level in the present study.

This work is organized as follows. In Sec. II we present the basic formulation of thermodynamics routes for hypersphere systems. Section III gives, within the PY theory, the five-dimensional RDF for a coupled particle with arbitrary size, a key quantity for evaluating the EOS in the route. The technical details are elaborated in the Appendixes. The results are presented in Sec. IV. Finally, in Sec. V we offer some conclusions.

Ii Framework

Figure 1: (Color online) Cartoon of the charging process of the solute from a point particle () to a particle equivalent to any solvent particle ().

The EOS for a single component system of hard -sphere particles can be expressed as a relationship between the compressibility factor (where is the pressure, is the average particle density, is the Boltzmann constant, and is the temperature) and the packing fraction , with the diameter of the particles and the volume of a -dimensional sphere of unit diameter. In the route, the EOS is given by Santos (2012); Santos and Rohrmann (2013); Santos (2014); Beltrán-Heredia and Santos (2014)




where is the contact value of the solute-solvent RDF, which depends on both and the coupling parameter . The latter is defined as the ratio between the solute and the solvent diameters, so that the minimum possible distance between solute and solvent particles is


Thus, regulates the strength of the interaction between the test particle and the rest of the fluid. When the test particle is a point that cannot penetrate the solvent particles, while when the test particle is indistinguishable from any particle of the bulk fluid. This charging process is schematically illustrated by Fig. 1.

For comparison, the compressibility factor in the virial and compressibility routes are expressed by Rohrmann and Santos (2007)




where is the isothermal susceptibility. Since hard spheres are athermal, the energy route to the EOS becomes useless.

In order to evaluate within the RFA approach Rohrmann and Santos (2011a), it is convenient to introduce the Laplace functional defined by


where is the RDF of the pair () and is the reverse Bessel polynomial of order Rohrmann and Santos (2007). This functional is directly related to the static structure factors of a multicomponent fluid,


where is the wavenumber, is the mole fraction of species , and i is the imaginary unit. The functional provides us with all the necessary information about the structure and thermodynamics of the fluid state. In particular, the contact values are Rohrmann and Santos (2011a)


being the contact distance of the pair ().

Iii PY approach

The exact solution of the PY approximation for -odd AHS mixtures with any number of components has been described in detail in Ref. Rohrmann and Santos (2011a) within the framework of the RFA method. For the sake of consistency, the main expressions particularized to Five-dimensional () binary mixtures are presented in Appendix A. From them one can take the limit of solute infinite dilution to obtain the solute-solvent contact value . The details are worked out in Appendix B and here we only quote the final result:




In the special case one recovers the solvent-solvent value, namely.


In the opposite limit one obtains the exact results Santos (2012) and . Figure 2 shows the -dependence of for the representative packing fraction .

Figure 2: (Color online) Plot of the solute-solvent contact value versus the solute-to-solvent size ratio at a packing fraction . The values of and its slope at are exact, but the full curve is the approximate prediction of the PY theory [see Eq. (9)].

Once has been analytically determined, the integral in Eq. (2) can be evaluated with the explicit result


where . The compressibility factor in the route can be easily evaluated from Eq. (1) by numerical integration.

Equation (11) can be used to evaluate the virial route to the EOS from Eq. (4). Finally, the isothermal susceptibility of hyperspheres in the PY approach is Freasier and Isbister (1981); Rohrmann and Santos (2007)


From here one can obtain the compressibility route to the EOS via Eq. (5) as


Iv Results

iv.1 Virial coefficients

We start this section by considering the virial expansion


for the fluid as predicted by the three routes in the PY approximation. The virial coefficients and corresponding to the virial and compressibility routes are obtained from Eqs. (11) [combined with Eq. (4)] and (14), respectively. As for the route, Eq. (1) implies that


where the coefficients , which are defined by the expansion


are easily obtained from Eq. (12).

Table 1: First virial coefficients obtained from the PY theory according to the virial (), compressibility (), and () routes.
Table 2: Comparison between the numerical values of the virial coefficients for as obtained from several routes in the PY approximation and the exact values () for Lyberg (2005) and Zhang and Pettitt (2014).
Figure 3: (Color online) Normalized differences between the virial coefficients in the virial (squares), compressibility (circles) and (triangles) routes, with respect to the exact values Lyberg (2005); Zhang and Pettitt (2014).

Table 1 contains the first virial coefficients obtained from the three routes in the PY theory. As is well known, the PY theory yields the exact virial coefficients up to third order, regardless of the thermodynamical route. However, discrepancies among results from different routes appear for upper order coefficients (). As a peculiar feature, it may be observed that, in contrast to the three-dimensional (3D) case Santos (2012), the route yields irrational virial coefficients (due to the logarithmic term) for . On the other hand, the virial and compressibility routes yield rational numbers for the virial coefficients of hypersphere systems in all odd dimensions Rohrmann et al. (2008).

Table 2 shows the three sets of PY virial coefficients in numerical format and compares them with the exact result for Lyberg (2005) and with recent accurate values for Zhang and Pettitt (2014). The relative deviations are displayed in Fig. 3. It is quite clear that the compressibility route gives the closest agreement with the exact values, while the values calculated from the and virial routes are rather similar, with a slight improvement of over . It is interesting to remark that the three PY routes capture the alternating sign change between and , while a negative sign of is wrongly anticipated by the virial and routes. In the PY case, the alternating character is related to the existence of a branch point singularity on the negative real axis (at ), which determines the radius of convergence of the virial series Rohrmann et al. (2008).

Regarding the closeness between and , it turns out to be higher in five dimensions than in three dimensions. While in the 3D case the ratio monotonically decreases from () to () Santos (2012), in the 5D case the ratio first increases from () to a maximum value () and then (from ) monotonically decreases towards an asymptotic value . It can then be speculated that the general similarity between and tends to increase with increasing dimensionality. We will return to this point at the end of Sec. IV.2.

iv.2 Equation of state at finite density

As said before, the compressibility factor of hyperspheres in the PY route can be obtained through numerical integration from Eq. (1) by making use of the explicit expression (12). The virial and compressibility routes yield Eqs. (4), combined with Eq. (11), and (14), respectively.

Figure 4: (Color online) Compressibility factors evaluated from computer simulations Lue et al. (2010) (symbols) and from an accurate Padé approximant (dash-dotted line) Bishop and Whitlock (2005); Lue et al. (2010) are compared with PY results in the virial (long dashed line), compressibility (short dashed line), and (solid line) routes. The dotted vertical lines mark the freezing and melting packing fractions van Meel et al. (2009).

Figure 4 shows the density dependence of according to the PY theory in the virial (), compressibility (), and () routes. Predictions of computer simulations Lue et al. (2010) are denoted by symbols, while the dash-dotted line represents a Padé approximant Bishop and Whitlock (2005); Lue et al. (2010) based on high precision calculations of the first ten virial coefficients Clisby and McCoy (2006). The freezing (, ) and melting (, ) densities calculated for the lattice van Meel et al. (2009) are also indicated. As expected, all three PY routes converge to the exact results at very low density. On the other hand, as the packing is increased, the compressibility route overestimates , whereas the virial and routes underestimate it, the values from the route being slightly better than those from the virial route.

Figure 5: (Color online) Relative differences between the compressibility factor evaluated from computer simulations () Lue et al. (2010) or from the Padé approximant () Bishop and Whitlock (2005); Lue et al. (2010) and the PY results in several routes (filled symbols). Open symbols correspond to the interpolations (squares) and (triangles) from Eqs. (18) and (19), with and , respectively.

The relative deviations of , , and from the exact compressibility factor are plotted in Fig. 5. It can be observed that the deviations of the compressibility and routes from computer simulations are very similar (but with different sign) within the fluid phase (). For instance, at , the deviations between simulation or Padé-approximant results and the PY theory become %, %, and % for the virial, compressibility, and routes, respectively, while they increase to %, %, and %, respectively, when the packing fraction is reached.

Figure 6: (Color online) Ratio with (squares) and (triangles), using the compressibility factor obtained either by simulations Lue et al. (2010) (filled symbols) or from the Padé approximant Bishop and Whitlock (2005); Lue et al. (2010). As reference, horizontal lines denote the values , , and , and vertical lines refer to the freezing and melting packing fractions van Meel et al. (2009).

Looking for a better agreement with numerical experiments from PY solutions, one can construct new EOSs through Carnahan–Starling-like interpolations of the form


To test this possibility, one can define density-dependent weight functions and to check to which extent they are weakly dependent on . This is done in Fig. 6, which shows evaluations of the parameter from available data. Variations observed at the lowest densities () are due to the ratio of very small values; actually, all three routes yield accurate results for very dilute gas conditions. On the other hand, it is clear from Fig. 6 than there is not a unique value of or in the fluid phase region (). However, it is to be noted that and may be reasonable choices for and , respectively. It can be observed from Fig. 6 that the value in Eq. (18) is better in the region than the simpler value proposed in Ref. Santos (2000). Data obtained from the hybrid compressibility factors (with ) and (with ) are also included in Fig. 5. Both EOSs have discrepancies with respect to simulation results lower than % for .

The improvement of the route over the virial route observed in the virial coefficients (see Table 2 and Fig. 3) as well as in the EOS at finite densities (see Figs. 4 and 5) can be rationalized by the following heuristic argument Santos (2012). Comparison between the exact statistical-mechanical formulas (1) and (4) shows that, while the route is related to the (weighted) average of the contact value in the range , the virial route is directly related to the local value . Since both and its first derivative are given exactly by the PY approximation at (see Fig. 2), it seems reasonable that the average value is better estimated than the end point by the PY approximation. In this respect, it is also worth noting that the weight function in Eq. (2) increases with (except in the one-dimensional case, ), so that the average value is more influenced by the approximate values of near than by the quasiexact values near . This bias strongly increases with increasing dimensionality [in fact, in the limit ], which explains why the discrepancy between the PY and virial routes is much smaller with than with .

V Final remarks

In this paper, we have considered an application of the Kirkwood coupling parameter method (or route) to the thermodynamics of classical fluids. Specifically, the EOS of the -sphere fluid has been derived from the route by exploiting the knowledge of the exact solution of the PY theory for this system. The basic quantity required for this implementation is the contact value of the RDF for a partially coupled particle, which corresponds to a binary mixture with one infinitely dilute component. Although the PY approximation had been solved for the -sphere system by Freasier and Isbister Freasier and Isbister (1981) and by Leutheusser Leutheusser (1984), their methods do not provide the RDF of a binary mixture required here. For this purpose we have used the RFA method Rohrmann and Santos (2011a) which exactly solves the PY equation for hypersphere mixtures in odd dimensionalities. The RDF obtained in this way has an analytical expression and the resulting EOS in the route is derived from a numerical integration over density.

By examining the first few coefficients in the density expansion for the EOS and comparing them with available exact or highly accurate values, it is seen that, although the PY theory is exact up to the third coefficient, it rapidly worsens for increasing density order. Specifically, the compressibility equation gives in general the most satisfactory values, followed by the route, which becomes somewhat better than the virial one. Contrary to the case of -sphere systems Santos (2012), linear interpolations of the PY EOS from different routes do not improve the evaluation of the virial coefficients.

Regarding the behavior at finite densities, the pressure predicted by the route is observed to lie below to simulation data within the fluid phase. Interestingly, the deviations of the compressibility and routes with respect to the numerical experiments are roughly of the same magnitude (but of opposite signs) in this region and both are lower than those from the virial route. An accurate fit to simulation data in the fluid regime () is obtained if the compressibility and routes are added together with equal weights. An alternative accurate fit to numerical experiments may be derived by a similar interpolation between the virial and compressibility routes with appropriate weights. These results show a rather regular behavior of the three PY routes around the exact EOS within the whole fluid region.

The work of R.D.R. has been supported by the Consejo Nacional de Investigaciones Científicas y Técnicas (CONICET, Argentina) through Grant No. PIP 114-201101-00208. A.S. acknowledges support from the Spanish government through Grant No. FIS2013-42840-P and by the Regional Government of Extremadura (Spain) through Grant No. GR15104 (partially financed by ERDF funds).

Appendix A RDF for a general five-dimensional binary mixture

For an AHS binary mixture at , the exact PY solution of the functional (6) may be expressed by Rohrmann and Santos (2011a)


where () and are matrices. More specifically,


where is the unit matrix and () are diagonal matrices with elements






Here, is the th diameter moment and is the total packing fraction. As for the matrix , it obeys the quadratic equation




In general, the solution to Eq. (27) must be obtained numerically Rohrmann and Santos (2011a). However, an analytical form is possible for a binary mixture with an infinitely diluted component (the solute), as shown in Appendix B.

The asymptotic behavior of the functional required in the evaluation of the RDF contact value [Eq. (8)] can be derived from (20) as


where we have called . It was shown in Ref. Rohrmann and Santos (2011a) that




Appendix B Test particle limit

The expressions in Appendix A apply for any five-dimensional AHS binary mixture in the PY approximation. Now, without loss of generality, we adopt the convention and for solvent and solute particles, respectively, and denote by and the respective diameters. Next, we take the test particle limit for the solute, i.e., , , so that . This drastically simplifies the matrices (28)–(A).

Let us write the matrix as


The element of Eq. (27) yields a closed quadratic equation for . The physical solution is identified by the condition with the result


where is defined by Eq. (10). Since is a coefficient related to the bulk fluid, it does not depend on and its expression is equivalent to that of the one-component fluid [as given by Eq. (E10) of Ref. Rohrmann and Santos (2007)]. As for the remaining elements of the matrix , they are given from Eq. (27) by




It can be easily checked that if the solute is equivalent to a solvent particle (i.e., ), then the elements , , and coincide with , as expected.

Now we determine the contact values. Taking into account that and using Eqs. (24), (26), (39), and (40), one obtains


Insertion of this into Eq. (38) yields


Its inverse matrix is


Finally, from (37) one has