Davide ; Pini

Alberto Parola

Jader Colombo

Luciano Reatto

An investigation of the SCOZA for narrow square-well potentials and in the sticky limit

Keywords: square-well potential, sticky hard spheres, liquid-vapor transition, SCOZA

Narrow attractive interactions are ubiquitous in the realm of colloidal dispersions. A most relevant instance is represented by depletion-induced forces in mixtures of highly asymmetric hard spheres, a field that has received a fundamental contribution by Bob Evans [1]. Although these investigations have shown that the structure of effective depletion interactions is in general more complex than just a single attractive tail, a narrow attractive contribution at short distance is still to be expected in most cases, and its relevance to the modelization of both the equilibrium and non equilibrium phase behavior of colloidal systems has long been acknowledged.

In view of this, it is somewhat disappointing that the difficulties experienced by most liquid-state theories when phase separation occurs, the foremost of which is lack of convergence in a (not necessarily small) neighborhood of the critical region, only become worse when the range of the attractive potential decreases to a small fraction of the particle size.

However, it may be argued that this is not a fundamental difficulty, inasmuch as this class of interactions can be dealt with by the sticky hard-sphere (SHS) model [2]. This is obtained from the prototypical attractive interaction, i.e., the hard-core plus square-well (SW) potential v SW (r) given by

v SW (r) =          +∞ r < σ , -ǫ σ ≤ r ≤ (1 + δ)σ , 0 r > (1 + δ)σ , (1)
by taking the limits δ → 0, ǫ → ∞ so that the second virial coefficient B 2 stays finite and non vanishing:

βv SHS (r) =            +∞ r < σ , -ln 1 12τ (1 + 1 δ ) σ ≤ r ≤ (1 + δ)σ δ → 0 , 0 r > (1 + δ)σ . (2)
The stickiness parameter τ is related to the second virial coefficient of the SHS model B SHS 2 by the relation

τ = B HS 2 4 (B HS 2 -B SHS 2 ) , (3)
where B HS 2 = 2πσ 3 /3 is the second virial coefficient of the hard-sphere fluid. As shown by Baxter [2], this model be can solved exactly within the Percus-Yevick (PY) integral equation [START_REF] Hansen | Theory of Simple Liquids[END_REF]. Because of the way the limit (2) is taken, the solution does not depend on the temperature, but rather on the stickiness τ or, equivalently, on B SHS

2

. The connection between the SHS model and an actual fluid at a certain temperature T can be made by a generalized principle of corresponding states that, although not rigorous, has been found to work remarkably well. According to this principle, usually referred to as Noro-Frenkel (NF)

scaling [START_REF] Noro | [END_REF], a pair potential v(r) consisting of a hard-core repulsion with diameter σ and an attractive tail w(r) can be mapped into the SW potential (1), provided that the width δ of the attractive well is determined so that its second virial coefficient B 2 (T ) will be the same as that of the original interaction at the same temperature T , assuming that in both cases T is measured in units of the well depth. On the other hand, when δ is small, by following the same procedure the SW potential can in turn be mapped into the SHS model. In this way, NF scaling can be combined with Baxter analytical solution of the SHS model to provide a unified description of fluids with narrow attractive interactions. For instance, the critical temperature T c of such a fluid can be evaluated by that at which its B 2 (T ) is equal to the critical value of B SHS

2

.

Although the SHS model has a perfectly regular solution within the PY approximation and has also been the subject of intensive simulation studies [5][6][7][8], the singularity embedded into its very definition should instigate a mild feeling of suspicion. In fact, twenty-three years after the appearance of the model, Stell [9] discovered that the latter from a rigorous statistical mechanical point of view is thermodynamically unstable whenever the number of particles is greater than eleven. However, the state of affairs is less serious than it might seem: the instability arises from a number of rather fancy particle configurations, that are readily destroyed by the slightest amount of polidispersity in the system. Therefore, the SHS model and its approximate analytical solution via the PY equation [2] still remain a valuable tool for the interpretation of the phase diagram of colloidal dispersions subject to short-range interactions.

There is, however, a significant limitation of this solution that negatively affects its predictive power, i.e., its lack of thermodynamic consistency. This is particularly evident when the liquid-vapor transition of the model is considered: the coexistence curves determined by the compressibility and internal energy routes described in Sec. II differ markedly from each other, the discrepancy between the critical densities amounting to nearly a factor three, and neither of them agrees quantitatively with the simulation results [7]. In this respect, it is also worth recalling that the PY approximation upon which the analytical solution of the SHS is based, is expected to become less and less accurate as the density is increased, while on the other hand the critical density increases as the range of the interaction decreases, and for δ → 0 the limiting value predicted by the simulations [7,10] is considerably larger than that expected for Lennard-Jones like fluids.

Recently, a different approach has emerged, which is based upon enforcing consistency between the aforementioned energy and compressibility routes to the thermodynamics, namely, the self-consistent Ornstein-Zernike approximation (SCOZA) [11][12][13][14][15]. This has been applied especially to the case in which v(r) consists of a hard core and a Yukawa (HCY) attractive tail:

v HCY (r) =    +∞ r < σ ,
-ǫ σ r e -z(r/σ-1) r > σ ,

where z is the inverse-range parameter of the interaction. For the HCY potential, the SCOZA has shown a remarkable robustness even for very short-range tails, while still remaining accurate [14,15]. On the other hand, one might argue that accuracy cannot be expected to persist for attractions of arbitrarily short range, since after all, as recalled in Sec. II, the SCOZA assumes that the direct correlation function c(r) of the fluid is linear in the attractive part w(r) of the interaction, and its low-density expansion does not even give back the correct expression of B 2 (T ). In [16] it was suggested that for narrow attractive interactions it could be more appropriate to have a c(r) which depends linearly not on w(r), but rather on its Mayer function e -βw(r) -1, where β = 1/(k B T ) is the inverse temperature.

If this prescription were adopted in the SCOZA, one would obtain a "non-linear SCOZA" bound to become exact at low density. Now, the SW potential is quite special in this respect, because both w(r) and its Mayer function are constant in the well region and vanish elsewhere, and therefore they have the same functional form as far as they dependence on r is concerned, the only difference being their (state-dependent) amplitude. On the other hand, in the SCOZA this amplitude is determined by an exact thermodynamic consistency condition. As a consequence, for the SW potential the usual SCOZA must coincide with its non-linear counterpart, and hence predict the exact B 2 (T ) at low density, unlike for a generic tail potential. It is then tempting to investigate how the SCOZA will behave for narrow SW potentials: will the quite special status of this potential enable it to recover a finite B 2 (T c ) in the limit δ → 0? This possibility is quite enticing, because it would allow one to fix the lack of thermodynamic consistency that is the main defect of the PY solution of the SHS fluid [2], possibly obtaining a better agreement with the simulation results for the phase diagram [7].

With this in mind, we have performed a study of the SCOZA for narrow SW interactions.

Such a task requires a fully numerical solution of the SCOZA, because one cannot take advantage of the wealth of analytical results available for the HCY potentials [17]. In fact, this problem was already tackled in [18], where the numerical solution of the SCOZA for a generic tail interaction w(r) was worked out and applied to the SW potential. The investigation performed there showed that on narrowing the well, the results become more and more dependent on the choice of the upper boundary ρ 0 of the density interval where the SCOZA is solved, at which a high-density boundary condition needs to be specified. In order to get rid of such a dependence, one needs to move ρ 0 to higher and higher values, but eventually this is prevented by a lack of convergence of the numerical algorithm. On the basis of the conclusions of that work, we have then reconsidered the numerical solution of the SCOZA paying special attention to the problem of obtaining convergence at high density, and developed a numerical algorithm that remains convergent up to ρ 0 ≃ 1.4σ -3 , which is near the density at close packing. By this, we have managed to obtain results for wells down to δ = 0.1, which however is not small enough to represent reliably the behavior of the SCOZA in the limit δ → 0.

We have then resorted to a different strategy, namely, we have taken the sticky limit (2) of the SW potential from the outset, and imposed thermodynamic consistency on the top of this. Again, the results turn out to be strongly dependent on the high-density boundary condition, but in this case the simpler form of the SCOZA equation allows one to circumvent this problem by resorting to a solution technique which requires only the behavior at low density to be specified. Unfortunately, according to this solution the SCOZA fails to have a liquid-vapor transition in the sticky limit (2).

However disappointing, the full picture appears to be consistent enough to indicate that the SCOZA critical temperature will decrease too rapidly to give a finite value of B 2 (T c ) for δ → 0. On the other hand, the approximation which is introduced in the SCOZA, namely, that the contribution to the direct correlation function c(r) due to the tail interaction w(r)

has the same range as w(r) itself, and therefore vanishes wherever w(r) does; is made also in the PY equation, which moreover is thermodynamically inconsistent. This suggests The plan of the paper is as follows: in Sec. II, the algorithm employed for the fully numerical solution in the case of a generic attractive tail w(r) is described; in Sec. III our results for SW potentials of several amplitudes δ are presented, and compared with numerical simulations; the sticky limit of the SCOZA is considered in Sec. IV; our conclusions are presented in Sec. V. At the end of the paper, three Appendixes deal with some technical details pertinent to Sec. II: Appendix A presents the finite-difference algorithm used to integrate the SCOZA partial differential equation (PDE); Appendix B concerns the minimization algorithm which we used to optimize c(r) inside the repulsive core of the interaction; Appendix C describes the technique adopted for the treatment of the hard-sphere part of the correlations at high density.

II. THE SCOZA FOR A GENERIC ATTRACTIVE TAIL POTENTIAL

We are interested in a simple model fluid of particles interacting via a two-body, spherically symmetric potential v(r) that consists of a hard core with diameter σ and an attractive tail w(r):

v(r) =    +∞ r < σ , w(r) r > σ . (5)
As is customary in integral-equation theories, the SCOZA is based on an approximate closure of the exact Ornstein-Zernike (OZ) equation that relates the two-body correlation function h(r) and the direct correlation function c(r):

h(r) = c(r) + ρ d 3 r ′ c(|r -r ′ |) h(r ′ ) , (6)
where ρ ≡ N/V is the number density of the fluid, V being the volume, and N the particle number. The SCOZA closure reads

   h(r) = -1 r < σ , c(r) = c HS (r) -K(ρ, β) w(r) r > σ , (7)
where K(ρ, β) is a state-dependent function of the density ρ and the inverse temperature β = 1/(k B T ), and c HS (r) is the direct correlation function of the hard-sphere fluid. The latter is assumed to be known, e.g. by the Verlet-Weis [19] or the Waisman [20] parametrization. Here the Waisman parametrization has been used. Equation ( 7) is similar to the well-known optimized random phase approximation (ORPA) [START_REF] Hansen | Theory of Simple Liquids[END_REF]: in particular, the requirement on h(r)

for r < σ, usually referred to as the core condition, is exact because of the hard core of the potential (5), while the expression of c(r) is clearly an approximation, since the off-core contribution to c(r) is assumed to depend linearly on the interaction w(r), hence having the same range as w(r) itself -the so-called OZ ansatz. However, in the ORPA the amplitude of w(r) is set to K = β, while in the SCOZA K is a priori unknown, and must be determined so as to achieve consistency between the compressibility and internal energy routes to the thermodynamics. This requirement is expressed by the following condition on the isothermal reduced compressibility χ red and the excess internal energy per particle U:

∂ ∂β

1 χ red = ρ ∂ 2 (ρU) ∂ρ 2 , ( 8
)
where it is understood that χ red and U are determined respectively by the compressibility and internal energy routes, namely:

1 χ red = 1 -ρ d 3 r c(r) , (9)
U = 1 2 ρ d 3 r [h(r) + 1]w(r) . (10)
Equation ( 8) would in fact be a trivial identity in a hypothetical exact description of the system, but this is not true anymore once an approximate closure such as Eq. ( 7) is introduced. Together with Eqs. ( 7), ( 9), (10), Eq. ( 8) yields a PDE for K(ρ, β), that must be solved numerically. As in [12,18], Eq. ( 8) is rewritten using the internal energy per unit volume u = ρU as the unknown function:

D(ρ, u) ∂u ∂β = ρ ∂ 2 u ∂ρ 2 , (11)
where we have set

D(ρ, u) = ∂ ∂u 1 χ red ρ . ( 12
)
The finite-difference algorithm used for the numerical integration of Eq. ( 11) has been described in detail in Appendix A. Here we observe that, irrespective of the specific discretization of the PDE, a necessary step of the integration scheme consists in solving the OZ equation ( 6) supplemented with the closure (7). This amounts to solving the ORPA for an effective inverse temperature β eff = K. In several previous applications of the SCOZA, w(r) was chosen as a Yukawa tail [11][12][13][14], a linear combination of Yukawa tails [21,22], or other similar forms such as the Sogami-ESE [23] or related potentials [24]. If c HS (r) for r > σ is assumed to be vanishing as in the PY equation for the hard-sphere fluid [START_REF] Hansen | Theory of Simple Liquids[END_REF], or is also given by a Yukawa tail as in the Waisman parametrization, these expressions of w(r) allow for an almost fully analytical solution of Eqs. ( 6), (7). Although the integration of the PDE must always be carried out numerically, this leads nevertheless to a considerable simplification of the whole procedure.

For the SW potential considered here, on the other hand, it is necessary to solve also Eqs. ( 6), (7) in a fully numerical fashion. Moreover, as already found in [18] and recalled in the Introduction, for short-range wells the solution of the PDE is very sensitive to the position of the upper boundary ρ 0 of the density interval, unless this is moved to very high values. Therefore, it is of paramount importance that the numerical algorithm employed remains convergent in this regime. In [18] Eqs. ( 6), (7) were solved by the Labik, Malijevsky, and Vonka algorithm [25]. However, in that work it was found that the algorithm failed to converge at high density and low temperature, so that for the shortest-range wells investigated there, namely for 0.25 ≤ δ ≤ 0.5, the high-density boundary ρ 0 was not moved beyond ρ 0 σ 3 = 1.15. Therefore, here we have adopted a different method, which was originally developed in [26,27] for the ORPA, and again applied to the ORPA for narrow SW and square-shoulder potentials in [28]. This is based on the well-known property [START_REF] Hansen | Theory of Simple Liquids[END_REF] that imposing the core condition in the ORPA is equivalent to requiring the Helmholtz free energy A RPA given by the random phase approximation (RPA) [START_REF] Hansen | Theory of Simple Liquids[END_REF] to be stationary with respect to variations inside the hard core of the quantity φ(r) defined as

φ(r) ≡ c(r) -c HS (r) . (13)
We recall that outside the core one has φ(r) = -βw(r) in both the RPA and the ORPA.

Specifically, the expression of

A RPA is - βA RPA N = - βA HS N + f HTA + 1 2 f Ring , (14)
where A HS is the Helmholtz free energy of the hard-sphere fluid, and f HTA and f Ring are given by where h HS (r) and S HS (k) are respectively the two-body correlation function and the structure factor of the hard-sphere fluid, and the hats denote Fourier transforms. In the equations above, f HTA refers to to the high-temperature approximation, whereby the two-body correlation function h(r) is replaced by that of the hard-sphere fluid, while f Ring refers to the sum of the ring diagrams in the expansion of A RPA in powers of φ(r) [START_REF] Hansen | Theory of Simple Liquids[END_REF]. If f Ring is regarded as a functional of φ(r) inside the hard core, one finds

f HTA = 1 2 ρ d 3 r φ(r)[1 + h HS (r)] (15)
f Ring = - 1 ρ d 3 k (2π) 3 ρS HS (k) φ(k) + ln 1 -ρS HS (k) φ(k) , (16) 8
δf Ring δφ(r) = ρ ∆h(r) , (17)
δ 2 f Ring δφ(r)δφ(r ′ ) = ρ d 3 k (2π) 3 e ik•(r-r ′ ) S 2 (k) ( 18
)
where we have set ∆h(r) = h(r) -h HS (r), and S(k) is the RPA structure factor given by

S(k) = S HS (k) 1 -ρS HS (k) φ(k) . ( 19
)
Note that, since h HS (r) satisfies exactly the core condition, f HTA does not depend on φ(r)

for r < σ, and therefore does not contribute to the functional derivatives of A RPA for r < σ.

Equation (17) shows that the core condition is satisfied if and only if the RPA free energy is stationary with respect to variations of φ(r) inside the core, and Eq. (18) shows that f Ring is convex, so that this stationary point is unique and corresponds to a minimum of f Ring .

As a consequence, one can solve the ORPA by minimizing f Ring numerically with respect to φ(r) in the core region [26]. Since Eqs. (17), (18) hold irrespective of the form of φ(r)

outside the core, this procedure can be equally well applied to the present case where (see Eq. ( 7)) φ(r) for r > σ is given by -K(ρ, β)w(r). Therefore, in order to solve the OZ equation ( 6) with the closure (7), we can, for any given K, minimize the functional f Ring of Eq. ( 16) with φ(r) = -Kw(r) for r > σ. The difference with respect to the ORPA consists in the fact that in the ORPA, the minimum of the functional gives the Helmholtz free energy of the system as predicted by the energy route, while this is not true anymore for a generic K. In fact, one finds

∂f Ring ∂K = -ρ d 3 r ∆h(r)w(r) = -2 (U -U HTA ) , ( 20
)
where U is given by Eq. ( 10), and U HTA is the corresponding quantity in the high-temperature approximation: For K = β, Eq. ( 20) coincides with the usual thermodynamic relation ∂(βA/N)/∂β = U with A given by Eq. ( 14), so that A is indeed the Helmholtz free energy obtained by integrating the internal energy U with respect to β. On the other hand, this is not the case if K has a nontrivial dependence on β. Therefore, unlike in the ORPA and the RPA, in the SCOZA Eqs. ( 14), ( 15), ( 16) do not give the Helmholtz free energy. Nevertheless, Eq. ( 20)

U HTA = 1 2 ρ d 3 r [h HS (r) + 1]w(r) . (21) 9
turns out to be very useful to overcome a difficulty related to the numerical solution of the SCOZA. Specifically, the calculation of the quantity D(ρ, u) defined in Eq. ( 12) requires that Eqs. ( 6), ( 7) must be solved at fixed ρ and U, rather than at fixed ρ and K, i.e., for any given density one needs to find the value of K such that the internal energy per particle takes a given value U = Ū . One way to do this would be to change K by trial and error, solve Eqs. ( 6), ( 7) with respect to φ(r) inside the core for each K, and find the corresponding U by Eq. ( 10), until the condition U = Ū is met up to a prescribed accuracy. However, we found that one can solve Eqs. ( 6), (7) simultaneously with respect to φ(r) and K in a single optimization run. To this end, let us observe that one has

∂ 2 f Ring ∂K 2 = ρ 2 d 3 k (2π) 3 [S(k)w(k)] 2 , (22)
where S(k) is given by Eq. ( 19). We now introduce the Legendre transform of f Ring

S = f Ring -K ∂f Ring ∂K = f Ring + 2K(U -U HTA ) . (23)
In the case of the ORPA in which K = β, S coincides with twice the excess entropy per particle with respect to the hard-sphere fluid. Equations ( 20), (22) imply that, if 2(U -U HTA ) is fixed at some given value ∆U, S is a convex function of K such that

∂S ∂K = -ρ d 3 r ∆h(r)w(r) + ∆U , (24)
so that its minimum is obtained for that K which, when inserted into Eq. ( 7), will yield ∆U via Eqs. ( 10), ( 21) for the internal energy. Therefore, in order to solve the OZ equation ( 6)

with the closure (7) for a fixed density ρ and internal energy U, it is sufficient to minimize S, regarded as a function of K and a functional of φ(r) for r < σ. The numerical procedure adopted for the minimization is described in detail in Appendix B.

A difficulty related to the aforementioned request that Eqs. ( 6), ( 7) must be solved also at very high density, is that in this regime the structure factor of the fluid develops very high and narrow, "pseudo-Bragg" peaks. These do not have a direct physical meaning, since they occur in a density range where the stable phase would actually be the solid one. Nevertheless, they must be taken into account when solving Eqs. ( 6), ( 7) numerically. In fact, this problem appears already for the simple hard-sphere fluid: if one starts from the (analytical) ĉHS (k)

given by the Waisman parametrization, finds ĥHS (k) by the OZ equation ( 6) and obtains h HS (r) by performing directly the inverse Fourier transform of ĥHS (k), it is found that for ρ * 1.2, the core condition is very poorly satisfied. The reason for this is that many of the peaks of ĥHS (k) will simply be missed by the discrete sampling of the wave vector, unless the number of point used in the transform is so large (about 2 18 ), that the numerical computation becomes unpractical. In order to cope with this problem, the contribution of these peaks to the inverse Fourier transform has been determined analytically by the method described in Appendix C.

The initial condition required by the PDE ( 11) at β = 0 and the boundary conditions at ρ = 0 and on the spinodal curve, i.e., the locus of diverging compressibility below the critical temperature, have been detailed elsewhere [12]. The high-density boundary condition is more problematic in this context, since this is not known exactly while, as stated in Sec. I, the results for narrow SW are quite sensitive to both the position of the high-density boundary ρ 0 , and the kind of approximation used for u at ρ = ρ 0 . This point will be considered in more detail in Sec. III; here we remark that in this investigation the high-density boundary condition has not been imposed on the second derivative of u with respect to the density ∂ 2 u/∂ρ 2 as in a number of previous works [12][13][14]18], but rather directly on u, because we found that the results obtained in this way were less sensitive to the choice of u(ρ 0 , β).

The vapor-liquid coexistence curve is obtained by imposing the conditions of thermodynamic equilibrium, i.e., by finding the densities ρ v , ρ l which give the same pressure P and the same chemical potential µ on the liquid and vapor branches of the sub-critical isotherm. P and µ are determined by integrating with respect to β the exact identities ∂(βP )/∂β = -u + ρ∂u/∂ρ, ∂(βµ)/∂β = ∂u/∂ρ, which are fulfilled also by the SCOZA because of the consistency condition (8).

Concerning the parameters used in the numerical calculation, we typically set the density step at ∆ρ * = 10 -3 , and the initial inverse-temperature step at ∆β * 0 = 10 -3 . In order to locate the critical point accurately, ∆β is decreased as the critical point is approached, and then gradually expanded back to ∆β 0 below the critical temperature. We found that the results are not very sensitive to the density or temperature mesh, e.g., even a substantial coarsening of the mesh such as ∆ρ * = 10 -2 , ∆β * 0 = 10 -2 leaves the results nearly unaffected, except for the coexistence curve in the neighborhood of the critical point. In this case, a very high accuracy in the numerical calculation is required to fulfill the conditions of thermodynamic equilibrium, so that a coarse mesh results in a large gap at the top of the coexistence curve. On the other hand, the numerical output turns out to be much more sensitive to the discretization used in the numerical Fourier transform. We then used a rather small step in real space, ∆r * = 5 × 10 -4 , and a large number of mesh points, N = 2 15 . We employed the fast Fourier transform routines developed by Takuya Ooura at the Research Institute of Mathematical Sciences, Kyoto University, and made freely available online.

Whenever the function to be transformed is discontinuous, such as ∆h(r) or ∆c(r), it is best to subtract off the discontinuity, and add the slowly decaying transform of the discontinuous part determined analytically at the end. Here we implemented this procedure on both the function and its derivative, so that the numerical transform is done on a smooth function that is not only continuous, but also differentiable. If Df (r 0 ), Df ′ (r 0 ) denote the discontinuities respectively of f and its derivative f ′ at r = r 0 , i.e., Df (r 0 ) = lim r→r - 0 f (r) -

lim r→r + 0 f (r), Df ′ (r 0 ) = lim r→r - 0 f ′ (r) -lim r→r + 0 f ′ (r), the function f that is transformed numerically is given by f(r) = f (r) -[Df (r 0 ) + (r -r 0 )Df ′ (r 0 )] [1 -Θ(r -r 0 )] , (25)
where Θ(r) is the Heaviside step function defined by Θ(r) = 1 for r > 0, Θ(r) = 0 for r < 0.

For the SW potential studied here, the discontinuities of ∆h(r) and ∆c(r) are located at

r 0 = σ and r 0 = σ(1 + δ).
The accuracy within which the core condition is fulfilled can be determined a posteriori by calculating the radial distribution function g(r) ≡ h(r)+1 inside the repulsive core r < σ.

Typically, we found that for T ≃ T c and high densities such as ρ ≃ 0.9σ -3 , the order of magnitude of g(r) is g(r) ∼ 10 -8 -10 -7 for 0.1σ r < σ and g(r) ∼ 10 -6 -10 -3 for r 0.1σ.

III. RESULTS

We have employed the numerical solution of the SCOZA described in Sec. II to study SW potentials of several widths δ. For large δ, our results are basically identical to those obtained in [18]. As can be inferred from the comparison between Tab. I shown further in this Section and Tab. I of [18], the critical points for δ ≥ 0.5 differ at most by ∼ 0.2%. Here we discuss in detail our data for δ ≤ 0.25, which is the smallest value studied in [18]. In this regime, the liquid-vapor transition considered here is actually expected to be metastable with respect to freezing [29].

As already found in [18] and anticipated in Secs. I, II, we also have observed a strong influence on the results of the high-density boundary condition, that becomes more and more important as δ decreases. This is especially true for the critical density ρ c and the highdensity branches of the spinodal and coexistence curves. In order to assess this sensitivity, the high-density boundary ρ 0 was initially set at ρ * 0 = 1, and then moved forward to values as high as ρ * 0 = 1.4, which is very close to the close-packing density ρ * 0 = √ 2. Moreover, for each ρ 0 several different approaches for the boundary condition were tried. Specifically, the internal energy u(ρ 0 , β) has been determined by the HTA; the ORPA; the "non-linear

ORPA" [16], whereby β is replaced by the Mayer function e β * -1 as the amplitude of the potential in c(r), i.e., one sets K * = e β * -1 in Eq. ( 7); the EXP [START_REF] Hansen | Theory of Simple Liquids[END_REF], which is obtained by setting g EXP (r) = g HS (r) exp[∆h ORPA (r)], where g(r) is the radial distribution function, and ∆h ORPA (r) = h ORPA (r) -h HS (r); the linearized version of EXP known as LIN [START_REF] Hansen | Theory of Simple Liquids[END_REF], where one has g LIN (r) = g HS (r)[1 + ∆h ORPA (r)]; and, for very narrow wells, the solution of the PY equation in the SHS limit [2]. For the values of δ discussed below, the results were found to be enough stable with respect to a further increase of ρ 0 or a change of the approximation used for u(ρ 0 , β) to warrant comparison with other predictions.

In Fig. 1 we have reported the SCOZA coexistence and spinodal curves for δ = 0.25 obtained with ρ * 0 = 1.2 and u(ρ 0 , β) given by the non-linear ORPA. The critical temperature and density predicted by the SCOZA are T * c = 0.769, ρ * c = 0.355. These values are close, but not identical to those found in [18] for the same value of δ, namely T * c = 0.761, ρ * c = 0.343. The difference is to be traced back to both the somewhat lower value of ρ 0 used there, ρ * 0 = 1.15, and the fact that the boundary condition at ρ 0 was imposed on ∂ 2 u/∂ρ 2 rather than u. As observed in Sec. II, we found that this enhanced the sensitivity of the results to the choice of the approximation used at ρ 0 , namely, the HTA in the case of [18]. The main difference between the two coexistence curves is that the liquid branch of the present one is slightly shifted to higher densities with respect to that obtained in [18]. Nevertheless, the two curves are close to each other, and the comparison with the simulation results [29][30][31][32][33][34] shown in Fig. 1 confirms the trend already found there, i.e., the liquid branch of the Triangles: simulation data by Liu et al. [29] Inverted triangles: simulation data by del Río at al. [32]. Circles: simulation data by Orea et al. [34]. Crosses: simulation data by Pagan et al. [33].

The SCOZA spinodal has also been shown by the dashed line.

SCOZA coexistence curve underestimates the simulation results. The different simulations reported in the figure are enough consistent with one another to indicate that this is indeed an inaccuracy of the SCOZA. The critical temperature is in good agreement with most simulation results [29,30,[32][33][34], while the critical density is underestimated with respect to all of them, although it lies within the error-bar given in [30], and the discrepancy is generally smaller than that observed on the liquid branch of the coexistence curve.

Figures 2,3 show the SCOZA spinodal and coexistence curves respectively for δ = 0.15 and δ = 0.1, together with simulation data for the coexistence curve and the critical point [10,29,[33][34][35][36]. The SCOZA results were again obtained using the non-linear ORPA as the boundary condition at ρ 0 , and setting ρ * 0 = 1. [START_REF] Noro | [END_REF]. Unlike what found for δ = 0.25, the coexistence curve of the SCOZA is wider than most simulation predictions. This behavior is largely due to the fact that now the SCOZA critical temperature is appreciably higher than the simulation ones. As a general remark, the overall agreement with the phase diagram predicted by the SCOZA appears to be worse than that found for HCY potentials of comparable range [13][14][15], although for δ = 0.15 and δ = 0.1 a quantitative assessment of the performance of the SCOZA is difficult, because of the rather large discrepancies shown by different simulations.

A perhaps more conclusive comparison can be made with the recent simulation results for the critical point obtained in [10], some of which have been reported in Tab. I together with the SCOZA predictions. These indicate that the SCOZA for the SW fluid underestimates ρ c , the discrepancy on ρ c being larger than that on T c . In a previous investigation of the SCOZA for narrow HCY fluids [15], ρ c was again found to deviate from simulations more than T c , although in that case ρ c was overestimated by the SCOZA.

A qualitative feature common to the SCOZA phase diagrams shown in Figs. 1, 2, 3 is that the coexistence curve is much wider then the spinodal curve, so that even at densities near the critical one, and a fortiori at low or high densities, there is a sizable temperature interval where the fluid may survive in a single phase as a metastable state, before spinodal decomposition sets in. According to this, in the experimental practice one should be circumspect before using spinodal decomposition to identify equilibrium phase separation on the assumption that the coexistence and spinodal curves run close to each other. In fact, according to the SCOZA they do not, especially at the low volume fractions that are frequently considered in experiments on colloidal dispersions.

For wells such that δ ≤ 0.05, the dependence of the SCOZA on the high-density boundary condition becomes so strong, that the results do not show any clear trend towards saturation in the region ρ * 0 ≤ 1.4. It is possible that this could be achieved by moving ρ 0 to even higher values, bigger than that corresponding to the physical close packing. However, with the fully numerical algorithm considered here we failed to obtain solutions of Eqs. ( 6), (7) beyond close packing, because in this regime the correlation functions become very singular, and even the prescription described in Appendix A proved insufficient, at least for the number of points in the numerical Fourier transform that can be handled in practice. Therefore, the results that we have obtained for δ ≤ 0.05 cannot be considered meaningful.

It is worthwhile pointing out that the sensitivity on the boundary condition at high density is much stronger for the SW potential considered here than for the HCY potential.

In order to compare the two cases, it is useful to resort to the NF mapping [START_REF] Noro | [END_REF] According to the estimate given in [10] for the case δ = 0.05 (not reported here), the errors in the MC critical parameters are about ±0.14% for T * c and ±1.56% for ρ * c .

Sec. I, according to which a fluid of hard-core particles of diameter σ with a given attractive tail potential at reduced temperature T * can be mapped into a SW fluid with the same σ, and a range determined so that its second virial coefficient B 2 (T ) coincides with that of the original fluid for T = T * . We remark that here this mapping is used without any assumption as to the behavior of the system for interactions of vanishing range, so its use is fully legitimate also within the SCOZA. A HCY fluid with inverse-range parameter z = 5.5

at T * = 0.48, which is the critical value according to the SCOZA, should then be equivalent to a SW fluid with δ = 0.1 at the same T * . The SCOZA critical temperature of this SW,

T * c = 0.49, is indeed close to that of the HCY, in agreement with the prediction from the mapping, while the difference between the critical densities, ρ * c = 0.459 for the SW and ρ * c = 0.427 for the HCY, although larger, is still below 10%. However, for the HCY fluid setting ρ * 0 = 1 is sufficient to obtain a stable critical point, while for the SW ρ * 0 must be set at a substantially higher value, ρ * 0 = 1.4. Further insight into the difference between the two potentials is obtained by checking the behavior of the effective range δ of the HCY below T c : we find that, as the temperature decreases, the effective range also decreases. As a consequence, for T < T c it should be possible to map the densities at coexistence of the HCY fluid into those of a SW fluid at the same T * , but with a lower T * c , i.e., at a higher rescaled temperature T /T c . This implies that, according to the NF mapping, the coexistence curve of the HCY fluid in the T -ρ plane should be narrower than that of the "equivalent" SW at T * c . That this is actually the case is shown in the upper panel of Fig. 4, where the coexistence and spinodal curves of the HCY fluid with z = 5.5 and the SW fluid with δ = 0.1 have been plotted using rescaled temperatures T /T c and densities ρ/ρ c . It is evident that the SW coexistence curve is appreciably wider, especially when the liquid branch is considered. It is also worthwhile considering the behavior of the coexistence and spinodal curves when the temperature is replaced by B 2 (T ), or equivalently by the stickiness parameter τ obtained by replacing B SHS The dotted line is a guide to the eye.

τ c = 0.115 for the SW and τ c = 0.107 for the HCY. If the rescaled temperature T /T c is replaced by the rescaled stickiness τ /τ c , as in the lower panel of Fig. 4, the vapor branches of the coexistence and spinodal curves become nearly coincident, and the discrepancy between the liquid branches of the coexistence curve becomes substantially smaller, although the SW still gives a wider curve. Quite curiously, the situation with the liquid branches of the spinodals is reversed with respect to that shown in the upper panel, i.e., as a function of τ the SW spinodal is narrower than the HCY spinodal, while as a function of T it is wider. On the basis of the above considerations, the sensitivity of the SW potential to the high-density boundary condition may depend at least in part on the fact that its phase diagram is shifted towards higher densities than those of a HCY fluid of equivalent range.

The behavior of B 2 (T c ), shown in Fig. 5, is in semi quantitative agreement with the simulation results [10]. In both cases, the B 2 (T c ) increases as δ decreases. However, because of the impossibility to obtain meaningful results for the critical point of narrow wells such that δ ≤ 0.05, we are not able to say what the fate of the SCOZA B 2 (T c ) will be in the limit δ → 0, namely, if it will tend to a finite value as predicted by the simulations [7,10].

The arguments put forth in Sec. IV rather suggest that the SCOZA would not yield a finite value of B 2 (T c ) for vanishing attraction range.

We now briefly consider some thermodynamic properties of the SW fluid. Figure 6 displays the compressibility factor βP/ρ and the internal energy per particle U * for δ = 0.1 along the isotherms T * = 2 and T * = 1, compared with simulation results [37]. The agreement is very satisfactory both at low and high density, although it is fair to point out that the temperatures at which the simulations were performed can be considered quite high, since even T * = 1 is still more than twice the expected critical value reported in Tab. I. A comparison encompassing temperatures closer to T c is shown in Fig. 7, where βP/ρ and U * for δ = 0.25 are shown along the isochore ρ * = 0.42 as a function of the inverse temperature β, and compared with simulation results [38]. The agreement is again quite good down to the lowest temperature attainable by the SCOZA, i.e., that at which the isochore hits the high-density branch of the spinodal curve. Inspection of Fig. 2 of [38] suggests that the simulation data at the lowest temperatures, roughly for β 1.33, would correspond to unstable states in the thermodynamic limit.

It is also worthwhile examining how the present theory performs for the correlations of the fluid. Figure 8 shows the structure factor S(k) for a SW with δ = 0.1 and two different states, one at moderate density ρ * = 0.5 and relatively high temperature T * = 2, and another at lower density ρ * = 0.229 and temperature T * = 0.5 near the critical value.

In both cases, the SCOZA reproduces quite accurately the simulation results [39]. The correlations in real space are considered in Fig. 9, where the radial distribution function al. [39] extracted from their Figs. 6 and8.

closely with the simulation, both inside and outside the well region. This is not surprising since, as pointed out in Sec. I, the SCOZA for the SW potential becomes exact at low density, while this is not the case for other tail interactions. On the other hand, as ρ is increased the SCOZA overestimates both the contact value g(σ + ) of g(r), and the amplitude of the discontinuity at the well boundary ξ ≡ (1 + δ)σ. This also testifies the quite special status of the SW potential within the SCOZA: in fact, for the HCY potential [13] the SCOZA was found to underestimate g(σ + ), as one would expect from a theory whose c(r)

is linear in the perturbation potential. In the present case, the behavior at high density is instead reminiscent of that of the "non-linear ORPA" studied in [16]. Another feature worth attention is the hump displayed by the simulation g(r) for ρ * = 0.8 and r * ≃ 1.8, which is the precursor of the discontinuity for r * = √ 3 displayed by g(r) in the sticky limit [8,41]. This is not reproduced by the SCOZA, most likely because of the assumption implied by Eq. ( 7) that c(r) has the same range as the potential.

A closer look at the g(r) at contact and at the well boundary is given in Fig. 10, where g(σ + ), g(ξ -), g(ξ + ) for δ = 0.1 are shown along the isotherms T * = 0.5 and T * = 1. As one may expect, the aforementioned discrepancy with the simulation data that sets out at large ρ becomes less severe as T is increased, the results for T * = 1 being in good agreement with the simulation [40] even at the highest ρ considered.

IV. THE STICKY LIMIT OF THE SCOZA

As discussed in Sec. III, the numerical solution of the SCOZA for the SW potential developed here becomes unpredictive for wells such that δ ≤ 0.05, still considerably wider than the HCY potentials that were considered in [14,15]. However, in the limit δ → 0 it is possible to study the behavior of the SCOZA by adopting a different standpoint, namely, by imposing thermodynamic consistency on the solution to Eqs. ( 6), ( 7) after the limit δ → 0 has been taken. The clue to this procedure lies in the observation made in [42] that in the 24 limit δ → 0, any OZ closure like Eq. ( 7), such that one has h(r) = -1 for r < σ and c(r) = 0 for r > (1 + δ)σ, gives a fixed form for c(r):

c(r) = - ηλ r * -a 2 + η 6(a + b) 2 -aλ r * - ηa 2 2 r * 3 Θ(1 -r * ) + λ 12 δ(r * -1) , (26)
where Θ(x) and δ(x) are respectively the Heaviside and Dirac functions of argument x, and we have set η ≡ πρσ 3 /6. The quantities a and b are given by

a = 1 + 2η (1 -η) 2 - λη 1 -η , ( 27
) b = - 3η 2(1 -η) 2 + λη 2(1 -η) . ( 28
)
Therefore, in the limit δ → 0 all the specific features of the OZ closure considered are conveyed into a single, state-dependent parameter λ. The latter is related to the stickiness parameter τ defined in Eq. ( 3) and to the contact value of the cavity function y(r) = g(r) exp [βv(r)] by the expression

λ = y(σ) τ . (29)
By applying the compressibility rule (9) to Eq. ( 26), one obtains

1 χ red = a 2 , (30)
whereas the energy route gives [43] βA

N = βA HS N -η +∞ τ dτ ′ λ(η, τ ′ ) τ ′ . ( 31
)
If the form of λ(η, τ ) obtained by the PY closure is substituted into Eqs. ( 26)-( 31), the aforementioned solution of the SHS model due to Baxter [2] is recovered. On the other hand, as discussed in detail in [42], different choices for λ are also possible. Since one of the most problematic features of the original SHS model is its lack of thermodynamic consistency, it is tempting to determine λ via the SCOZA, thus enforcing consistency between the compressibility and energy routes given by Eqs. ( 30), (31). If we introduce the quantities

Λ = η 2 λ , (32)
t = 1 τ , (33)
the SCOZA equation reads where D(η, t) is given by

D(η, t) t ∂Λ ∂t = ∂ 2 Λ ∂η 2 , (34) 25
D(η, t) = 2 η (1 + 2η) -(1 -η)Λ [ η (1 -η)] 3 . ( 35
)
The PDE (34) can be solved numerically by the finite-difference scheme described in the Appendix of [44], or the similar one given here by Eq. (A1); unlike in the case of non vanishing δ, the diffusion coefficient is given explicitly by Eq. ( 35), so there is no need to evaluate it by finite differences as in Eq. (A2). As usual, for the problem to be completely specified, Eq. ( 34) must be complemented with the initial and boundary conditions. As for the initial condition, we observe that the numerical integration cannot be started exactly from t = 0, because for this value of t the diffusion coefficient 1/(tD) would become infinite.

Instead, the integration was started from a small, but non null t 0 > 0, using as the initial condition for Λ the leading term in an expansion in powers of t 0

Λ(η, t 0 ) ≃ t 0 ϕ(η) . (36)
By substituting Eq. ( 36) into Eq. ( 34) we find that the η-dependent term ϕ(η) must satisfy the ordinary differential equation 1

d 2 ϕ dη 2 = 2 1 + 2η η 2 (1 -η) 3 ϕ . (37)
This is solved numerically with the conditions

ϕ(0) = dϕ dη η=0 = 0 , (38)
d 2 ϕ dη 2 η=0 = 2 , (39)
that follow from Eqs. ( 29), ( 32), (36) and the property [42] lim

η→0 y(σ) = 1 . (40)
We have verified a posteriori that the numerical solution of Eq. ( 34) is virtually independent of the particular choice of t 0 , as long as the latter is sufficiently small. Equations ( 29), (32), also imply for the low-density boundary condition Λ(η = 0, t) = 0 for every t.

1 Although for τ → ∞ the contact value of the radial distribution function coincides with the PY result, y(σ) does not because of the double limiting procedure δ → 0, τ → ∞ in the prefactor e βv(σ + ) . As observed in Sec. II, the choice of the high-density boundary condition entails a certain amount of arbitrariness. A natural possibility is using for λ the expression λ Bax obtained by the analytical solution of the SHS model within the PY approximation [2]:

26

Λ(η 0 , t) = η 2 0 λ Bax (η 0 , t) for every t, (42)
where η 0 denotes the position of the high-density boundary. Quite disappointingly, although perhaps not surprisingly, the outcome of the numerical integration of Eq. ( 34) with the initial and boundary conditions given by Eqs. ( 36), ( 41), ( 42) suffers from the same problem experienced in Sec. III for wells of non vanishing width, namely, a strong dependence of the results on the position of η 0 . This is illustrated in Fig. 11, where the inverse reduced compressibility 1/χ red obtained for different choices of η 0 has been plotted along the isotherm τ = 0.11, which is close to the critical one according to numerical simulations [7]. For the solution to be deemed reliable, the various curves, at least for η 0 greater than a certain threshold, ought to lie near one another -a requirement that is manifestly not satisfied.

Some of the isotherms are not even consistent with τ = 0.11 being close to the expected critical value. We have also verified that this state of affairs is not changed by different choices of the high-density boundary condition.

We can, however, resort to a different method to solve Eq. (34), that does not require the high-density boundary to be specified. This is made possible by the fact that, unlike in the case of wells of non vanishing width discussed in Secs. II, III, the diffusion coefficient of Eq. ( 34) is known explicitly as a function of η and τ . Let us assume that Λ can be expressed as a power series in t

Λ(η, t) = ∞ n=1 t n ω n (η) , (43)
where the zeroth-order term has been omitted, since Λ must vanish for t → 0. By inserting Eq. ( 43) into Eq. ( 34) and equating the terms of the same degree in t, Eq. ( 34) is replaced by an infinite set of ordinary differential equations:

d 2 ω n dη 2 = 2n 1 + 2η η 2 (1 -η) 3 ω n - 2 η 3 (1 -η) 2 n-1 j=1 j ω j ω n-j , (44)
which gives back Eq. ( 37) for n = 1. Since this set contains only backward references, it can be integrated numerically to any desired order n. Only the low-density behavior of ω n is required. Specifically, one needs the values of ω n , dω n /dη, and d 2 ω n /dη 2 at η = 0 for n ≤ n. 34) with a high-density boundary condition computed in the PY approximation [2]. The different lines refer to the boundary packing fractions η 0 = πρ * 0 /6

given in the legend. The stickiness parameter is set at τ = 0.11, which is close to the critical value according to numerical simulations [7].

If it is assumed that ω n can be expanded in power series of η, it is found that for n > 6 all the above quantities are in fact vanishing, while for n ≤ 6 they can be obtained from exact diagrammatic expansions available in the literature [42,45]. Unfortunately, the radius of convergence t R of the series (43) corresponding to most values of the packing fraction η comes out to be too small for it to be of direct usefulness. This is apparent from Fig. 12, which shows how t R rapidly decreases as the density is increased. In particular, at the critical value of t predicted by simulations [7], ρ * should be well below 0.2 to obtain convergence. This is again disappointing, yet not uncommon: a power series obtained as a formal solution of a differential equation may easily turn out to be divergent. In those cases, one can try a resummation technique in order to provide an analytic meaning to such a divergent series.

To this end, we have employed Borel summation method [START_REF] Kawai | Algebraic analysis of singular perturbation theory[END_REF], whereby one starts from the exact equality and formally interchanges the series and the integral operators to obtain

Λ(η, t) = ∞ n=1 n! t n n! ω n (η) = 1 t ∞ n=1 +∞ 0 ds e -s/t s n n! ω n (η) (45) 28
Λ B (η, t) = 1 t +∞ 0 ds e -s/t Λ(η, s) , (46)
where Λ(η, s) is the regularized power series

Λ(η, s) = ∞ n=1 s n n! ω n (η) , (47)
that has an improved radius of convergence because of the factorial in the denominator.

Where the original series is convergent, the two functions Λ and Λ B are in fact identical.

The latter, however, yields a finite value in a much larger domain, providing the analytic continuation of the function defined by the original series. If one is willing to trust the Borel sum Λ B in this enlarged region, then a puzzling result is found: the liquid-vapor phase transition disappears, at least for reasonable values of the packing fraction and temperature.

In by means of Borel summation method gives instead a finite χ red even if the temperature is lowered further. This result is qualitatively similar to that found by direct numerical integration of Eq. ( 34) when the high-density boundary η 0 is pushed to (unphysically) high values, but is clearly at odds not only with the PY solution of the SHS model [2], but also with the simulation results [7]. The present investigation supports the conclusion that in the limit of vanishing well width δ → 0, the SCOZA critical temperature decreases too rapidly to give a finite value of B 2 (T c ), so that B 2 (T ) or, equivalently, the stickiness τ , do not represent a convenient way to measure the temperature in this limit. Imposing thermodynamic consistency on the SHS model does not lead to a better quantitative agreement with simulation compared to the highly inconsistent PY solution, but merely makes the liquid-vapor transition disappear. We have performed an investigation of the SCOZA for SW potentials, focused on the regime of narrow wells which is the most relevant for the modelization of colloidal particles.

This study was prompted by the lack of a theory able to account quantitatively for the liquid-vapor phase transition in a system of sticky hard spheres: in the well known and widespreadily used PY approximation [2], the various routes to the thermodynamics of the fluid are inconsistent with one another, and agree only qualitatively with the results of computer simulation [7]. On the other hand, for HCY potentials the SCOZA can be solved semi-analytically, and has already been usefully applied to narrow-range tails [14,15]. The rationale for turning to the SW potential within the SCOZA is that this interaction has the unique feature of having the SCOZA become exact at low density, thereby recovering the correct B 2 (T ) which plays such an important role for narrow interactions.

However, the results presented here show that this peculiarity turns more into a liability than an asset. For wells of non vanishing width δ, after extending the study already made in [18] to (slightly) smaller δ, we cannot but confirm the conclusions of that work, namely: i) the results are found to be very sensitive to the boundary condition at high density that is needed to integrate the SCOZA PDE, actually much more so than in the case of HCY potentials of comparable range. Such an undesired feature can be tamed by pushing the density of the boundary to very high values, but this has the effect of making the numerical solution of the PDE unfeasible for δ 0.05. ii) In the range of δ for which meaningful results can be extracted, the overall agreement with simulation data, at least for the phase diagram, is actually worse than what was found for the HCY potential [13][14][15], although a sharp assessment of the performance of the SCOZA is made difficult by the considerable discrepancies between different simulations.

If the sticky limit is taken from the outset, and subsequently thermodynamic consistency Stating that the SCOZA is inadequate for SW interactions in the sticky limit is merely an observation, and does not shed any light on the reason behind the fact. In this respect, one should recall the result of Stell [9] mentioned in Sec. I, who showed that the SHS model is strictly speaking ill-defined if a class of problematic configurations of particles is not ruled out, for instance by introducing a small amount of polidispersity into the system. Can it be that the (possibly more accurate) SCOZA approach exposes the malicious effect of the diverging clusters, which are simply neglected in the PY approximation? We regard this as unlikely, since the two approximations are quite similar in their handling of the correlation functions of the fluid, giving the same functional form for c(r) in the sticky limit [42]. What we deem most probable, instead, is that the failure is a combined effect of the approximate closure (7) and the consistency requirement (8). Indeed, the closure is likely to miss a number of singular features, i.e., discontinuities and Dirac-δ peaks, that show up in the real correlation function [8,41]: the isothermal compressibility, being dependent on the volume integral of the same function, suffers from the same defect. On the contrary, the internal energy of the fluid depends on the contact value of the cavity function alone:

when thermodynamic consistency is imposed, the latter has to compensate for the missing contributions, unnaturally pushing away the liquid-vapor phase boundary. We stress that at the moment this is a mere, untested hypothesis.

One may hope that the strong dependence of the results on the high-density boundary condition found for narrow SW potentials will become less drastic turning to a different form of the interaction. In fact, as observed above, for the HCY potential this unwelcome feature, albeit still there, is nevertheless much less pronounced. If this is actually the case, then the fully numerical solution of the SCOZA which has been developed here and the one previously developed in [18], could still be usefully applied to different forms of short-range potentials, not to mention the instance of interactions that are not of short range [START_REF] Betancourt-Cardenas | [END_REF]. Moreover, this solution could be helpful in the context of the renormalization-group based hierarchical reference theory (HRT) [48], specifically for enforcing the core condition on g(r) exactly (of course within the limit of numerical accuracy). In fact, in the smooth cut-off formulation of the HRT [49], similarly to the SCOZA, this condition can be (and has been) implemented analytically for a HCY potential, but for a generic tail potential a numerical solution would again be needed. Also, this is always the case for the sharp-cut off HRT, irrespective of the form of the tail potential. In the latter theory, the core condition was implemented by introducing a number of approximations [50], which however affect the accuracy of this approach for short-range interactions [13]. A better procedure such as that developed here would certainly be more adequate.

As for the sticky limit, the results presented in this work suggest that, if thermodynamic consistency is to be employed as a way to improve upon the accuracy of the PY solution for the SHS potential (2), one should first get rid of the approximation common to the SCOZA and the PY equation: namely, that the contribution to the direct correlation function c(r) due to the tail interaction w(r) has the same range as w(r) itself, and therefore vanishes wherever w(r) does. In this respect, it is worthwhile recalling that the same approximation is used in conjunction with thermodynamic consistency not only in the SCOZA, but also in the HRT [49]. Therefore, if the combination of these ingredients is bound to become dangerous for vanishing attraction range, it is likely that this would affect the HRT as well, and that both the SCOZA and the HRT would benefit from a solution to this problem.

(ρ j , u n+1/2 j ) = 1 χ red n+1 j - 1 χ red n j u n+1 j -u n j , (A2)
where subscripts refer to the density ρ and superscripts to the inverse temperature β in the (ρ, β) grid, so for instance u n j ≡ u(ρ j , β n ), and ∆β, ∆ρ are respectively the inverse temperature and density steps used in the discretization. Equation (A1) has been solved by an implicit algorithm as follows: assuming that u n j is known at a certain inverse temperature β n for every ρ j , a guess for u n+1 j at the "new" temperature β n+1 is given. This is determined by solving Eqs. ( 6), (7) with the approximation with the well depth. Note that, for λ to be the same dimensionless quantity in Eqs. (B1) and (B2), the differentiation of S has to be performed with respect to K * and √ ρφ(r).

K n+1 ≃ K n +
To implement the algorithm, one starts from a guess for K and φ(r) for r < σ, and determines φ(r) for r > σ by the relation φ(r) = -Kw(r). In the numerical solution of the SCOZA, the guess at the inverse temperature β n+1 and density ρ j was given by φ(r, ρ j , β n+1 ) ≃ φ(r, ρ j , β n ), K(ρ j , β n+1 ) ≃ K(ρ j , β n ) + ∆β. By switching to reciprocal space, the Fourier transform of ∆h(r) is determined via the OZ equation as

∆ ĥ(k) = S HS (k)S(k) φ(k) , (B4)
where S(k) is given by Eq. (19). Equation (B4) is then inverse-transformed to give ∆h(r).

This in turn is used to determine ∆V by Eqs. ( 10), ( 21), (B3), and the ∆h(r) and ∆V thus obtained are used in Eqs. (B1), (B2) to update φ(r) and K iteratively, until convergence is achieved.

As observed in [26], a valuable improvement over the steepest descent is the conjugate gradient algorithm, which we have adopted here. In this method, the direction of descent at the step l + 1 does not coincide with that of the gradient of S at step l, but is instead given by a linear combination of the gradients at all previous steps according to

φ l+1 (r) = φ l (r) -λ γ l (r) r < σ , (B5)
K * l+1 = K * l -λ Γ l , (B6)
where γ l (r) and Γ l are defined by recurrence as [26] γ l (r) = ∆h l (r) + α l γ l-1 (r) r < σ , (B7)

Γ l = ∆V l + α l Γ l-1 , (B8)
with α 0 = 0 and α l for l ≥ 1 given by

α l = ρ r<1 d 3 r ∆h l (r) [∆h l (r) -∆h l-1 (r)] + ∆V l (∆V l -∆V l-1 ) ρ r<1 d 3 r [∆h l-1 (r)] 2 + (∆V l-1 ) 2 . ( B9
)
The numerator and denominator of Eq. (B9) are scalar products of the form v l |v l -v l-1 , v l-1 |v l-1 for the present case in which the "vector" v l consists of both the function √ ρ∆h l (r) and the scalar ∆V . In the implementation of the algorithm, special attention has been paid to the choice of the parameter λ. In both the steepest descent and conjugate gradient methods, the optimal value of λ at each step is determined by line minimization [26]. In the case of the conjugate gradient, this amounts to minimizing S[φ l (r) -λγ l (r), K * l -λΓ l ] as a function of λ for fixed φ l (r), γ l (r), K * l , Γ l . A possible way would be to use some method based on repeated evaluation of the functional S along the line φ l (r) -λγ l (r), K * l -λΓ l . We found that this procedure is not recommended here for a number of reasons: first, it is not very efficient, since it requires many evaluation of S at each step. Second, at low temperature (see Sec. II) there is a density interval where the argument of the logarithm in Eq. ( 16) is not positive definite as a function of k, so great care has to be taken to prevent the functional S from being evaluated for these states [26]. Finally, there is a subtle, yet important reason related to the fact that, for the hard-core plus tail potential of Eq. ( 5), both c(r) and h(r) are discontinuous at r = σ, and display additional discontinuities if w(r) is itself discontinuous, as in the case of the SW potential. In order to determine with good accuracy the Fourier transform of the correlation functions needed in the iterative procedure, it is then advisable to subtract off the discontinuous contribution before performing the numerical Fourier transform, and add to the result its Fourier transform determined analytically, as illustrated in Sec. II. However, when the functional f Ring of Eq. ( 16) is discretized in the numerical calculation, the fundamental relations (17), (18) will hold rigorously if the numerical Fourier transforms of φ(r) and ∆h(r) are performed on the full functions, rather than just on their continuous part.

As a consequence, the gain in accuracy achieved by splitting the functions to be transformed entails as a drawback that the solution of the Euler-Lagrange equation δf Ring /δφ(r) = 0 for r < σ will not rigorously be anymore a minimum of f Ring as defined by Eq. ( 16), and it is not obvious to us whether Eq. ( 16) can be modified so as to recover this property. This is one more reason why we chose to perform the minimization with respect to λ by an algorithm which would not rely on the direct evaluation of S, and hence of f Ring . Instead, we set ψ(λ) = S[φ l (r) -λγ l (r), K * l -λΓ l ], and determined the line minimum by solving the equation dψ/dλ = 0 with respect to λ. This was done by the Raphson-Newton method, where we have understood that γ(r) is set identically to zero for r > σ, and for brevity we have omitted the index l in φ(r), γ(r), K * , Γ. All the derivatives of S are evaluated at φ l (r) -λ k γ l (r), K * l -λ k Γ l , and are easily calculated by Eqs. ( 17), ( 18), ( 20), ( 22), (23). The quantities required are those already needed in Eq. (B4) for ∆ ĥ(k). Specifically, Eqs. In practice, we have found that, for each conjugate gradient step (B5), (B6), (B7), (B8), a single Raphson-Newton step (B10), (B13), (B14) for λ is sufficient. Iterating Eq. (B10) so as to get closer to the minimum of ψ(λ) did not lead to any significant improvement in the overall speed of convergence of the algorithm.

We notice that at the critical point and on the spinodal curve, where the compressibility χ red diverges, the integrals that appear in Eq. (B14) diverge as well, because S 2 (k) develops a non integrable singularity at k = 0, see Eq. (19). In order to account for this behavior, for large values of χ red the contribution to the integrals in the small-k domain is determined analytically by replacing S(k) with its Ornstein-Zernike form S(k) ≃ 1/(χ -1 red + bk 2 ), and evaluating the remainder of the integrand at k = 0. where ĉ(k) has been subtracted off from ĥ(k) in order to make the integrand decrease more rapidly at large k. The function ζ(k) = k ĉ2 HS (k)/[1 -ρĉ HS (k)], which is known analytically for the Waisman parametrization, is then sampled with a step ∆k ≃ 5 × 10 -5 σ -1 , much smaller than that used in the numerical Fourier transform, so that the "pseudo-Bragg" peaks can be included in the sampling. Each peak such that the structure factor exceeds the threshold value S HS (k) = 3 is modeled by a Lorentzian A/[1 + α 2 (k -k 0 ) 2 ], where the parameters k 0 , A, α are determined so as to reproduce the position of the peak, its height, and its curvature. The sum of the Lorentzians is subtracted off ζ(k), resulting in a smooth function that can be transformed numerically without substantial loss of information, while the Fourier transform of the Lorentzians is determined analytically. In the present case in which k 0 = 0 this cannot be done exactly, but we found that for narrow peaks such that α ≫ 1 like those considered here, a satisfactory approximation is given by where γ = 0.5772 . . . is Euler constant.

[

  consistency as a way to improve the PY solution for the SHS potential (2) would require to give up the ansatz on the range of c(r).

11

 11

13

 13

  FIG. 1: Liquid-vapor coexistence curve and critical point of the SW fluid of width δ = 0.25 in the density-temperature plane. Solid line: SCOZA. Full circle: SCOZA critical point. Squares with error bars: simulation data by Vega et al.[30]. Diamonds: simulation data by Elliott et al.[31].

14 PageFIG. 2 :

 142 FIG.2: Same as Fig.1for δ = 0.15. Solid line: SCOZA. Full circle: SCOZA critical point. Circles: simulation data by Orea et al.[34]. Crosses: simulation data by Pagan et al.[33]. Triangles: simulation data by Liu et al.[29]. Dashed line: SCOZA spinodal.

15 Page 15 of 41 URLFIG. 3 :

 15413 FIG. 3: Same as Fig. 1 for δ = 0.1. Solid line: SCOZA. Full circle: SCOZA critical point. Squares with error bars: simulation data by Duda [36]. Circle with error bars: simulation critical point by Skibinsky et al. [35]. Diamond: simulation critical point by Largo et al. [10]. Dashed line: SCOZA spinodal.

mentioned in 16 Page

 16 Critical temperature and density of the SW fluid for several values of the well width δ, as obtained by the SCOZA and Monte Carlo (MC) numerical simulations performed in[10].

17

 17

FIG. 4 :

 4 FIG. 4: Solid line: SCOZA coexistence (outer) and spinodal (inner) curves for a SW fluid of width δ = 0.1. Dashed line: same for a HCY fluid with inverse-range parameter z = 5.5. The upper panel shows the curves in the temperature-density plane, with both variables rescaled by their critical values. In the lower panel the rescaled temperature has been replaced by the rescaled stickiness parameter τ /τ c [2]. Note how the curves for the two potentials get closer in the latter case.

2 with B 2 FIG. 5 :

 225 FIG. 5: Reduced second virial coefficient of the SW fluid at the critical temperature B * 2 (T c ) as a function of the well width δ. Squares: SCOZA. Circles: simulation results by Largo et al. [10].

FIG. 6 :

 6 FIG. 6: Compressibility factor βP/ρ (upper panel) and excess internal energy per particle U * (lower panel) of the SW fluid with δ = 0.1. Solid line: SCOZA for T * = 1. Dashed line: SCOZA for T * = 2. Circles: simulation results at the same temperatures [37].

20 Page 20 of 41 URLFIG. 7 :

 20417 FIG. 7: Compressibility factor βP/ρ (upper panel) and excess internal energy per particle U * (lower panel) of the SW fluid with δ = 0.25 and ρ * = 0.42. Solid line: SCOZA. Circles: simulation results [38].

g 21 Page 21 of 41 URLFIG. 8 :

 21418 FIG. 8: Structure factor of the SW fluid with δ = 0.1 for T * = 2, ρ * = 0.5 (upper panel) and T * = 0.5, ρ * = 0.229 (lower panel). Solid line: SCOZA. Circles: simulation data by Shukla et

22

 22

FIG. 9 :

 9 FIG. 9: Radial distribution function of the SW fluid with δ = 0.1 for T * = 0.5 and ρ * = 0.2 (top panel), ρ * = 0.4 (central panel), ρ * = 0.8 (bottom panel). Solid line: SCOZA. Circles: simulation data by Largo et al. [40] extracted from their Fig. 6.

23 Page 23 of 41 URLFIG. 10 :

 234110 FIG.10: Radial distribution function of the SW fluid with δ = 0.1 for r → σ + , r → ξ -, r → ξ + , where ξ ≡ (1 + δ)σ is the value of r beyond which the potential vanishes. Squares, triangles, and circles represent the SCOZA results for g(σ + ), g(ξ -), and g(ξ + ) respectively. Crosses, pluses, and starbursts represent the simulation data obtained by Largo et al.[40] for the same quantities. Dotted lines are a guide to the eye. Upper panel: T * = 1 isotherm. Lower panel: T * = 0.5 isotherm.

27

 27

8 FIG. 11 :

 811 FIG.11: Inverse reduced compressibility 1/χ red of the SHS fluid as a function of the reduced density ρ * as found via the SCOZA equation (34) with a high-density boundary condition computed in

FIG. 12 :

 12 FIG. 12: Radius of convergence t R of the power series(43), estimated by means of the ratio convergence test, as a function of the reduced density ρ * . The horizontal dashed line marks the value of the inverse stickiness t corresponding to the critical point according to numerical simulations[7].

Fig. 13 the 29 Page 29 of 41 URLFIG. 13 :

 13294113 FIG.13: Inverse reduced compressibility 1/χ red of the SHS fluid as a function of the inverse stickiness t along the η = 0.25 isochore. The compressibility route of the PY approximation (dashed line) predicts the compressibility to diverge at t ≃ 12. According to the solution of Eq. (34) evaluated by means of the Borel summation method (solid line), the isochore does not cross the spinodal at all.

30 Page 30 of 41 URL

 3041

  is imposed, the SCOZA PDE is obtained explicitly in closed form, but the high-density boundary problem is again found to be looming. If the PDE is turned into an infinite set of ordinary differential equations and the low-density behavior alone is imposed, the liquid-vapor phase boundary disappears altogether. This supports the conclusion that as the attractive potential range vanishes, the SCOZA critical temperature will decrease more rapidly than what is necessary to yield a finite B 2 (T c ).

  ∆β. The values of the inverse reduced compressibility (1/χ red ) n+1 j corresponding to u n+1 j are then obtained by solving Eqs. (6), (7) at fixed density ρ and internal energy u via the method described in Sec. II and Appendix B. These values are used in Eq. (A2) to obtain a trial value of D(ρ j , u n+1/2 j ), which is inserted into Eq. (A1) to give a closed linear system of equations for u n+1 j at β n+1 . The values of u n+1 j thus obtained replace the initial guess, and the procedure is iterated until convergence is obtained. The u n+1 j at convergence are fed back into the algorithm to obtain u n+2 j , and so on.

35

 35

dλ 2 = d 3 r d 3 r ′ δ 2 S

 22 a guess value λ 0 and updating λ recursively by the sequenceλ k+1 =δφ(r)δφ(r ′ ) γ(r)γ(r ′ ) + 2Γ d 3 r ∂ ∂K * δS δφ(r) γ(r) + ∂ 2 S ∂K * 2 Γ 2 , (B12)

= ρ d 3 r

 3 ∆h(r)γ(r) + Γ∆V , 3 [S(k)γ(k)] 2 -2ρΓ d 3 k (2π) 3 S 2 (k) ŵ * (k)γ(k)+ρΓ 2 d 3 k (2π) 3 [S(k) ŵ * (k)] 2 . (B14)

37

 37

  NUMERICAL CALCULATION OF THE HARD-SPHERE COR-RELATION FUNCTION AT HIGH DENSITYWe start from the exact relationh(r) = c(r) + d 3 k (2π) 3 e ik•r [ ĥ(k) -ĉ(k)] = c(r)

2 -

 2 2 cos(k 0 r)Ci(k 0 r) -sin(k 0 r) π Si(k 0 r) ,where Si(x) and Ci(x) are respectively the sine and cosine integrals defined as:

  1] B. Götzelmann, R. Evans, and S. Dietrich, Phys. Rev. E 57, 6785 (1998); M. Dijkstra,

	F o r F o r F o r
	P P P
	e e e
	e r e r e r
	R R R
	e e e
	v i e v i e v i e
	w w w
	O n l O n l O n l
	y y y
	38

R.

van

Roij, and R. Evans, Phys. Rev. Lett. 81, 2268 (1998); M. Dijkstra, R. van Roij, and R. Evans, Phys. Rev. Lett. 82, 117 (1999); M. Dijkstra, R. van Roij, and R. Evans,

Acknowledgments

We would like to dedicate this paper to Bob Evans on his 65 th birthday. We thank Pedro Orea for sending us simulation data for the phase diagram of the SW fluid and the Italian Ministry of Education and Research (MIUR) for PRIN 2008 funding. D. P. would also like to thank Paolo Giaquinta for stimulating the investigation of the behavior of the SCOZA for narrow attractive interactions and for his unflagging interest in this problem.

APPENDIX A: FINITE-DIFFERENCE DISCRETIZATION OF THE SCOZA PDE Equation ( 11) is a nonlinear, diffusive-like PDE, where the quantity D(ρ, u) given by Eq. ( 12) plays the role of the reciprocal of the diffusion coefficient. This PDE is discretized in the following way:

APPENDIX B: NUMERICAL OPTIMIZATION OF THE FREE ENERGY

FUNCTIONAL

The most straightforward method to minimize the functional S defined by Eq. ( 23) is the steepest descent [26], in which φ(r) and K are updated iteratively by moving downhill in the direction opposite to that of the gradient of S, namely

Here the indexes l, l + 1 refer to the values of φ(r) and K at two consecutive iterations, λ

(not to be confused with the state-dependent function λ(η, τ ) introduced in Sec. IV) is a parameter determining the step size, and we have set

The asterisks in Eqs. (B2), (B3) denote reduced quantities K * = Kǫ, ∆U * = ∆U/ǫ, w * (r) = w(r)/ǫ, where ǫ is the energy scale of w(r). For the SW potential considered here, ǫ coincides 34