The issue of Frontiers in Physics this paper belongs to is devoted to commemorating the celebration of fifty consecutive annual Winter Meetings in Statistical Physics in Mexico. Therefore, we have chosen to write on a subject that has been present in these meetings from the beginning; namely, the thermodynamic properties of fluids that we are persuaded can still offer some interesting results.

We begin by recalling that, in an attempt to prove the validity of the thermodynamic limit of classical statistical mechanics, van Hove [1] introduced in 1949 a potential ϕ ( r ) consisting of a hard core of radius r 0 and a finite-range attractive tail. The actual form of this so-called van Hove potential is ϕ ( r ) = { ∞ , 0 ≤ r ≤ r 0 , < 0 , r 0 < r ≤ b > − ε 0 , r 0 < r ≤ b , 0 , r > b , (1) where r is the distance, b corresponds to the range, and − ε 0 corresponds to the lower bound of the attractive tail, whose form is rather arbitrary. It should be pointed out that two popular models of intermolecular potentials used in liquid state physics, namely, the triangle-well potential and the square-well potential, fulfill the condition of being van Hove potentials and their thermodynamic properties have been thoroughly studied (see, for instance, Refs. [2–11] for the former model and Refs. [12–24] for the latter and references therein). Surprisingly, as far as we know, the parabolic-well potential, which is also a van Hove potential, has not been used for that purpose. The main aim of this paper is to contribute to partly remedying this situation.

We consider a parabolic-well fluid whose molecules interact with a potential of the form u ( x ) = { ∞ , 0 ≤ x ≤ 1 , ε [ ( x − 1 λ − 1 ) 2 − 1 ] , 1 < x ≤ λ , 0 , x > λ , (2) where x = r / σ is the reduced distance (r being the distance), σ is the diameter of the hard core, ε > 0 is the well depth, and λ > 1 is the potential range. As it occurs with other relatively simple models, the main asset of this model potential is probably that, despite being an idealized representation, it nevertheless contains the main features of true molecular interactions in fluids, namely, a repulsive hard-core and an attractive interaction that continuously goes to zero as the intermolecular distance increases. In this regard, it is interesting to recall what Widom [25] pointed out in the case of the square-well fluid: “Where I speak of the necessity to treat accurately the effects of the attractive or repulsive forces, I do not mean that it is important to know the corresponding part of ϕ ( r ) with quantitative accuracy. Indeed, even if the ϕ ( r ) of Figure 1 were idealized as a square-well potential, as in Figure 3, but the statistical mechanical consequences of such a potential were then determined without further approximation, there would undoubtedly result in an essentially correct description of all the macroscopic properties of matter throughout a vast region of the p and T plane, including the neighborhoods of the triple and critical points. Thus, what matters is not the quantitative accuracy of the assumed ϕ ( r ) , but rather the qualitative accuracy of the resulting spatial correlations of molecular positions; the triple and critical points are distinguished by having the relevant qualitative features of this correlation, and the nature of its propagation through the fluid, determined primarily by the short-range repulsive forces between molecules, or by the longer ranged attractive forces, respectively.” Something similar may be said about the parabolic-well potential. In fact, an interesting asset of this model is that its thermodynamic properties are readily amenable for treatment within the second-order thermodynamic perturbation theory of Barker and Henderson [26]. Within this approach, in order to derive the Helmholtz free energy of the parabolic-well fluid, two ingredients are required: on the one hand, one needs the Helmholtz free energy of the reference hard-sphere fluid of diameter σ. On the other hand, one also requires an expression for the radial distribution function g H S ( x ) of the hard-sphere fluid. In this work, we will profit from the availability of a method [27], the so-called rational function approximation (RFA) method, to (analytically) obtain an approximate g H S ( x ) which is thermodynamically consistent with the equation of state of the hard-sphere fluid and make use of this fact to derive the equation of state of the parabolic-well fluid.

The paper is organized as follows. In the next section, we recall the main aspects of the RFA method for the computation of g H S ( x ) and provide the explicit expression for this quantity in the first coordination shell. This is followed in Section 3 with the completely analytic derivation of the equation of state of the parabolic-well fluid within the second-order Barker-Henderson thermodynamic perturbation theory in the macroscopic compressibility approximation and taking the hard-sphere fluid as the reference system. Section 4 contains some illustrative results for the reduced pressure of the parabolic-well fluid and a comparison with our own Monte Carlo and Event-driven Molecular Dynamics simulation data. We close the paper in the final section with further discussion and some concluding remarks.

In this section, we provide the analytic result for the radial distribution function (rdf) of the hard-sphere fluid g H S ( x ) , as derived with the RFA method [27], and its explicit expression in the range 1 < x ≤ 2 . We begin by recalling two important relationships between the thermodynamic and structural properties of the hard-sphere fluid derived from statistical mechanics. On the one hand, the compressibility factor Z H S = p ρ k B T (where p is the pressure, ρ the number density, k B the Boltzmann constant, and T the absolute temperature) of the hard-sphere fluid is related to the contact value of the rdf g H S ( 1 + ) through Z H S = 1 + 4 η g H S ( 1 + ) , (3) where η = π ρ σ 3 / 6 is the packing fraction. On the other hand, the hard-sphere isothermal susceptibility χ H S ≡ [ d ( η   Z H S ( η ) ) d η ] − 1 is related to the rdf through χ H S = 1 + 24 η ∫ 0 ∞ d x   x 2   [ g H S ( x ) − 1 ] . (4)

In the RFA method [27], the Laplace transform of x   g H S ( x ) is taken to be given by G ( t ) = ℒ { x   g H S ( x ) } = t 12 η 1 1 − e t Φ ( t ) , (5) where Φ ( t ) = ( 1 + S 1 t + S 2 t 2 + S 3 t 3 + S 4 t 4 ) / ( 1 + L 1 t + L 2 t 2 ) , and the six coefficients S 1 , S 2 , S 3 , S 4 , L 1 , and L 2 (which depend on the packing fraction) may be evaluated in an algebraic form by imposing the following requirements: (i) χ H S must be finite and hence the first two integral moments of the total correlation function h ( x ) ≡ g ( x ) − 1 , i.e., ∫ 0 ∞ d x   x n h ( r ) with n = 1,2 , must be well defined; (ii) the approximation must be thermodynamically consistent with a prescribed equation of state; i.e., the thermodynamic relationship χ H S = [ d η   Z H S ( η ) d η ] − 1 must be satisfied. Using the first requirement, one finds that L 1 , S 1 , S 2 , and S 3 are linear functions of L 2 and S 4 . Imposing the requirement (ii) leads to explicit expressions for L 2 and S 4 in terms of χ H S and Z H S [27]. Finally, the expressions for all the coefficients are as follows: L 1 = 1 2 η + 12 η L 2 + 2 − 24 η S 4 2 η + 1 , (6) S 1 = 3 2 η − 1 + 4 L 2 − 8 S 4 2 η + 1 , (7) S 2 = − 1 2 − η + 8 η L 2 + 1 − 2 L 2 − 24 η S 4 2 η + 1 , (8) S 3 = 1 12 2 η − η 2 + 12 η 2 L 2 − 12 η L 2 − 1 − 72 η 2 S 4 ( 2 η + 1 ) η , (9) L 2 = − 3 ( Z H S − 1 ) S 4 , (10) S 4 = 1 − η 36 η ( Z H S − 1 / 3 ) [ 1 − [ 1 + Z H S − 1 / 3 Z H S − Z P Y ( χ H S χ P Y − 1 ) ] ] 1 / 2 . (11)

Here, Z P Y = 1 + 2 η + 3 η 2 ( 1 − η ) 2 and χ P Y = ( 1 − η ) 4 ( 1 + 2 η ) 2 are the compressibility factor and isothermal susceptibility arising in the Percus-Yevick theory. To close the problem, one has to give an expression for Z H S , so all the procedure is a function of this choice. For a given Z H S , the radial distribution function is given by g H S ( x ) = { 0 , 0 ≤ x < 1 , 1 12 η x ∑ n = 1 ∞ φ n ( x − n ) θ ( x − n ) , x ≥ 1 (12) with θ ( x − n ) being the Heaviside step function and φ n ( x ) = ℒ − 1 { − t [ Φ ( t ) ] − n } . (13)

Explicitly, using the residues theorem, φ n ( x ) = − ∑ n = 1 4 e t i x ∑ m = 1 n A m n ( t i ) ( n − m ) ! x n − m , (14) where A m n ( t i ) = lim t → t i 1 ( m − 1 ) ! ( d d t ) m − 1 ( t − t i ) t [ Φ ( t ) ] − n , (15) with t i being the four roots of 1 + S 1 t + S 2 t 2 + S 3 t 3 + S 4 t 4 = 0 ; namely, t 1 = − S 3 4 S 4 + y p − y n , (16) t 2 = − S 3 4 S 4 + y p + y n , (17) t 3 = − S 3 4 S 4 − y p − y m , (18) t 4 = − S 3 4 S 4 − y p + y m , (19) where y p = − 1 2 S 3 2 4 S 4 2 − 2 S 2 3 S 4 + y r + y s 3 S 4 , (21) y r = − S 2 2 − 3 S 1 S 3 + 12 S 4 3 S 4 y s , (22) y n = 1 2 S 3 2 4 S 4 2 − 4 S 2 3 S 4 − y r − y s 3 S 4 − S 3 3 8 y p S 4 3 + S 2 S 3 2 y p S 4 2 − S 1 y p S 4 , (23) y m = 1 2 S 3 2 4 S 4 2 − 4 S 2 3 S 4 − y r − y s 3 S 4 + S 3 3 8 y p S 4 3 − S 2 S 3 2 y p S 4 2 + S 1 y p S 4 , (24) y s = − 1 2 1 / 3 [ y t + − 4 ( S 2 2 − 3 S 1 S 3 + 12 S 4 ) 3 + y t 2 ] 1 / 3 , (25) y t = 2 S 2 3 + 9 S 1 S 2 S 3 + 27 S 3 2 + 27 S 1 2 S 4 − 72 S 2 S 4 . (26)

As we will indicate below, once Z H S ( η ) has been chosen, Eqs. 5–26 are all that is needed to evaluate the first- and second-order perturbation terms for the free energy of the parabolic-well fluid within the Barker-Henderson thermodynamic perturbation theory taking the hard-sphere fluid as the reference system. To close this section and for later use, we now write the explicit expression for the radial distribution function up to the first coordination shell which reads g H S ( x ) = { 0 , 0 ≤ x < 1 , ∑ i = 1 4 a i x e t i ( x − 1 ) , 1 ≤ x ≤ 2 (27) where a i = − t i ( 1 + L 1 t i + L 2 t i 2 ) 12 η ( S 1 + 2 S 2 t i + 3 S 3 t i 2 + 4 S 4 t i 3 )   ( i = 1,2,3,4 ) . (28)

Perturbation approaches for the computation of thermodynamic properties of fluids are well established theoretical tools [28, 29]. In the Barker-Henderson perturbation theory [26], one splits the potential into a hard-sphere part and a perturbation part; namely, u ( r ) = u H S ( r ) + u 1 ( r ) , where u H S ( x ) = { ∞ , 0 ≤ x ≤ 1 , 0 , x > 1 (29) and u 1 ( x ) = { 0 , 0 ≤ x ≤ 1 , ε [ ( x − 1 λ − 1 ) 2 − 1 ] , 1 < x ≤ λ , 0 , x > λ (30)

Once this separation has been made, the Helmholtz free energy per particle of the parabolic-well fluid is expressed as a power series in the inverse of the reduced temperature T ∗ = k B T / ε , which up to second order reads f N k B T = f H S N k B T + 1 T ∗ f 1 N k B T + 1 T ∗ 2 f 2 N k B T (31)

Here, N is the number of particles and f H S stands for the Helmholtz free energy of the reference hard-sphere fluid while f 1 and f 2 (this latter in the so-called macroscopic compressibility approximation) are given, respectively, by f 1 N k B T = 12 η ε ∫ 1 λ g H S ( x ) u 1 ( x ) x 2 d x (32) and f 2 N k B T = 6 η χ H S ε 2 ∫ 1 λ g H S ( x ) u 1 2 ( x ) x 2 d x . (33)

Note that we have made use of the fact that g H S ( x ) vanishes for 0 ≤ x < 1 and of the expression for u 1 ( x ) given in Eq. 30, respectively, to set the lower and upper limits of the integrals in Eqs. 32, 33. In turn, the equation of state of the parabolic-well fluid in this approximation is given by Z ≡ p ρ k B T = Z H S + 1 T ∗ η ∂ ∂ η ( f 1 N k B T ) + 1 T ∗ 2 η ∂ ∂ η ( f 2 N k B T ) . (34) And, the chemical potential may be readily obtained as μ k B T = f N k B T + Z , (35) where f N k B T is given in Eq. 31, together with Eqs. 32, 33, and Z is given in Eq. 34. So, provided we choose Z H S , which of course also determines χ H S and f H S , and take g H S to be the one computed with the RFA approach and such compressibility factor, the completely analytic formulation of the second-order Barker-Henderson thermodynamic perturbation theory in the macroscopic compressibility approximation for the parabolic-well fluid taking the hard-sphere fluid as the reference system has been derived. In our subsequent calculations, we will be restricted to relatively narrow wells ( 1 < λ ≤ 2 ) so that Eq. 27 for g H S ( r ) will be used. Furthermore, Z H S and χ H S will be chosen to be those corresponding to the Carnahan-Starling (CS) equation of state [30]; namely, Z H S ( η ) ≡ Z C S ( η ) = 1 + η + η 2 − η 3 ( 1 − η ) 3 (36) and χ H S ( η ) ≡ χ C S ( η ) = ( 1 − η ) 4 1 + 4 η + 4 η 2 − 4 η 3 + η 4 . (37)

Further, from the CS equation of state, it also follows that f H S ( η ) N k B T = f C S ( η ) N k B T = − 1 + ln 6 η π + ( 4 − 3 η ) η ( 1 − η ) 2 . (38)

The availability of the completely analytic (albeit approximate) forms of the Helmholtz free energy and the equation of state of the system (which are themselves not very illuminating and therefore will not be explicitly written down [31]) allows us in principle to compute, for a given value of λ, the compressibility factor using Eq. 34, the vapor-liquid coexistence curve from the equality of pressures, and chemical potentials of the two phases and also to obtain the critical point in the usual way.

Preliminary results for the isotherms will be presented in the following section, together with a comparison with our simulation data. But before presenting such results, we will take advantage of the simple form of the intermolecular potential of this fluid to compute its second virial coefficient. This is given by B 2 ( T ) ≡ − 2 π σ 3 ∫ 0 ∞ x 2 ( e u ( x ) k B T − 1 )   d x (39) = − π 6 σ 3 [ 6 k B T ε ( 1 − 2 e ε k B T + λ ) ( λ − 1 ) 2 − 4 λ 3 + 3 π k B T e ε k B T ε 3 / 2 ( 2 ε + k B T ( λ − 1 ) 2 ) ( λ − 1 ) E r f ( ε k B T ) ] .

Equation 39, which to the best of our knowledge has not been reported before, allows us to obtain the Boyle temperature T B of the parabolic-well fluid as a function of λ by equating B 2 ( T B ) to zero and solving numerically for T B . In Table 1, we show some particular values and, for comparison, we also include the values corresponding to triangle-well and square-well fluids with the same range λ.

To close this section, we will also take advantage of the knowledge of the second virial coefficient, to obtain estimates of the critical temperature according to the Vliegenthart and Lekkerkerker criterion [32], namely, from equating this coefficient with − 6 v m , where v m = π 6 σ 3 is the volume of the spherical core. The results for given values of the range are given in Table 2, where we have also included such estimates for the cases of the triangle-well and square-well fluids with the same range λ.

Note that, for all three model fluids, the values of both the reduced Boyle temperatures and the estimates of the reduced critical temperatures increase as the range λ is increased. Also note that the geometrical form of the well influences such values as reflected in the fact that, for the same value of the range, the ones corresponding to the triangle-well fluid are smaller than those of the parabolic-well fluid which, in turn, are smaller than those of the square-well fluid.

Now we return to our main aim. In order to assess the value of the thermodynamic perturbation theory approach presented in the previous section, we have carried out NVT Monte Carlo (MC) simulations to compute the pressure of parabolic-well fluids for various values of the range λ ≤ 2 and supercritical temperatures for later comparison with our theoretical results. The details of such simulations are as follows. The number of particles in our simulations is N = 1372 and we have considered a cubic box, of length L and with periodic boundary conditions. Reduced units are used, so that lengths are expressed in units of σ ( L = l σ , with l a pure number), the reduced temperature is T ∗ , the packing fraction is η = π 6 N V σ 3 (where V = L 3 is the volume), and the reduced pressure is p ∗ = p σ 3 ε − 1 .

For the sake of illustration, we report here the results of the simulations for λ = 1.25 and 1.75, along various isotherms. For each isotherm, eight different packing fraction values η = 0.05 to 0.5 with Δ η = 0.05 were simulated in order to compute the reduced pressure p ∗ . Each run was carried out using 1.5 ⋅ 10 6 Monte Carlo steps (MCS) discarding the first 10 6 MCS for equilibration, and the properties were measured every 20 MCS and averaged every 1000 MCS; furthermore, for each packing fraction, the values of p ∗ were averaged over 20 parallel simulations to obtain better statistics.

Finally, the pressure was calculated using the expression p ∗ = 6 π η T ∗ [ 1 + 4 η g P W ( 1 + ) − 4 η T ∗ ∫ 1 λ g P W ( x ) d u 1 ∗ ( x ) d x x 3 d x ] . (40) Here, u 1 ∗ ( x ) = u 1 ( x ) / ε , g P W ( x ) is the radial distribution function of the parabolic-well fluid (computed in the usual way [33] with the subscript PW standing for parabolic well), and the second term on the right-hand side of Eq. 40, obtained following a similar procedure to the one used by Rotenberg [12] in the case of the square-well fluid, accounts for the hard-core contribution to the parabolic-well potential.

In Figures 1‐3, we show the comparison between the results of the isotherms obtained with the thermodynamic perturbation theory and from simulation. Note the good agreement between theoretical and simulation results for all the values of T ∗ above the critical temperature that we considered.

On the other hand, the subcritical isotherms were obtained by means of Molecular Dynamics (MD) Event-driven simulations. We have performed event-driven simulations of N = 108000 elastic smooth spheres carried out with the DynamO software package [34]. The spheres interact by a stepped parabolic-well type potential [35–37], a discretized version consisting of a sequence of 15 steps of widths 0.05 σ , steps more than reasonable in most instances [37]. We have used σ, τ = m σ 2 / ε , m, and T ∗ as units of length, time, mass, and temperature, respectively. The MD event-driven simulations were performed for one value of the range of the potential, namely, λ = 1.75 , along the isotherms T ∗ = 1.1 ,   1.0 and 0.7. In the first stage, we have performed NVT simulations with an Andersen thermostat during 2 × 10 8 collisions, and after the equilibration, the second stage of NVE simulations with a duration of 8 × 10 8 collisions was performed [38]. The full pressure tensor for the system was determined by means of the expression P = P k i n e t i c + P int e r a c t i o n (41) where the kinetic pressure is given by P k i n e t i c = 1 2 V ∑ i N v i v i , (42) with v i v i being a dyadic product which yields a matrix result and the masses of the particles are set as m i = 1 . The contribution to the pressure due to interactions is given by P interaction = 1 V t s i m ∑ i j e v e n t Δ p i r i j , (43) where the summation is over each two-particle i , j event interaction, t s i m is the total simulation time, Δ p i is the momentum impulse on particle i, and r i j = r i − r j is the separation vector between the interacting particles. Finally, the hydrostatic pressure, which in this instance coincides with the reduced pressure, was computed from the trace of the tensor p ∗ = t r ( P ) / 3 = ( P x x + P y y + P z z ) / 3 . (44)

In Figure 2, we show the comparison of the subcritical isotherms obtained with the thermodynamic perturbation theory and from simulation for a range λ = 1.75 . The typical van der Waals loop is clearly seen for the theoretical isotherm with T ∗ = 0.7 (which grossly underestimates the simulation data) and is still present in the isotherms with T ∗ = 1 and T ∗ = 1.1 , respectively. On the other hand, we have checked both through Monte Carlo and Event-driven MD simulations and also through the outcome of the thermodynamic perturbation theory that the isotherm with T ∗ = 1.35 (not shown), which according to the Vliegenthart and Lekkerkerker criterion should be the critical one, is a supercritical isotherm. In fact, the simulation data indicate that the real critical isotherm for this value of the range lies above but close to the one corresponding to the theoretical curve for T ∗ = 1.1 .

While it is clear from Figures 1–3 that the qualitative trends observed in all the simulation results are correctly accounted for by the theory, a better perspective of its performance may be gained by looking at the quantitative differences. Therefore, in Table 3, we display the actual numerical values for a couple of isotherms. In both cases, it is clear that the good qualitative agreement seen in Figures 1, 2, respectively, is not accompanied by quantitative agreement. In fact, the first theoretical isotherm ( λ = 1.25 and T ∗ = 1.5 ), which is a supercritical isotherm, yields an underestimation of the reduced pressure when compared to the simulation values. On the other hand, for the second isotherm ( λ = 1.75 and T ∗ = 1.1 ), which is subcritical, the general overall trend is that the theoretical curve overestimates the value of the reduced pressure. As one would expect, in the case of the supercritical isotherms, the quantitative agreement is improved as the reduced temperature is increased.

In this paper, we have addressed the study of the thermodynamic properties of a fluid whose molecules interact through a parabolic-well potential. For this model, we obtained the exact second virial coefficient which in turn allowed us to compute the Boyle temperature and to estimate the critical temperature for arbitrary values of the potential range λ. The parabolic-well potential is in the same family as the triangle-well potential and the square-well potential, being in some sense intermediate between the other two. A reflection of this is the behavior of both the Boyle temperatures and the estimates of the critical temperatures in which, for a fixed range, the values corresponding to the parabolic-well potential lie between the ones corresponding to the other two. Whether this points out to a deeper relationship between the geometrical shape of the well and the location of the critical point in van Hove fluids is not clear to us at this stage but might be worth considering in the future.

In order to obtain further analytic results, we considered a thermodynamic perturbation theory approach for this fluid within the Barker-Henderson second-order macroscopic compressibility approximation and taking the hard-sphere fluid as the reference fluid. Restricting ourselves to values of the range in the interval 1 < λ ≤ 2 and evaluating the radial distribution function of the hard-sphere fluid according to the RFA method with the CS equation of state, we were able to derive (albeit approximate) fully analytic expressions for the Helmholtz free energy, the equation of state, and the chemical potential of the parabolic-well fluid. With such expressions, we were able to compute theoretically various isotherms for a given potential range. These were subsequently compared to our own Monte Carlo NVT and Event-driven MD simulation results. It must be emphasized that these simulation data are to our knowledge the only ones available in the literature for this system.

It should be clear that the calculations that we have presented in the previous section are still preliminary but we want to stress that further work on this subject is currently being carried out. Nevertheless, at this stage, a few additional comments are in order. We begin by pointing out that the qualitative agreement between the results for the isotherms above the critical one obtained from thermodynamic perturbation theory and those stemming out of NVT Monte Carlo simulations, as well as the improvement of the quantitative agreement as the reduced temperature is increased, although clearly rewarding, are not very surprising in view of the fact that our theoretical approximation relies on the convergence of the perturbation expansion for high temperatures. Also rewarding is the fact that the results of the Event-driven MD simulation for the isotherm T ∗ = 1.0 in the case in which the range is λ = 1.75 , which is a subcritical isotherm, are also well accounted for by the curve obtained using thermodynamic perturbation theory. The same happens with the isotherm with T ∗ = 1.1 . On the other hand, the gross underestimation of the theoretical curve for the subcritical isotherm with T ∗ = 0.7 and the same value of the range indicates that the convergence of the perturbation series is very poor for this reduced temperature. In any case, it is fair to say that the present theoretical approach provides a good starting point for the study of the thermodynamic properties of parabolic-well fluids. Future work with the same approach contemplates the computation of the critical point and the liquid-vapor coexistence curve of such fluids. Finally, since we are persuaded that the parabolic-well fluid may still offer some other insights on the thermodynamic behavior of fluids, we hope that the results of the present paper may also motivate others to conduct more studies using this model.