Cavity correlation and bridge functions at high density and near the critical point: a test of second-order Percus–Yevick theory

ABSTRACT Cavity correlation functions, pair correlation functions, and bridge functions for the Lennard-Jones fluid are calculated from first Percus–Yevick (PY) theory and second-order Percus– Yevick (PY2) theory, molecular dynamics, and grand canonical Monte Carlo techniques. We find that the PY2 theory is significantly more accurate than the PY theory, especially at high density and near the critical point. The pair correlation function near the critical point has the expected slowly decaying long-range behaviour. However, we do not observe any long-range behaviour in the bridge function for the state points near the critical point we have simulated. However, we do note that the bridge function, which is usually negative near r = 0, becomes positive as r → 0. This behaviour is seen for the bridge functions computed from both PY2 and molecular dynamics, but not from PY.


Introduction
This paper is dedicated to Yiping Tang and his family in honour of his important contributions to the statistical mechanics of fluids and in memory of his life that was cut short. Sadly, he was unable to obtain a research position. Much of his research was accomplished in his spare time. We cannot help but wonder what might have been accomplished had he been able to obtain an academic or research position and had he lived longer.
The study of the critical point and the region near the critical point is fascinating. In this region, the pair correlation function (PCF) has an exceedingly long tail. This results in interesting singularities in some of the thermodynamic functions at the critical point. Simple theories, which assume the thermodynamic functions are analytic CONTACT J. Karl Johnson karlj@pitt.edu Supplemental data for this article can be accessed at http://dx.doi.org/./...
at the critical point, predict singularities. However, the nature of the singularities seen in experiment and in exact and near exact studies is quite different from that predicted by such simple theories.
Although the nature of the PCF, g(r), in the region of the critical point is reasonably well understood through scaling and renormalisation group studies, little is known about the direct correlation function (DCF), c(r). Even less is known about the cavity correlation function (CCF), y(r), at small r, where r is the separation between the centres of a pair of molecules. The CCF is defined as where β = 1/kT (k and T being the Boltzmann constant and absolute temperature, respectively) and u(r) is the pair interaction potential of the fluid. The bridge function, B(r), is defined to relate g(r), c(r), and y(r), The bridge function is of interest because various integral equation theories of fluids are based on approximations, called closures, to B(r). For example, the hypernetted chain (HNC) approximation is obtained from the neglect of B(r), and the Percus-Yevick (PY) approximation is obtained from Note that Equation (4) means that the PY B(r) must always be negative. The PY and HNC closures have well-known deficiencies. These can be illustrated by considering a hard sphere fluid. For hard spheres, the PY approximation gives reasonable results for the PCF and DCF. However, the PY values of y(r) are much too small in magnitude near r = 0. For hard spheres, the PY theory gives where η = πρd 3 /6 (ρ = N/V is the number of hard spheres divided by the volume and d is the hard sphere diameter). This result is in poor agreement with the result The first equality in Equation (6) is exact. In contrast, the second equality is approximate and is based on the Carnahan-Starling empirical expression. The PY values of B(r) for hard spheres are fairly accurate because g(r) and c(r) computed from PY theory are fairly accurate for hard spheres. However, the accuracy of B(r) can only be the result of a cancellation of errors between y(r) and Equation (4). Near r = 0, the HNC closure is somewhat better. At least the HNC approximation gives a result for y(0) that is exponentially large. However, for r > d the HNC result for g(r) is worse than the PY result.
To be sure, the thermodynamic functions are determined largely by the value of g(r) near the minimum in the pair interaction. One could adopt the extreme view that the nature of y(r) and B(r) near r = 0 is unimportant. However, we do not take this view. We note that for many systems, the HNC theory gives better results than the PY theory for g(r) even outside the molecular core. In the case of such systems, the error in the PY y(r) at small r, which is universal, persists to the larger r in the vicinity of the minimum in the pair interaction. It is probable that the development of improved closures requires reasonable values of B(r) at all r.
The PY and HNC theories tend to be best at low densities. Off hand, one might expect these theories to be useful for the critical region since the critical density is fairly low. This is true in part; the PY theory locates the critical point quite well. However, Henderson and Murphy [1] have shown that the PY g(r) has the form suggested by the Ornstein-Zernike (OZ) theory of the critical point, whereas the studies of Fisher [2] suggest that the power of r in the denominator should be greater than 1. As a result, the nature of the singularities in the PY results does not match what is found in experiment. Studies of the HNC in the critical region [3,4] indicate that the HNC equation ceases to yield solutions before the critical point is reached. Whether this is only a numerical problem or a real deficiency of the approximation is not known. In any case, because of this problem, the HNC theory does not yield much insight into the critical region.
Another difficulty is the fact that each of these approximations yields different results for the thermodynamic functions when differing routes relating thermodynamics and correlation functions are employed. As a result, each theory can have more than one critical point. The critical point derived from the compressibility equation is probably the most fundamental, because it results from the long-range behaviour of the correlation function.
The usual approach to the development of a theory is to state a closure and to examine its consequences. This is an ad hoc procedure. A preferable procedure would be to employ an approximation that produces a systematically improvable bridge function. The second-order Percus-Yevick (PY2) theory developed by Attard [5] and examined by Henderson and Sokołowski [6][7][8] is such a procedure. The PY2 theory has been tested, mostly at high densities, for hard spheres and Lennard-Jones (LJ) molecules with good results.
Our goal here is to study the PY2 approximation at high densities and in the critical region. We do not expect this equation to describe the critical point singularities correctly. What we hope is that this theory will generate qualitatively reasonable bridge functions. The PY2 results for the correlation functions will be compared with simulation results. Hopefully, a knowledge of the bridge function will suggest new and better closures not only near the critical point but also for other regions.
It is possible to develop a second-order hypernetted chain (HNC2) approximation that is analogous to the PY2 approximation. Perhaps this would allow one to reach the critical point and be a significant step beyond the HNC approximation. We do not do this here because it is harder to achieve solutions with the HNC2 equation than with the PY2 equation. Even so, this might be a useful future study.
All our results are obtained with the LJ pair interaction, where ϵ is the well depth of the potential and σ is the diameter of the molecule, where u(r = σ ) is 0. The LJ fluid has been studied extensively. The nature of the critical point for a fluid is universal. Thus, little would be gained by studying other potentials.

Theory
The usual theories (PY and HNC, for example) are based on the combination of the OZ relation, with a closure. The second-order theories are based on the second-order Ornstein-Zernike (OZ2) relation, (10) The OZ2 relation comes from the theory of inhomogeneous fluids and reflects the fact that in an inhomogeneous fluid the density is not constant. Attard's idea is that even in a homogeneous fluid the density, is locally inhomogeneous and Equation (10) is applicable. The source of the inhomogeneity is one of the molecules of the fluid (we have given it the subscript 4). Note that the function ρ(r 3 ) is a pair function. The functions, are not pair functions but are triplet functions. An approximation results from the combination of the OZ2 relation with a closure. Here we use the PY closure and obtain the PY2 approximation. However, this is not enough. A relation between the pair and triplet functions is required. We use the Lovett-Mou-Buff-Wertheim equation [9,10] to provide this relation, The method used to solve the PY2 equations is outlined in the papers of Attard and Henderson and Sokołowski [5][6][7]. In calculating the thermodynamic properties from the PY2 equations, we use only the compressibility route, as only for this route are the singularities in the thermodynamic functions a direct result of the long-range nature of h(r) in the critical region.

Simulation Details
As a test of the PY2 equations, we have performed molecular dynamics (MD) and grand canonical Monte Carlo (GCMC) simulations to compute the PCF and CCF at various state points. The CCF was computed from Lee and Shing's adaptation of the method developed by Henderson [11,12]. The GCMC method was only used to compute y(r) at r = 0. The MD method was used to compute y(r) for 0 < r ࣘ σ . The PCF was also computed from the MD simulations, and y(r) for r > σ was computed from the PCF through Equation (1). No significant differences between the GCMC and MD simulations were seen, so only results from MD simulations are reported here. MD simulations were carried out using an in-house code and the LAMMPS software package [13] in the isothermal-isobaric (NpT) and canonical (NVT) ensembles. The PCF was computed from the LAMMPS trajectories using the VMD [14] plug-in utilising GPU (Graphics Processing Unit) acceleration. The bridge function was calculated using Equation (2) for MD and PY2, and Equation (4) for the PY results. The DCF was computed from the convolution theorem. MD simulations were first carried out in the NpT ensemble using 4000 atoms run for 250,000 time steps with a 10σ cut-off and long-range corrections. These small simulations were carried out in order to determine the vapour and liquid densities near the near critical point. Simulations were carried out for T * = 1.30 and 1.31 (where T * = kT/ϵ) over a range of pressures to identify densities that are close to the saturation conditions at those temperatures but clearly in the onephase region. The densities and pressures were compared with saturation densities and pressures as predicted from the LJ equation of state [15] as a check. Care was taken to avoid state points in or too close to the two-phase region because of spurious behaviour we noted in the correlation functions in simulations that were in or too close to the two-phase region. These densities were then used to carry out production runs for the PCF in the NVT ensemble for a system containing 500,000 atoms, also with a 10σ cut-off and long-range corrections. We verified that the cut-off was sufficiently large by comparing densities from NpT simulations and PCF values from NVT simulations with test runs using a cut-off of 20σ , where we found no significant differences. The production runs were equilibrated for 500,000 timestep followed by data collection for 1,000,000 time steps.

Results
As a test of the accuracy of the PY2 theory, we have computed the correlation functions at T * = 1.0 and ρ * = 0.845 (ρ * = ρσ 3 ), which is a high density liquid state point. Our PY2 code did not yield convergent results for lower values of T * or higher values of ρ * at T * = 1.0. The PCF and CCF computed from MD, PY2, and PY are plotted in Figure 1. We see that the PY2 theory is in good agreement with the MD simulation results, slightly underestimating the height of the first peak, but capturing the subsequent peaks with quantitative accuracy. In contrast, the PY theory overestimates the height of the first peak and qualitatively misses the shape of the second peak, having a first minimum past the first peak at a value of r that is too small, the second peak being both too large and at a value of r that is too large. In this work we have plotted the inverse hyperbolic sine function of the CCF, sinh −1 [y(r)], in order to show both the large and small values of the CCF with good resolution. We see from Figure 1(b) that PY dramatically underestimates the value of the CCF as r approaches 0. The PY2 theory also underestimates the CCF compared with the simulations, but is qualitatively correct. Figure 1 demonstrates that PY2 gives a major improvement in accuracy compared with PY for dense liquids.
We estimate that the critical point from the PY2 theory is T * c = 1.31, ρ * c = 0.27. This value was obtained by computing the bulk modulus or inverse compressibility from the compressibility equation in the critical point region and seeking a minimum in the resulting paraboliclike curves that occurred when the minimum of the bulk modulus that occurred for vanishing bulk modulus. Because a solution was not always found, analytic continuation was needed in one or two instances. The value of the critical point from the PY theory (obtained from the compressibility route) is T * c = 1.31 and ρ * c = 0.278, as reported by Henderson and Murphy [1]. The critical point of the LJ fluid has been estimated from simulations using mixed-field finite-size scaling calculations giving values of about T * c = 1.313, ρ * c = 0.317 [16,17]. The critical temperature estimated from PY2 is very close to the simulations, but the PY2 critical density is about 15% lower than that estimated from simulations.
MD simulations carried out at the PY2 critical point are in the two-phase region and therefore cannot provide accurate values of the correlation functions. Likewise, simulations at the critical point estimated from simulations will be subject to large fluctuations. We therefore computed the correlation functions at two temperatures close to criticality, T * = 1.30 and 1.31, but at densities that are in the stable vapour and liquid regions, as estimated from NpT simulations. The densities chosen were ρ * = 0.146 and 0.462 for T * = 1.30 and ρ * = 0.177 and 0.460 for T * = 1.31.
The vapour branches show the most interesting behaviour of these four state points. The PCFs, CCFs, and bridge functions are plotted in Figure 2 for the T * = 1.30, ρ * = 0.146 state point. Results from simulations, PY2, and PY are plotted. The PCFs computed from each of the three methods are in good agreement, as seen in Figure 2(a). The CCFs computed from MD and PY2 are in excellent agreement, while PY significantly underestimates the CCF for r * < 0.5, as seen from the inset of Figure 2(b). The behaviour of the bridge functions, shown in Figure 2(c), is the most interesting. We see that both MD and PY2 show that B(r) > 0 for r → 0, while the bridge function from PY is always negative, as required by Equation (4). Similar results can be seen for the correlation functions at the T * = 1.31, ρ * = 0.177 state point shown in Figure 3. The major difference between these two state points is that B(r) computed from MD is slightly negative rather than positive. We see from the insets of  Figure 4 and are given in Figure S1 of the supplemental information. The PY theory gives errors in the locations of the minima and maxima of the PCF, compared with both MD and PY2 data. Specifically, the PY theory gives locations of the first minimum after the first peak at values of r that are slightly too small and second peak positions that are at values of r that are slightly too large. We note that the PCFs for the low-density (vapour phase) state points at both T * = 1.30 and 1.31 (see Figures 2(a) and 3(a)) show the expected gradually decaying behaviour at large r. In contrast, g(r) for T * = 1.30, ρ * = 0.462 decays to zero at relatively small values of r (see Figure 4(a)).
The CCFs for liquid state points at T * = 1.30 and 1.31 are shown in Figures 4(b) and S1(b) in the supplemental information, respectively. As with the vapour density state points, the PY theory significantly underestimates y(r) as r → 0. At larger values of r, the MD, PY2, and PY results are virtually indistinguishable on the scale of the plots.
The bridge functions at T * = 1.30, ρ * = 0.462 and T * = 1.31, ρ * = 0.460 are shown in Figures 4(c) and S1(c) of the supplemental information, respectively. In contrast to the vapour densities at these temperatures (see Figures 2(c) and 3(c)), B(r) computed from PY2 theory is always less than 0 as r → 0. As with other state points, the PY theory predicts much more negative values of B(r) as r → 0.   We note that, in contrast to the PCF, which becomes long ranged near the critical point, we found no evidence of long-range behaviour in B(r), at least for the state points we were able to simulate. In order to test whether the PY2 theory predicts long-range behaviour in B(r) near the PY2 critical point, we computed the correlation functions for a state point of T * = 1.31 and ρ * = 0.25, which is the closest we could get our code to converge to the critical point. The correlation functions are shown in Figure S3 of the supplemental information. We found that while g(r) decays slowly, none of the other correlation functions, including B(r), are long ranged. This indicates that simple approximations for B(r) should be useful even near the critical point.
We have also computed the correlation functions for a supercritical state point, T * = 1.4 and ρ * = 0.28, in order to assess the influence of near criticality on the correlation functions. The PCF, CCF, and bridge functions for this state point are plotted in Figure S2 of the supplemental information. The major difference between the near and supercritical state points is that B(r) is clearly <0 for all r at the supercritical condition. The PCF and CCF T * = 1.4 are similar to the vapour state points at T * = 1.30 and 1.31. As with the other state points, the PY theory underestimates y(r) as r → 0 and predicts B(r) values that are significantly too negative in this same range.
There is some experimental evidence that y(0) might diverge near the critical point [18]. This conjecture is based upon the observation that the infinite dilution partial molar volume of naphthalene becomes large and negative as the critical point of carbon dioxide is approached [18]. This is somewhat analogous to the situation for computing y(0) in our case. For a soft sphere fluid βμ ∞ t = ln y(0), where μ ∞ t is the infinite dilution chemical potential of a double strength test particle, i.e. a particle having the same size but twice the interaction strength (well depth) as the solvent molecules. The infinite partial molar volume for this case is given bȳ where i is the test particle and j denotes the solvent. This would indicate that y(0) changes very rapidly with p near the critical point and that as p is increased beyond criticality (density increased) that y(0) decreases. While this conjecture may be valid, we note that we did not see any evidence for divergence of y(0) in our calculations.

Conclusions
We have presented a test of the second-order Percus-Yevick theory by comparing the CCFs and bridge functions computed from MD simulations with those computed from PY2 and PY theories at a high density state point (T * = 1.0, ρ * = 0.845) and at conditions near the critical point. The most striking feature is that B(r) > 0 when r is close to 0 for points near criticality, as computed from both PY2 theory and MD simulations. We have not found any evidence that the bridge function is long ranged near the critical point. This may be due to our limitations of carrying out simulations that are not too close to the critical point. We have found that the PY2 theory is significantly more accurate than the PY theory for dense liquids. Indeed, the PY theory fails to give qualitatively correct pair and CCFs compared with MD simulations for a state point corresponding to T/T c ∼ 0.76 and p/p c ∼ 3.