Virial expansions and augmented van der Waals approach: Application to Lennard-Jones-like Yukawa fluid

We argue that recently proposed [Melnyk et al., Fluid Phase Equilibr., 2009, Vol. 279, 1] a criterion to split the pair interaction energy into two parts, one of which is forced to be responsible the most accurate as possible for excluded volume energy in the system, results in expressions for the virial coefficients that improve the performance of the virial equation of state in general, and at subcritical temperatures, in particular. As an example, application to the Lennard-Jones-like hard-core attractive Yukawa fluid is discussed.


Introduction
Despite the significant progress in the development of modern tools in the statistical theory of liquids [1][2][3][4][5][6][7][8][9][10], there still are numerous studies where thermodynamic properties are expressed in the terms of virial expansions (e.g., see references [11][12][13] and references therein). The classical example of the virial expansion approach is the virial equation of state (EOS) [14,15] p k B T = ρ + B 2 (T )ρ 2 + B 3 (T )ρ 3 + B 4 (T )ρ 4 + . . . , (1.1) where p is the pressure, k B is the Boltzmann constant, T is the temperature, ρ is the number density and B n (T ), n = 2, 3, 4, . . . are the virial coefficients. The success of the virial expansion approach relies, first of all, on the knowledge of the virial coefficients. The first two virial coefficients, B 2 (T ) and B 3 (T ), by various techniques can be obtained experimentally, while in theoretical studies both can be relatively easily evaluated numerically (for some fluid models even analytically [16]) following their definition in terms of Mayer cluster integrals [15] B 2 (T ) = − 1 2V f (r 12 )dr 1  is the Mayer function and u(r ) is the pair interaction energy in a target model fluid. The expressions for higher order virial coefficients are much more complicated, especially due to a significant increase in the number of distinct integrals that are required to be evaluated [17]. Therefore, there are not so many potential functions u(r ), for which the virial coefficients B n (T ) of the order n > 2 are known.
Traditionally, the virial expansion approach is the most advanced for the model fluid composed of purely repulsive hard spheres of diameter σ. In this case, the virial coefficients B n are independent of temperature and have been calculated up to the twelfth order [18][19][20][21]. Being explored in the virial EOS, these coefficients lead to the pressure p hs of the hard-sphere (hs) fluid that is rather accurate in comparison with computer simulations data for densities up to the fluid-solid transition [22]. However, the success of the virial expansion approach is not so evident when apart from the hardsphere repulsion, the interaction potential u(r ) and, consequently, the exponential of the Mayer function in equations (1.2)-(1.4) both include the attractive interaction energy between molecules. In this case, the virial expansions approach tends to diverge when approaching the thermodynamic states associated with condensation. This fact imposes serious limitations on the applicability of the virial expansion approach to properly describe the vapour-liquid equilibrium in fluid systems. The problem of virial expansion divergence, in the region of condensation, turned out to be a long-standing issue in the case of Lennard-Jones (LJ) fluid (e.g., see recent papers by Ushcats [23][24][25][26] and references therein).
In this communication we wish to focus on another popular model system in the liquid state theory, namely, the LJ-like hard-core attractive Yukawa (HCAY) fluid model [6] figure 1 for details regarding the relation between the LJ-like HCAY interaction potential and the original LJ potential]. This model fluid has been studied intensively by computer experiment [7,27,28] as well as by other approaches, such as the mean-spherical integral equation theory (MSA) [6], the MSA-based first-order perturbation theory (FMSA) [8], the MSA-based high temperature expansions theory [9,10]. As for the virial expansion approach, to the best of our knowledge, there is only one paper by Naresh and Singh [13], where the virial coefficients up to the sixth order, i.e., B 2 (T ), B 3 (T ), . . ., B 6 (T ), for the LJlike HCAY fluid have been reported. After substituting these coefficients into the virial EOS, equation (1.1), the applicability of the latter in the case of the LJ-like HCAY fluid (see figure 6 for illustration) can be summarized as follows [13]: (i) being truncated by B 6 (T ), the virial EOS is rather accurate in the density range up to a reduced density ρσ 3 ≈ 0.5, and remains to be qualitatively correct for the entire density range, but only for the reduced temperatures T * = k B T / = 1.5 and 2 that are supercritical temperatures for this fluid model (the critical point temperature in this case T * c ≈ 1.2 [7,27,29]); (ii) the same virial EOS begins to fail already right after the density ρσ 3 ≈ 0.1 in the case of the subcritical temperatures (the results shown by the dashed line in figure 6 for the reduced temperature T * = 1).
By using the LJ-like HCAY model fluid as a pilot system, the purpose of this exploratory study is to show that even in the case of the truncated virial EOS, the performance of the latter in the wide range of density and temperature conditions, including the subcritical ones, can be substantially improved by implementing the ideas that were elaborated within the framework of the augmented van der Waals theory [29,30]. These ideas concern the issue of a split of the total interaction potential u(r ) into two terms. Namely, in contrast to presumably van der Waals's suggestion that the total potential energy is composed of the repulsion and attraction contributions, the "augmented" version of the van der Waals theory means that one term is representing the most accurate possible the full excluded volume energy in the system, that is the interaction energy between the neighbouring molecules, while the remaining part is responsible for the weak long-range attractive interaction energy, or the energy of cohesion, between the next-neighbouring molecules.
The remainder of this paper is organized as follows: in section 2, we provide an overview of the augmented van der Waals theory [29,30]. In section 3, we discuss how the ideas of this theory can be implemented within the virial expansion approach and we present the corresponding results for the LJ-like HCAY fluid in section 4. We conclude with section 5.

Augmented van der Waals theory
Recently, it has been shown that the thermodynamics as well as the vapour-liquid equilibrium in the LJ-like HCAY fluid can be rather accurately described within the augmented van der Waals theory [29,30]. In particular, within this theory, the EOS of the LJ-like HCAY model fluid reads, p = p 0 − ρ 2 a, (2.1) where, similar to the original van der Waals suggestion, coefficient a is related to the contribution from the attractive interaction energy between molecules, while pressure p 0 stands for the pressure due to the excluded volume energy. Following van der Waals, the excluded volume pressure p 0 originates from the fact that in the system of a volume V and composed of N molecules with a hard-core diameter σ, each molecule excludes an amount of volume v 0 from being allowed to explore by all other molecules of the system. Thus, the volume accessible for molecules is reduced to V − N v 0 . This phenomenon, that firstly 13501-3 was pointed out by van der Waals, results in the pressure, that is referred to as the excluded volume pressure. The excluded volume itself is uniquely defined upon the distance between each pair of the neighbouring molecules in the fluid (see figure 2). As the first approximation, van der Waals assumed that the excluded volume per molecule is a constant and equals fourfold the molecular volume, i.e., v 0 = b ≡ (2/3)πσ 3 . Indeed, it is the case for a dilute gaseous phase [see figure 2 (b)] when the mean distance 〈r 〉 between molecules is large (more precisely, when the mean distance 〈r 〉 between the centers of the pair of the neighbouring molecules is larger than 2σ). In the dense gaseous phase and, especially, in the liquid phase, the mean distance 〈r 〉 between the neighbouring molecules becomes shorter than 2σ, and, consequently, the excluded volume shells start to overlap [see figure 2 (c)], resulting in the excluded volume per molecule being smaller than fourfold the molecular volume, i.e., v 0 < b. Since it is rather evident that the mean intermolecular distance 〈r 〉 between the neighbouring molecules is affected by the number density ρ, it became necessary to incorporate the density dependence into the excluded volume v 0 as well. The most natural way to comply with this requirement was to utilize the hard-sphere fluid model for the evaluation of the excluded volume contribution to the whole spectrum of fluid properties and to the EOS, in particular. Such an assumption lies behind the perturbation theory of fluids due to Zwanzig [2] and was successfully exploited by Widom [3], Barker and Henderson [4], Weeks, Chandler and Andersen [5] and many others [1].
However, although it is less evident, but the mean intermolecular distance 〈r 〉 between two neighbouring molecules is affected by the temperature as well. Such a feature of the intermolecular distance is mediated by the energy of the short-range attractive interaction energy that together with the energy of the hard-sphere repulsion are present for the pair of neighbouring molecules. This observation suggests that the short-range repulsive and short-range attractive interaction energies between the neighbouring molecules should be incorporated into the scheme to evaluate the contribution of the excluded volume to the pressure.
As the first step to comply with this idea, let us follow the suggestion [29,30] and present the pair interaction energy u(r ) in the form, which is in contrast to the common practice [2][3][4][5] that prefers to utilize another form, u(r ) = u hs (r ) + u attr (r ) , (2.4) which assumes that the pair interaction energy u(r ) is separated into purely repulsive hard-sphere term u hs (r ) and attractive interaction energy u attr (r ) contributions. The interaction energy u nn (r ) in equation (2.3) represents the full interaction energy of a target molecule and its neighbouring (nn) counterpart. The neighbouring molecules and the corresponding interaction energy u nn (r ) are identified by means of the range (distance) criterion. According to this criterion, the excluded volume interaction energy u nn (r ) includes the full energy of the hard-core repulsion u hs (r ) and only a part of the full attrac-tion energy, namely, the part that is responsible for the interaction of a target molecule with its nearestneighbour molecule only, In fact, this is the interaction energy with the molecules that belong to the first coordination shell of a target molecule. Following such a definition, the attraction u sr attr (r ) incorporates the full attraction energy u attr (r ) at the contact distance r = σ between two molecules, but decays faster than the full attraction energy, in order not to exceed the radii of the first coordination shell. In reality, the range of the short-range attraction u sr attr (r ) is around one molecular hard-core diameter σ, and in the present case can be approximated by fixing the decay parameter at z 0 σ = 4. Then, the term u lr attr (r ) in equation (2.3) is determined as the difference, u(r ) − u nn (r ), and reads u lr The pair potential u lr attr (r ) corresponds to the interaction energy of the target molecule with any other molecule but from outside the first coordination shell. Figure 1 shows the total pair interaction energy u(r ), the excluded volume interaction energy u nn (r ), and the long-range attractive interaction energy u lr attr (r ), all according to their definitions by equations (1.5), (2.5) and (2.6), respectively. The shape of the long-range attractive interaction energy u lr attr (r ) [see figure 1 (c)] appears to be crucial [30] for the evaluation of the van der Waals coefficient a, that in general case is given by The function g 0 (r ) in this equation stands now for the radial distribution function of the system with the excluded volume interaction potential u nn (r ). We wish to stress, that only by using for excluded volume interaction energy u nn (r ) its definition according to equation (2.5), it is possible to justify the so-called mean-field assumption, g 0 (r ) = 1, in equation (2.7). In the case of the LJ-like HCAY fluid, this results in the simple expression, we have used for the radial distribution function g 0 (r ) the closed-form analytical equation [31,32] in the case of hard-sphere model, and the Monte Carlo simulation data [29] in the case of short-range attractive Yukawa model.
Obviously, the magnitude of the coefficient a is different for each model. However, the most intriguing insight from figure 3 comes from analysing the values of coefficient a within the same model but obtained from two different equations, equations (2.7) and (2.8), respectively. We can see that in the case of the hard-sphere model for excluded volume interaction [figure 3 (a)], the mean-field and exact values of coefficient a are quite different, both quantitatively and qualitatively. By contrast, in the case of short-range attractive Yukawa model, we could admit the tendency for coefficient a to be the same, independently of equations, (2.7) or (2.8), and, consequently, which is most important -to be independent of the density.

Augmented virial EOS
Augmented virial EOS, which is one of the goals of the present study, will be obtained as the series in powers of the density ρ of the augmented EOS (2.1). Since coefficient a does not depend on density within the augmented van der Waals theory, the series will concern exclusively the excluded volume pressure p 0 . Namely, where the virial coefficients B nn n (T ) are defined exactly as it is discussed in equations (1.2)-(1.4), but instead of the total interaction energy u(r ), the short-range excluded volume interaction energy u nn (r ) must be used in the exponential of the Mayer function.
The resulting augmented virial EOS is obtained by substituting the virial series for the excluded volume pressure, equation (3.1), into the augmented van der Waals EOS, equation (2.1), From the first glance at equation (3.2), one immediately notices two important features concerning the role that the long-range attraction energy u lr attr (r ) plays within the augmented virial expansion approach. First of all, since coefficient a does not depend on the density, it follows that the long-range attraction energy u lr attr (r ) contributes to the second virial coefficient only,  In what follows, we apply the augmented virial EOS, equation (3.2), to calculate the compressibility factor pV /(N k B T ) of the LJ-like HCAY fluid.

Results and discussions
Up to date, only the first five virial coefficients, B nn 2 (T ), . . ., B nn 6 (T ), for the excluded volume interaction energy u nn (r ) are known [13]. Figure 4 shows a set of data for the excluded volume pressure, p 0 ≡ p nn , that result from the virial EOS, equation (3.1), (solid lines) truncated at the sixth virial coefficient as well as those that were obtained from computer experiment (symbols) to compare. There are three isotherms, namely, T * = 1, 1.5 and 2 that correspond to the excluded volume pressure within the short-range attraction Yukawa model, while the forth isotherm represents the pressure of the hardsphere fluid, i.e. corresponds to the case when T * → ∞. The most important conclusions that follow from the results presented in figure 4 concern the accuracy and, perhaps, even more generally -applicability of the virial expansions approach in the case of excluded volume interactions. First of all, we can see that virial EOS, equation  interaction energy, u(r ), is separated into two parts following a common practice [2][3][4][5], i.e., it consists of purely repulsive hard-sphere energy, u hs (r ), and full attractive interaction energy, u attr (r ).
data from computer experiment in the density range 0 < ρσ 3 < 0.6. We note, that this observation prac-  figure 5 (a). In particular, the Boyle temperature [the temperature at which B 2 (T ) assumes zero value] of the LJ-like HCAY fluid in both cases is the same, being fixed at approximately T * B ∼ 2.7. We note, that this feature of the second virial coefficient is sensitive to the way how the total pair interaction u(r ) is split into the excluded volume and long-range contributions. To illustrate this point, figure 5 (a) shows the results for the augmented B 2 (T ) in the case when the u(r ) is split in accordance to the common practice [2][3][4][5] given by equations (2.4), i.e., when the nearest-neighbour interaction potential u nn (r ) consists of the hard-core repulsion u hs (r ) only. We can see, that second virial coefficient in this case is quite different from a rigorous one.
The results for compressibility factor, pV / (N k B T ), of the LJ-like HCAY fluid, that follow from the augmented virial EOS (3.2), are shown in figure 6 (the thick solid lines) to be compared against computer experiment data [7] as well as against the rigorous virial EOS (1.1). We can see that for all three temperatures that include both the supercritical and subcritical conditions, there are notable improvements in 2), truncated at the sixth virial coefficient, while symbols correspond to the computer experiment data by Shukla [7]. The dashed lines show the results of the truncated conventional virial EOS, equation (1.1), that are taken from the study by Naresh and Singh [13]. The temperature conditions are specified in the figure. We note that critical point temperature for the LJ-like HCAY fluid is estimated to be around T * c ≈ 1.2 [7,27,29]. The curves have been shifted for clarity.
searching for an agreement with computer experiment data. However, the most valuable result concerns the performance of the augmented virial EOS (3.2) at the subcritical temperature T * = 1, where the original virial EOS equations (1.1) fails. Some discrepancies between the augmented virial EOS and computer experiment data, that still are observed for densities ρσ 3 > 0.6, are pretty similar to those that we already discussed in figure 4, and they should be attributed to the truncation of the excluded volume virial EOS (3.1) at the sixth virial coefficient.

Conclusions
In the present study, the issue of the performance of the virial expansion approach in the liquid state theory is discussed by using as an example the LJ-like HCAY fluid model. More precisely, we were devoted here to discuss the issue of the divergence of the conventional virial EOS (1.1) at subcritical conditions that for the LJ-like HCAY fluid recently was reported by Naresh and Singh [13].
To deal with this issue, the recent advances [29,30] in the augmented version of the van der Waals theory have been explored. The essence of the van der Waals theory lies in the non-trivial split of total interaction energy into two parts which are forced to be responsible for the excluded volume energy and the cohesive energy, respectively. Traditionally, in what is called "van der Waals picture of liquids" [1][2][3][4][5] it is assumed that excluded volume part is well represented by the hard-core repulsion energy, while the cohesive part is associated with the remaining full attractive interaction energy. By contrast, within the augmented van der Waals theory [29,30], it considers that not only the repulsive part, but the

13501-9
full energy of interaction between the pair of neighbouring molecules, must be treated as the excluded volume energy, and what remains is representing the cohesion energy.
By applying the virial expansion approach to the augmented van der Waals EOS (2.1), it is obtained that the cohesion energy, which is the long-range part of the total interaction energy, contributes to the second virial coefficient, B 2 (T ), only. All other virial coefficients result from the short-range excluded volume interaction energy, u nn (r ), which, however, consists of both the hard-core repulsion energy and short-range attraction energy between the neighbouring molecules.
To define the excluded volume interaction energy, u nn (r ), we have used in this study the distance criterion, requiring that the range of excluded volume interaction should not exceed the molecular hardcore diameter σ. However, from the analysis of the second virial coefficient in figure 5 (a) we may conclude that as the criterion for how much of attraction energy in the total pair interaction energy, u(r ), should be included into the excluded volume term, u nn (r ), might be the requirement of the equality between the augmented second virial coefficient, equations (3.3), and the rigorous one, equation (1.2). Namely, the Boyle temperature that follows from the augmented van der Waals theory should be pretty close to the true one, which follows from the rigorous second virial coefficient; otherwise the physics of two systems could differ as well.
The range of attractive interaction energy seems to be an important issue for the convergence of the virial expansion approach. For example, in the limiting case of the hard-sphere repulsion, u hs (r ), when there is no attractive tail at all, the virial EOS shows no sign of divergence (e.g., see reference [12] and discussion therein). Very similar conclusions can be drawn for the virial EOS (3.1) truncated at the sixth virial coefficient in the case of model fluid defined by the interaction potential u nn (r ). Although this has not been proved in the present study, we still are suggesting that virial expansions for the excluded volume pressure, equation (3.1), do not diverge in the range of temperatures that are of interest for the parent fluid, i.e., the LJ-like HCAY fluid in this case, including temperatures that are below its critical point temperature but higher than its triple point temperature. Some discrepancies between the virial EOS and computer experiment data that are observed in figure 4 for densities ρσ 3 > 0.6, should be attributed to the truncation of the virial EOS (3.1) at the sixth virial coefficient. The resulting augmented virial EOS, equation (3.2), has been tested for the LJ-like HCAY fluid in the wide density range and for temperature conditions that were studied in the literature so far, including those where the conventional virial EOS, equation (1.1), exhibits difficulties [13]. Being rather accurate in the range of densities up to around ρσ 3 = 0.6 at the supercritical temperatures, the augmented virial EOS, equation (3.2), remains qualitatively correct at subcritical temperatures as well, showing no sign for divergence at the temperature as low as T * = 1. Nevertheless, for making final conclusion regarding the performance of the augmented virial EOS in the case of the LJ-like HCAY fluid, as well as for the potential application of this approach to investigate more complex and/or realistic class of fluid models, the evaluation of the higher order virial coefficients, namely, B nn 7 (T ), B nn 8 (T ), . . . , B nn 10 (T ) for the excluded volume interaction energy u nn (r ), is highly desirable.
In general, the excluded volume pressure could be obtained by means of equation (2.2), supposing that molecular excluded volume v 0 as function of both the density and temperature is known. Unfortunately, function v 0 (ρ, T ) is not available in general case. On other hand, the pressure p 0 can be obtained from the knowledge of the forces that are responsible for excluded volume. Namely, similarly to the case of the original van der Waals theory, when excluded volume pressure p 0 was identified with the pressure p hs of the fluid system with a hard-sphere repulsion u hs (r ), the excluded volume pressure within the augmented van der Waals theory can be obtained as the pressure p nn of the fluid system with interaction potential u nn (r ). These data can be extracted, for instance, from computer simulation experiment [7] or within the integral equation theories [9,10]. Such a route has been already explored [29,30], resulting in the augmented van der Waals EOS for the LJ-like HCAY fluid. The other possibility might be to utilize the excluded volume pressure p 0 in the framework of the perturbed virial EOS approach [33,34].  1] критерiй розбиття потенцiалу парної взаємодiї на двi частини, одна з яких описує як можна точнiше виключений об'єм в системi, приводить до виразiв для вiрiальних коефiцiєнтiв, якi суттєво покращують точнiсть вiрiального рiвняння стану в цiлому та для температур нижчих за критичну, зокрема. Як приклад, розглянуто застосування до Леннард-Джонсiвської моделi твердих сфер з притяганням Юкави.