Rational-function approximation for fluids interacting via piece-wise constant potentials

The structural properties of fluids whose molecules interact via potentials with a hard-core plus n piece-wise constant sections of different widths and heights are derived using a (semi-analytical) rational-function approximation method. The results are illustrated for the cases of a square-shoulder plus square-well potential and a shifted square-well potential and compared both with simulation data and with those that follow from the (numerical) solutions of the Percus-Yevick integral equation.

While the phase diagram and the thermodynamic and transport properties of the discrete-potential fluids have been thoroughly examined, studies of their structural properties, either theoretical or from simulations, are more limited. In fact, for the case where the potential is built as a combination of square wells and shoulders, we are only aware of the rather recent work reported in references [32,38,39].
In previous papers [22,40], following a methodology that, although approximate, has proven successful for many other systems [41], the structural properties of the square-well and the square-shoulder fluids were derived. The main aim of this paper is to use a similar methodology, referred to as the method of rational-function approximation (RFA), to generalize the previous results and derive the structural properties of fluids whose molecules interact via potentials with a hard-core plus piece-wise constant sections of different widths and heights.
It should be emphasized that the RFA is not a new integral equation but rather an alternative approach to derive the structural properties of fluids in an analytic or semi-analytic way. Instead of proposing a closure relation for the Ornstein-Zernike equation, it deals directly with the radial distribution function in Laplace space and involves some coefficients which are determined by imposing (basic) specific physical conditions (see references [22,40,41] for details). An interesting feature of the method is that, when applied to hard-sphere systems in odd dimensions [42][43][44][45], its simplest implementation coincides with the exact solution of the Percus-Yevick (PY) integral equation [46][47][48][49][50][51][52][53]. In the case of even dimensions, the mathematical problem becomes much more complicated, which is connected with the absence of a simple relationship between the Laplace transform (which is a key element of the RFA approach) and the Fourier transform (related to the static structure factor) in even dimensions. In fact, the PY equation for hard-sphere fluids has no analytical solution in those dimensions.
The paper is organized as follows. In section 2 we provide the background material required for the subsequent development. This is followed in section 3 by two different proposals for the RFA, both for the general n-step case and for the particular case of a two-step potential. Section 4 contains the results of our calculations for illustrative cases and their comparison with both simulation data [39] and with those that stem out of the numerical solution of the PY equation for this system. We close the paper in section 5 with further discussion and some concluding remarks. The consistency of the present approach with the results for both the square-well and square-shoulder results is proven in an appendix.

Fundamental relations
The radial distribution function g (r ) of a fluid of particles interacting via a potential ϕ(r ) is directly related to the probability of finding two particles separated by a distance r [54]. It can be measured from neutron or x-ray diffraction experiments through the static structure factor S(q), which is related to the Fourier transform of g (r ) − 1 by where ρ is the number density and is the Laplace transform of r g (r ). The isothermal compressibility of the fluid, κ T = ρ −1 ∂ρ/∂p T , where p is the pressure and T is the temperature, is directly related to the long-wavelength limit of the structure factor: where k B is the Boltzmann constant. Thus, all the physically relevant information about the equilibrium state of the fluid is contained in g (r ) or, equivalently, in G(s).
Let us consider now the piece-wise constant potential This potential is characterized by a hard core of diameter σ and n steps of (positive or negative) heights ǫ j and widths (λ j − λ j −1 )σ (where the convention λ 0 = 1 is understood). Thus, λ n σ denotes the total range of ϕ(r ). The sign of ǫ j defines whether the j -th step is either a "shoulder" (ǫ j > 0) or a "well" (ǫ j < 0). The interaction potential at r = λ j σ ( j = 1, . . . , n) is repulsive if ǫ j > ǫ j +1 and attractive if ǫ j < ǫ j +1 (where, by convention, ǫ n+1 = 0). As usual, the density is measured by the packing fraction η ≡ π 6 ρσ 3 . Henceforth, we will take the hard-core diameter σ = 1 as the length unit.
It is convenient to define an auxiliary function F (s) directly related to the Laplace transform G(s)

23602-2
Piece-wise constant potentials Laplace inversion of equation (5) provides a useful representation of the radial distribution function: where f m (r ) is the inverse Laplace transform of s[F (s)] m and Θ(r ) is Heaviside's step function. Thus, the knowledge of g (r ) [or G(s) or S(q)] is fully equivalent to the knowledge of the auxiliary function F (s). The representation (6) reflects the fact that, due to the hard core at r = 1, the radial distribution function possesses singularities (of decreasing order) at r = 1, 2, 3, . . .. In particular, the value of g (r ) at contact, g (1 + ), is given by the asymptotic behaviour of F (s) for large s: Since g (1 + ) should be finite and different from zero, we get the condition From equation (1), it turns out that the behaviour of G(s) for small s determines the value of S(0): Insertion of equation (9) into the first equality of equation (5) yields the first five terms in the expansion of F (s) in powers of s [40,55,56], and expresses S(0) in terms of the coefficients of s 5 and s 6 , namely Equations (8) and (10) provide the behaviours for large s and small s, respectively, that the auxiliary function F (s) should necessarily satisfy. Let us now decompose F (s) as where, as said before, λ 0 = 1. The aim of this decomposition is to reflect the discontinuities of g (r ) at the points r = λ j where the potential is discontinuous. Let us denote by ξ j (r ) the inverse Laplace transform of sR j (s). Thus, Equation (12) implies that the inverse Laplace transform of sF (s) is From now on we will assume, for the sake of concreteness, that λ n 2 (although the case λ n > 2 can also be dealt with in a similar way). Therefore, equation (6) yields

23602-3
As a consequence, On the other hand, the cavity function y(r ) ≡ g (r )e ϕ(r )/k B T and its first derivative y ′ (r ) are continuous at r = λ j [57]. This means that or, according to equation (15), In the low-density limit, one has g (r ) → e −βϕ(r ) , so that where, apart from ǫ n+1 = 0 and λ 0 = 1, it is understood that ǫ 0 = ∞. Comparison with equation (12) yields

The general n-step case
According to the scheme presented in the preceding section, the full knowledge of the radial distribution function g (r ) for the n-step potential (4) amounts to a prescription for the functions R j (s), such that the function F (s) obtained through equation (12) fulfills the conditions (8) and (10). Now we assume the following rational-function approximation (RFA) for R j (s): Since the degree difference between the numerator and denominator of R 0 (s) is equal to 2, the form (22) ensures the consistency with equation (8). In fact, according to equation (13), Insertion into equation (16) yields

RFA1
As a first proposal to determine the 2n coefficients A j and B j with j = 1, . . . , n, we may discard the density dependence of A j ( j = 1, . . . , n). As a consequence, those coefficients are fixed at their zero-density values, namely Insertion of equation (38) into equation (27) implies that A 0 is also given by its zero-density expression, equation (34). The coefficients B j ( j = 1, . . . , n) are determined from the n constraints (18) stemming from the continuity of y(r ). Making use of equations (23) and (26) yields These conditions should be enforced numerically. We will refer to this variant as the RFA1 method. In the one-step case (n = 1), it coincides with the one previously proposed for square-well [11,40,41,[58][59][60] and square-shoulder [22] fluids. The RFA1 reduces to the PY solution for hard spheres and sticky hard spheres in the appropriate limits [22,40].

RFA2
A more sophisticated version (here referred to as RFA2) consists in replacing the simple prescription (38) by the enforcement of the continuity of y ′ (r ) via equation (19). Again, taking into account equations (23) and (26), one gets Now the problem requires the numerical solution of the 2n coupled transcendental equations (39) and (40), instead of the n equations (39) required in the RFA1.
While the RFA2 is internally more consistent than the RFA1, it is not necessarily more accurate. As we will illustrate in section 4, the requirement of a rather subtle continuity condition of y ′ (r ) at r = λ j may force a radial distribution function with features less realistic than those found.
In either case, once the Laplace transform G(s) is fully determined, its numerical Laplace inversion [61] yields g (r ).

Particularization to two-step potentials
Let us now adapt the previous scheme to the case n = 2. Then, where the three functions R 0 (s), R 1 (s) and R 2 (s) have the rational-function form (22). There are 9 coefficients to be determined. In the simpler RFA1 approach, A j are fixed by their zero-density values, i.e., according to equation (38), Next, S 1 , S 2 and S 3 are given as linear combinations of B j by equations (28)-(30), where Moreover, equation ( These two transcendental equations close the problem. In a more sophisticated RFA2 version of the approximation, equation (42) is replaced by and equation (40), i.e., These two equations should be solved in conjunction with equations (44) and (45).
It is shown in the appendix that the RFA for the two-step potential reduces to the one for the one-step potential in the limit ǫ 1 = ǫ 2 as well as in the limit ǫ 1 → ∞.

Both potentials have recently been analyzed by
Bárcenas et al. [39] by means of exchange replica Monte Carlo (MC) simulations, with special emphasis on the coexistence curves between the vapour and the condensed phase (liquid or solid) and the structural properties of the condensed phase at coexistence.
The chosen state points are ρσ 3 = 0.421, k B T /|ǫ 2 | = 0.64 for the SS+SW system and ρσ 3 = 0.427, k B T /|ǫ 2 | = 0.76 for the sSW system, both lying on the condensed branch of the respective coexistence curve [39]. The radial distribution function g (r ) at a subcritical temperature and at the density of the coexisting liquid represents a rather stringent test of the theoretical approach developed in this paper. To put the results in a proper context, we have also numerically solved the PY integral equation by an iterative method [22]. The comparison between RFA and PY is especially relevant since both theories coincide in the hard-sphere and sticky-hard-sphere limits [22,40]. Figures 1 and 2 show the radial distribution function and the cavity function for the SS+SW and sSW systems, respectively, as obtained from MC simulations [39], the PY numerical solution and our RFA1 and RFA2 approaches. We observe from the figures that the RFA1 and RFA2 curves are practically indistinguishable in the region r > λ 2 , where they describe the main oscillating trends of the true distribution, at least in a semi-quantitative way. Inside the well (λ 1 r λ 2 ), however, both versions of the RFA differ. In the SS+SW case (figure 1), the RFA2 prediction is excellent, visibly differing from the MC data only near the inner radius of the well, r = λ 1 , while the RFA1 exhibits an accused convex shape; as a consequence, in the shoulder region (1 r λ 1 ) the RFA1 predicts a g (r ) with positive slope. In the sSW case (figure 2), on the other hand, it is the RFA2 that presents a convex shape inside the well (λ 1 r λ 2 ) and hence exhibits a positive slope in the shift region (1 r λ 1 ). By contrast, the RFA1, which (as clearly seen in the bottom panel of figure 2) has a discontinuous first derivative at r = λ 2 , presents an almost linear form in the regions 1 r λ 1 and λ 1 r λ 2 with a qualitatively correct behaviour. Note, however, that the RFA2 is very accurate in the region inside the well and near the outer radius (1.4 r λ 2 ).
As for the PY theory, it is obvious from figures 1 and 2 that it performs worse than any of the RFA results for r λ 1 and than either RFA2 (SS+SW system) or RFA1 (sSW system) for 1 r λ 1 . This remarkable result (given that the PY solution requires full numerical work, in contrast to the semi-analytical character of the RFA) was already observed for square-shoulder fluids [22].

Conclusions
In summary, in this paper we have proposed a method of deriving the structural properties of a particular kind of discrete-potential fluids, namely the ones in which molecules interact via a potential consisting of a hard-core plus n-step constant sections of different heights and widths. The method is based on assuming rational-function forms for n + 1 functions R j (s) defined in Laplace space, the coefficients being constrained by physical consistency conditions. Two approximations were considered: one in which some unknown constants are fixed at their zero density value (RFA1) and one which enforces the continuity of the derivative of the cavity function (RFA2). In the former case, one has in the end to solve (numerically) n coupled transcendental equations while in the latter, which is internally more consistent, the set to be solved is made of 2n coupled transcendental equations.
The illustrative cases of the SS+SW and sSW potentials that we examined indicate that our approach, being semi-analytical in nature, is a reasonable compromise between simplicity and accuracy. For the SS+SW potential, the RFA2 approximation gives the best results, whereas the RFA1 is superior in the sSW potential. On the other hand, in these illustrative cases the present development also represents a clear improvement over the solutions of the corresponding PY integral equations, which require harder numerical labor.
These further examples of the adequacy of the RFA method as compared with the results of the PY equation reinforce the notion that such methodology is a valuable alternative to the integral equation approach for the derivation of the structural properties of fluid systems. This conclusion is of course based on a limited analysis. One could easily argue that the comparison with the results stemming out of the Ornstein-Zernike equation with either the Rogers-Young or the hypernetted chain closures would be a more solid test ground. However, these approximations require much harder numerical work, while our approach yields the explicit s-dependence of the Laplace transform G(s). In any case, we find it interesting
For simplicity, henceforth we will omit the arguments of the barred and unbarred quantities. From equation (43),

23602-9
where use has been made of equations (53)- (57). This ensures that equations (28)- (31) are satisfied by the one-step system, provided they are satisfied by the two-step system in the limit ǫ 1 → ǫ 2 . Next, the transcendental equations (44) and (47) are trivially satisfied by equations (53) and (56). Finally, the transcendental equations (45) and (48) become the same for the two-and one-step cases.