Table Of Content

A liquid-state theory that remains successful 8 9 in the critical region 9 1 n D. Pini, G. Stell a J Department of Chemistry, State University of New York 7 at Stony Brook, Stony Brook, New York 11794–3400, U.S.A. 2 ] N. B. Wilding h c Department of Physics and Astronomy, University of Edinburgh e m Edinburgh EH9 3JZ, U.K. - t a t s . Abstract t a m A thermodynamically self-consistent Ornstein-Zernike approximation (SCOZA) - is applied to a fluid of spherical particles with a pair potential given by a hard-core d n repulsion and a Yukawa attractive tail w(r) = exp[ z(r 1)]/r. This potential − − − o allows one to take advantage of the known analytical properties of the solution to c [ the Ornstein-Zernike equation for the case in which the direct correlation function outside the repulsive core is given by a linear combination of two Yukawa tails 1 v and the radial distribution function g(r) satisfies the exact core condition g(r) = 0 7 for r < 1. The predictions for the thermodynamics, the critical point, and the 7 2 coexistence curve are compared here to other theories and to simulation results. 1 In order to unambiguously assess the ability of the SCOZA to locate the critical 0 point and the phase boundary of the system, a new set of simulations has also 8 9 been performed. The method adopted combines Monte Carlo and finite-size scaling / t techniques and is especially adapted to deal with critical fluctuations and phase a m separation. It is found that the version of the SCOZA considered here provides very good overall thermodynamics and a remarkably accurate critical point and - d coexistence curve. For the interaction range considered here, given by z = 1.8, the n o critical density and temperature predicted by the theory agree with the simulation c results to about 0.6%. : v i X r a Prepared for the John Barker festschrift issue of Molecular Physics. State University at Stony Brook College of Eng. and Appl. Sci. Report No. 754, Jan. 1998. 1 1 Introduction After applying their version of thermodynamic perturbation theory to square-well and Lennard-Jones fluids, John Barker and Doug Henderson characterized it as a “successful theory of liquids” [1]. And so it was. When tested against simulation results it proved to be impressively accurate at liquid-state densities and temperatures, unlike some versions of thermodynamic perturbation theory that had preceded it. And it bypassed the trou- bling lack of thermodynamic self-consistency associated with the direct use of the radial distribution functions obtained from the integral-equation theories then available, as well as yielding thermodynamic results as good or better than the best results obtainable from such integral equations. These positive features became hallmarks of successful thermodynamic perturbation theories for simple fluids and were shared by the versions [2] that followed the Barker and Henderson work as well as an alternative perturbative approach set forth somewhat earlier by Hauge and Hemmer [3] that was based on using the inverse range of the attractive interaction rather than its strength as a perturbation parameter. Integral-equation approaches with improved self-consistency were also developed subsequently to yield accurate liquid-state thermodynamics [4]. Unfortunately, the accuracy of all these approaches begins to decrease substantially as one leaves the liquid-state region located slightly above the triple point in temperature and follows the liquid-gas coexistence curve in the density–temperature plane up to the criticalregion. Inparticular, theshape ofthe coexistence curve andlocationofthe critical point are not accurately reproduced, nor are related critical parameters. In the case of the perturbation theories, it is not hard to understand why this is so. All of them are mean-field-like in nature, associated with coexistence curves that are quadratic close to the critical point, whereas the true coexistence curve is very nearly cubic. That is, in these theories one finds near the critical point a coexistence curve of the form T T A ρ ρ x, x = 2, (1) c c − ≈ | − | where ρ and T are the critical values of number density ρ and absolute temperature c c T, and A is a constant. In contrast, in an exact treatment, one would expect to find x close to 3. In these theories the resulting T is usually more than 5% too high and the c critical compressibility factor (P/ρk T) is usually more than 10% too high. Here P is B c the pressure, and k is the Boltzmann constant. B The thermodynamics associated with the radial distribution function g(r) obtained form various integral-equation approaches cannot be so neatly categorized. However, in the cases in which there are substantial discrepancies between the several paths available for obtaining thermodynamics from g(r), the most reliable and accurate coexistence behavior is often obtained from evaluating the thermodynamics through the excess internal energy expressed in terms of an integral over the pair potential w(r)weighted by g(r). For continuum-fluid models the resulting critical behavior is typically mean-field like in the cases that we have studied, and thus subject to the same deficiencies as one approaches the critical region. In some integral-equation approaches that have been developed in order to insure a certain degree of thermodynamic consistency, the description of the critical 2 region and of the phase diagram appears to be more problematic: for instance, the mod- ified hypernetted chain (MHNC) theory [5] is indeed able to predict quite satisfactorily the liquid and the vapor branches of the coexistence curve of a simple fluid at low enough temperature, but it fails to converge close to the critical point, so that the two branches remain unconnected, and the position of the critical point is not given directly by the theory, but must be determined by extrapolation [6, 7]. The same kind of behavior [6, 8] is found also for the HMSA integral equation [the acronym coming from the fact that the theory [9] interpolates between the hypernetted chain (HNC) and the soft mean spherical approximation (SMSA)]. The self-consistent Ornstein-Zernike approximation (SCOZA) we consider here is not mean-field-like, and it remains highly accurate as one goes from liquid-state conditions to critical-point conditions. In particular the power x in Eq. (1) was recently shown analytically to be given in the SCOZA by exactly 20/7 [10]. And as we discuss in this paper, in the hard-core Yukawa fluid (HCYF) T appears to be within 0.6% of its value c as estimated by our simulation results. (Similarly, in recent three-dimensional lattice-gas studies [11, 12] the SCOZA T was found to be within 0.2% of its estimated exact value). c As described elsewhere [10, 12] the scaling behavior of the SCOZA thermodynamics is somewhat different from the simple scaling one expects to see in the exact thermodynamics, although those differences only begin to appear clearly when ρ and T are within less than1%oftheircriticalvalues. Closer tothecriticalpoint, theeffective exponents defined above T approach spherical-model values as the critical point is approached, whereas the c exponents defined below T do not. The exponents are discussed in Sec. 3. c The SCOZA was proposed some time ago by Høye and Stell [13, 14] but fast and accurate algorithms for evaluating its thermodynamic predictions were developed only recently [11, 12, 15]. A sharp assessment of its accuracy for the HCYF could not be made on the basis of existing simulations, and for that reason our study here includes new Monte Carlo (MC) results exploiting finite-size scaling (FSS) techniques [16]. We have chosen the HCYF pair potential as the first of the continuum-fluid potentials to be considered in our studies of the SCOZA for several reasons. First, it embodies the two key features one requires in an off-lattice pair potential in order to consider both the liquid state and liquid-gas criticality—a highly repulsive core and an attractive well. Second, the HCYF proves to be particularly convenient to analyze using the SCOZA (the square-well fluid is far less convenient in this regard). Third, the functional form of the hard-core Yukawa potential makes it appropriate as a generic solvent-averaged interaction potential between polyelectrolytes and colloids as well as a generic simple- fluid pair potential. For this reason it seems particularly useful to have an accurate theory for both the structure and thermodynamics of the HCYF, which has already been the subject of a number of previous studies. We shall make contact with several of those here. The paper is organized as follows: in Sec. 2 we describe the theory and present some details of the method for the system under study, in Sec. 3 our results are shown and a comparison with other theories and simulation results is made, and in Sec. 4 our con- clusions are drawn. The treatment of the hard-sphere gas and the main features of the MC-FSSsimulation method aresummarized respectively inAppendix AandAppendix B. 3 2 Theory Here we consider a fluid of spherical particles interacting via a two-body potential v(r) which is the sum of a singular repulsive hard-sphere contribution and an attractive tail w(r) < 0. The expression for v(r) is then + r < 1 ∞ v(r) = (2)   w(r) r > 1, where the hard-sphere diameter hasbeen set equal equal to unity. As is customary in integral equation theories of fluids, the present approach introduces an approximate closure relation for the direct correlation function c(r) which, once supplemented with the exact Ornstein-Zernike equation involving c(r) and the radial distribution function g(r), yields a closed theory for the thermodynamics and the correlations of the system under study. The basic requirement we want to incorporate in the SCOZA is the consistency between the compressibility and internal energy route to the thermodynamics. According to the compressibility route, the thermodynamics stems from the reduced compressibility χ as determined by the sum rule red 1 χ = , (3) red 1 ρc(k = 0) − where c(k) is the Fourier transform of the direct correlation function and ρ is the number b density of the system. In the internal energy route the key to the thermodynamics is insteadbprovided by the excess internal energy as given by the integral of the interaction weighted by the radial distribution function: +∞ u = 2πρ2 drr2w(r)g(r), (4) Z1 where u is the excess internal energy per unit volume and we have taken into account that g(r) vanishes for r < 1 due to the hard-core repulsion. In the following we will refer to the “excess internal energy” simply as the “internal energy”. If χ and u come from red a unique Helmholtz free energy it is straightforward to find that one must have ∂ 1 ∂2u = ρ , (5) ∂β χred! ∂ρ2 where β = 1/(k T), T being the absolute temperature, and k the Boltzmann constant. B B While this relation is of course satisfied by the exact compressibility and internal energy, this is not the case with those predicted by most integral equation theories. In order to cope with this lack of thermodynamic consistency, we consider the following closure to the Ornstein-Zernike equation: g(r) = 0 r < 1, (6)   c(r) = cHS(r)+K(ρ,β)w(r) r > 1,  4 where c (r) is the direct correlation function of the hard-sphere fluid, and K(ρ,β) is a HS function of the thermodynamic state of the system. In Eq. (6) the approximation clearly lies in the simple form of c(r) outside the repulsive core. The closure above resembles the one used in the approximation known as both the lowest-order gamma-ordered approximation (LOGA) [17] and the optimized random phase approximation (ORPA) [18]. However, while in the LOGA/ORPA one has K(ρ,β) β, in Eq. (6) K(ρ,β) is not ≡ − fixed a priori, but instead must be determined so that the thermodynamic consistency condition (5) is satisfied. This gives rise to a partial differential equation (PDE) for the function K(ρ,β), provided an expression for the hard-sphere direct correlation function c (r) is given. The most popular parameterization for c (r) in the fluid region is due HS HS to Verlet and Weis [19]. Another choice that yields the same thermodynamics as Verlet- Weis, and that we find convenient in view of the calculations performed in this work, is originally due to Waisman [20]. It was subsequently extended analytically by Høye and Stell [21] and explored in some detail by Henderson and coworkers [22]. It amounts to assuming that the function c (r) outside the repulsive core has a one-Yukawa form, so HS that for the hard-sphere system we have: g (r) = 0 r < 1, HS  (7)  c (r) = K exp[−z1(r−1)] r > 1. HS 1 r The Ornstein-Zernike equation supplemented by Eq. (7) can be solved analytically in terms of the amplitude K and the inverse range z of c (r). These can be in turn 1 1 HS determined as a function of the density by requiring, as in the Verlet-Weis parameterization, that both the compressibility and the virial route to the thermodynamics give the Carnahan-Starling equation of state. The basic features of the calculation are recalled in Appendix A. Aconsiderable, althoughpurelytechnical, simplificationintheclosurescheme outlined above based on Eqs. (5), (6) occurs when also the attractive potential w(r) in Eq. (2) is given by a Yukawa function, i.e. when one has exp[ z(r 1)] w(r) = − − , (8) − r z being the inverse range of the potential. From Eq. (7) it is then immediately seen that Eq. (6) becomes g(r) = 0 r < 1,  (9)  c(r) = K exp[−z1(r −1)] +K exp[−z2(r −1)] r > 1, 1 2 r r where K andz are the quantities referred to as K and z in Eq. (6), (8), and K , z are 2 2 1 1 known function of the density. It is now possible to take advantage of the fact that for theOrnstein-Zernike equationsupplemented by theclosure (9)extensive analyticalresults 5 have been determined [23, 24, 14]. If both K and K are given, as in the LOGA/ORPA, 1 2 this enables one to solve Eq. (9) altogether [25, 26, 27]. More generally, irrespective of the form of K and K , a prescription can be found to determine the reduced compressibility 1 2 χ as a function of the density ρ and the internal energy per unit volume u, which can red be used in Eq. (5) to obtain a closed PDE. A similar procedure for the same potential considered here was adopted in a previous work [15], where however the hard-sphere contribution to the direct correlation function c (r) outside the core was not taken into HS account, so that c(r) was given by a simple one-Yukawa tail. This further simplifies the theory, but implies that the description of the hard-sphere fluid coincides with that of the PY approximation, which as is well known is not very satisfactory at high density. This defect becomes more and more severe as the range of the attractive interaction decreases, and can considerably affect the phase diagram predicted by the theory, unless some more- or-less ad hoc procedure is adopted to correct the hard-sphere thermodynamics. In order to incorporate a better treatment of the hard-sphere fluid into the theory one can turn to the two-Yukawa form for c(r) of Eq. (9), whose use in the consistency condition (5) we are now going to illustrate in some detail. In the following we will exploit the results determined in Refs.[23, 24, 14], which will be respectively referred to as I, II, III. Let us introduce the packing fraction ξ = πρ/6 and the quantity 1 f = (1 ξ) , (10) − sχred which is the square root of the quantity referred to as A in I, II, III. Eq. (5) becomes 2f ∂f ∂u ∂2u = ρ . (11) (1 ξ)2 ∂u! ∂β! ∂ρ2! − ρ ρ β To obtain a PDE for u we need to express f as a function of ρ and u in Eq. (11). From Eq. (II.14) it is found that f can be written as (z2 z2)+4√q(γ γ ) z2 z2 γ γ (γ γ ) f = 1 − 2 2 − 1 1 − 2 1 2 2 − 1 , (12) − 4[(z /z )γ (z /z )γ ] − z z [(z /z )γ (z /z )γ ]2 1 2 2 − 2 1 1 1 2 1 2 2− 2 1 1 where we have set (1+2ξ)2 q = . (13) (1 ξ)2 − The quantities γ and γ are given by Eq. (II.5) 1 2 U 1 γ = 2 √q , (14) 1 − − U 0 W 1 γ = 2 √q . (15) 2 − − W 0 The ratios U /U and W /W depend on the integrals 1 0 1 0 +∞ I = 4πρ drrexp[ z (r 1)]g(r) (i = 1,2). (16) i i − − Z1 6 ¿From Eq. (I.35) it is found in fact W 4+2z z2 τ I 1 1 = 2 − 2 2 2 − , (17) W 2(2+z ) σ I 1 0 2 2 2 − and the corresponding relation with W /W replaced by U /U and the index 2 changed 1 0 1 0 to 1. The quantities τ and σ depend only on z and are given by Eq. (I.34): i i i 1 z 2 i σ = − +exp( z ) , (18) i i 2z z +2 − i (cid:20) i (cid:21) 1 z2 +2z 4 τ = i i − +exp( z ) , (19) i 2zi "4+2zi −zi2 − i # with i = 1,2. From the expression of the potential (8) it is readily seen that I is directly 2 related to the internal energy per unit volume u given by Eq. (4): 1 u = ρI . (20) 2 −2 Eqs. (15), (17), (20) allow then to express γ explicitly as a function of ρ and u: 2 4+2z z2 2τ u+ρ γ = 2 √q 2 − 2 2 . (21) 2 − − 2(2+z ) 2σ u+ρ 2 2 We now need γ as a function of ρ and u. This is less straightforward than for γ , since 1 2 the integral I does not have any direct thermodynamic meaning, the exponential in I 1 1 being related to the tail of the direct correlation function of the hard-sphere gas. We have then to make use of some further results determined in I–III. ¿From Eq. (I.36) it is found that the amplitudes K , K of the Yukawa functions in the closure (9) can be expressed 1 2 in terms of the above introduced quantities U , U , W , W . One has 0 1 0 1 2(z +2)2σ2 U 2 K = 1 1 U 1 α , (22) 1 3ξz2 0 U − 1 1 (cid:20) 0 (cid:21) where α is given by Eq. (I.37): 1 (4+2z z2)τ α = 1 − 1 1 , (23) 1 2(2+z )σ 1 1 and the corresponding equations with the index 1 replaced by 2 and U , U replaced by 0 1 W , W . Let us now introduce the quantities x, y given by 0 1 z2 x = √q 1 , (24) − 4γ 1 z2 y = √q 2 . (25) − 4γ 2 7 ¿From Eq. (III.30) one has 4 U = p(√q x)2, (26) 0 z2 − 1 4 W = s(√q y)2. (27) 0 z2 − 2 where p and s must satisfy Eq. (II.39) (in the notation of II one has x u /u , y q1 q0 ≡ ≡ w /w , p u , s w ): q1 q0 q0 q0 ≡ ≡ 4s p+s+ (y x)2 = 1z2 x2, z2 z2 − 4 1 −  p+s− z124−pz22(y −x)2 = 41z22 −y2. (28) 1 − 2 Eq. (28) is readily solved for p and s to give z2 z2 p = 1 − 2 4z2(y x)2 16y2(y x)2 (z2 z2) z2 z2 +4(y2 x2) , (29) −64(y x)4 2 − − − − 1 − 2 1 − 2 − − n h io and the expression for s is obtained by exchanging z , z and x, y in the r.h.s. of Eq. (29). 1 2 If Eqs. (14), (24), (26), (29) are used in Eq. (22) we finally obtain 2 4(2 √q α )(√q x) z2 4z2(y x)2 16y2(y x)2 − − 1 − − 1 2 − − − h i n 384ξz4 (z2 z2) z2 z2 +4(y2 x2) = 1 K (y x)4, (30) − 1 − 2 1 − 2 − −(z +2)2(z2 z2)σ2 1 − h io 1 1 − 2 1 and a similar equation obtained by exchanging the indices 1 and 2 and the quantities x, y. We recall that in Eq. (30) K , z , σ , and α are known functions of the density ρ which 1 1 1 1 refer to the hard-sphere system. For given values of ρ and u, Eqs. (21), (25) allow one to determine y. Eq. (30) can then be solved numerically with respect to x to obtain γ 1 via Eq. (24). This solves the problem of determining γ in terms of ρ and u. The partial 1 derivative (∂f/∂u) that appears in Eq. (11) can then be determined as ρ ∂f ∂f ∂γ ∂f ∂γ 1 2 = + , (31) ∂u!ρ ∂γ1!ρ ∂u !ρ ∂γ2!ρ ∂u !ρ where (∂γ /∂u) is calculated explicitly by Eq. (21), while (∂γ /∂u) must be determined 2 ρ 1 ρ as the derivative of the function implicitly defined by Eq. (30). If we write Eq. (30) as F(x,y,ρ) = 0, it is found straightforwardly that Eq. (11) takes the form ∂u ∂2u B(ρ,u) = C(ρ,u) , (32) ∂β ∂ρ2 where the functions B(ρ,u) and C(ρ,u) are given by the following expressions: 2f ∂γ ∂f ∂F ∂x ∂f ∂F ∂y 2 B(ρ,u) = , (33) (1 ξ)2 ∂u "∂γ2 ∂x ∂γ1 − ∂γ1 ∂y ∂γ2# − ∂F ∂x C(ρ,u) = ρ . (34) ∂x ∂γ 1 8 All the partial derivatives in Eqs. (33), (34) are calculated at constant ρ and can be determined by Eqs. (12), (21), (24), (25), (30). The resulting expressions are then evaluated as a function of ρ and u via the procedure described above. The same procedure also allows one to determine the reduced compressibility as 1/χ = f2/(1 ξ)2 once f has been red − obtained from Eq. (12). The PDE (32) is a non-linear diffusion equation that must be integrated numerically. To prevent the occurrence of any numerical instability, especially in the critical and sub-critical region, we have adopted an implicit finite-differences algo- rithm [28] tailored to equations that, although globally non-linear, depend on the partial derivatives of the unknown function in a linear fashion like Eq. (32). The integration with respect to β starts at β = 0 and goes down to lower and lower temperatures. Before each integration step Eq. (30) is solved numerically and the coefficients B(ρ,u), C(ρ,u) are determined. The density ρ ranges in a finite interval (0,ρ ), whose high-density bound- 0 ary has been typically set at ρ = 1. The initial condition can be determined by taking 0 into account that at β = 0 the radial distribution function coincides with that of the hard-sphere gas. From Eqs. (4) and (8) one has then +∞ u(ρ,β = 0) = 2πρ2 drrexp[ z (r 1)]g (r) for every ρ, (35) 2 HS − − − Z0 where g (r) is obtained in the present scheme by the closure (7). For such a closure, HS as shown in Appendix A, both U and U in Eq. (14) can be determined analytically as 0 1 a function of ρ, thus providing γ (ρ) at β = 0. This allows one to obtain u in Eq. (35) 1 analytically as well: in fact, one can solve Eq. (12) for γ as a function of γ , f, and ρ, 2 1 where f is readily obtained by using the Carnahan-Starling expression of χ in Eq. (10). red Once γ is known, Eq. (21) is solved with respect to u. It must be noted that solving 2 Eq. (12) for γ gives two branches, so attention must be paid in order to single out the 2 branch that actually corresponds to the physical solution. We also need two boundary conditions at ρ = 0 and ρ = ρ . ¿From Eq. (4) one has immediately 0 u(ρ = 0,β) = 0 for every β. (36) At high density we instead make use of the so-called high-temperature approximation (HTA),accordingtowhich theexcess Helmholtzfreeenergyperunitvolumeisdetermined via Eq. (35) for every temperature. In the fluid region of the phase diagram this of course is not exact unless β = 0, but it becomes more and more accurate as the density of the system is increased [29], so we expect that for a given sweep along the β-axis the results will not differ appreciably from what would be obtained using an hypothetical exact boundary condition, provided the boundary ρ is located at sufficiently high density. We 0 used the HTA at ρ = ρ for the reduced compressibility. This yields via Eq. (5) the 0 boundary condition ∂2u ∂2u (ρ ,β) = (ρ ,β=0) for every β. (37) ∂ρ2 0 ∂ρ2 0 We have checked that the output of the numerical integration of Eq. (32) is quite in- sensitive to the specific choice of the high-density boundary condition. Moreover, for 9 ρ 1 moving the boundary condition to higher densities also leaves the results un- 0 ≃ affected. Eq. (10) shows that to be physically meaningful, the quantity f has to be non-negative. On the other hand, below the critical temperature the solution of Eq. (32) does not satisfy this condition along the whole density interval (0,ρ ), but only outside a 0 certain temperature-dependent region (ρ (β),ρ (β)). For ρ = ρ (β) or ρ = ρ (β) the s1 s2 s1 s2 quantity f vanishes, and consequently the compressibility diverges. As β changes, ρ (β) s1 and ρ (β) give then respectively the low- and the high- density branch of the spinodal s2 curve predicted by the theory. The fact that f becomes negative for ρ (β) < ρ < ρ (β) s1 s2 not only implies that the theory behaves unphysically in this interval, but it also gives rise to an analytical instability which would make the numerical integration of the PDE (32) impossible, if one tried to determine the solution over the whole interval (0,ρ ) even below 0 thecriticaltemperature. Therefore, theregionboundedbythespinodalhasbeenexcluded from the integration of Eq. (32). Specifically, as soon as it is found that f changes sign, so that for a certain density ρ one has f(ρ,β) < 0, the integration is restricted to the interval (0,ρ ∆ρ) or (ρ + ∆ρ,ρ ) respectively for ρ < ρ or ρ > ρ , where ∆ρ is the 0 c c − spacing of the density grid. Weithin the preecision of the numerical discretization, one has ρ = ρ ∆eρ (or ρ = ρe+∆ρ) and the further bounedary conditeions s1 s2 − e ue(ρsi,β) = uS(ρsi) i = 1,2, β > βc, (38) where β is the critical inverse temperature and u (ρ) is the value of the internal energy c S per unit volume when the compressibility at density ρ diverges. This can be determined by setting f = 0 in Eq. (12) and solving for γ as a function of ρ and γ . If Eqs. (24) 1 2 and (25) are substituted into Eq. (30), an equation for γ is obtained that allows one to 2 determine the value of γ when 1/χ = 0 for a certain ρ. Solving Eq. (21) with respect 2 red to u then yields u (ρ). S Once the internal energy per unit volume u has been determined from Eq. (32), the pressure P and the chemical potential µ are obtained by integration with respect to β via the relations ∂(βP)/∂β = u + ρ∂u/∂ρ, ∂(βµ)/∂β = ∂u/∂ρ. Thanks to the self- − consistency of the theory, this route to the thermodynamics is equivalent to integrating the inverse compressibility with respect to ρ, but it does not require one to circumvent the forbidden region in order to reach the high-density branch of the subcritical isotherms. 3 Results The numerical integration of the PDE (32) with the initial condition (35) and the boundary conditions (36)–(38) has been performed on a density grid with ∆ρ = 10−3–10−4. At the beginning of the integration the temperature step ∆β was usually set at ∆β = 2 10−5–10−5. As the temperature approaches its critical value, ∆β can be further de- × creased if one wishes to get very close to the critical point, and then gradually expanded back. The integration was usually carried down to β 2.4β . The inverse range pa- c ≃ rameter of the attractive tail in Eq. (8) has been set at z = 1.8. For this value of z several simulations [30, 31, 32] and theoretical [6, 26, 27] predictions have already been reported in the literature. Fig. 1 shows the SCOZA results for the compressibility factor 10